ifft convert frequency domain to time domain
Show older comments
Hello everyone,
I have a function X in the frequency domain that I would like to convert to the time domain and get a plot in the time domain. However, the time domain singal has half of the oscillation period it should have. Would really appreciate any help! Here is the code i have:
m=1;
k=1;
c=0.01;
f=0:0.001:3;
omega=2.*f*pi;
X=1./(-omega.^2*m+1i*c.*omega+k);
% The Nyquist frequency is twice the maximum frequency in your signal
Fs = 2 * max(f) ;
N = length(f);
T = N / Fs; % Total time duration
t = 0:1/Fs:(N-1)/Fs;
% Compute the inverse Fourier transform
x = ifft(X, N);
% Plot the time-domain signal
figure
plot(t, x, 'b');
2 Comments
Star Strider
on 13 Feb 2024
If you have a one-sided complex Fourier transform (the magnitude and phase must be converted into the real and imaginary parts first), and the associated frequency vector, you will need to duplicate that as the complex conjugate, flip the complex conjugate, and concatenate it to the end of the original complex vector in order to do the inverse Fourier transform. Then, it should have the same length as the original signal.
The essential idea (assuming a row vector for simplicity) is:
reconstructed_fft = [original_complex_vector flip(conj(original_complex_vector))]
time_domain_signal = ifft(reconstructed_complex_vector)
In practice, it may be a bit more complicated than that.
The time vector then extends from 0 to length(time_domain_signal)-1 in steps of 1/(2*Fn) where ‘Fn’ is the Nyquist frequency. That should get you started.
Angie
on 13 Feb 2024
Accepted Answer
Paul
on 12 Feb 2024
1 vote
Hi Angie,
The input to ifft should cover all frequencies from 0 to Fs - df, where df is the frequency step. In this case, max(f) is the Nyquist frequency, so X is only half the points needed for the ifft. You'll have to figure out way to populate the second half of X, which you may be able to do based on what you can assume about the time domain sequence, x.
5 Comments
Angie
on 13 Feb 2024
Paul
on 14 Feb 2024
That's good. Upon further reflection, you should be able to get the result without actually populating the second half of X. If you want to pursue further, feel free to post your updated code with your solution and we can take a look.
Angie
on 14 Feb 2024
Hi Angie,
Let's take a look at your code
m=1;
k=1;
c=0.01;
f=0:0.001:5;
omega=2.*f*pi;
X=1./(-omega.^2*m+1i*c.*omega+k);
reconstructed_fft = [X flip(conj(X))];
time_domain_signal = ifft(reconstructed_fft);
Fn = 2 * max(f) ; % this suggests that max(f) is the Nyquist frequency
N = length(time_domain_signal);
t=0:1/(Fn): (N-1)/(Fn) ;
figure
% changed this line to plot the real and imaginary parts
plot(t,real(time_domain_signal),t,imag(time_domain_signal))
Sure enough, the imaginary part is relatively large, i.e., not just due to numerical roundoff error.
The problem is with this line, which is incorrect.
reconstructed_fft = [X flip(conj(X))];
To see why, let's consider a simpler problem. Suppose we have a 10-point sequence
rng(100);
z = rand(1,10)
and the sampling frequency, Fs, is 1Hz
Fs = 1;
f = (0:9)/10*Fs; % frequency vector
The fft of z is
Z = fft(z);
[f;Z].' % tanspose here for easier viewing
If you look carefully, you'll see that
a) the last four elements of Z have a nice relationship with elements 2-5 of Z.
b) element 6 (=10/2 + 1) of f is the Nyquist frequency (Fs/2).
Now, in your problem, you have the portion of the f vector from 0 to Fs/2, or
fangie = f(1:6)
and you have the corresponding elements of Z
Zangie = Z(1:6)
Given Zangie, can you reconstruct Z from which you should be able to reconstruct z using ifft(Z)?
If you can get this to work for this simple example, you should be able to get it to work for your actual problem.
A key assumption here is that max(f) is the Nyquist frequency, which implies that the original signal has an even number of elements. If you think that the original signal has an odd number of elements then we'll need to make a slight modification, though if you understand how the even case works it should be easy to figure out how the odd case works with a simple example as done above.
Keep in mind that, even if done correctly, the imaginary part of time_domain_signal might still be small and non-zero due to numerical rounding errors. Feel free to post back if you get this far and want to address the small imaginary part.
Although the solution to this problem can be found by identifying the governing pattern in simple examples, it's always better to understand the underlying math.
Let x[n], n = 0:N-1, be an N-point sequence and X[k], k = 0:N-1, be its Discrete Fourier Transform (DFT), which is what would be computed by fft.
where
Based only on this definition, we see three important features of the DFT
1) for any two values k1 and k2, X[k1] = conj(X[k2]) if theta_k1 == -theta_k2 IF x[n] is real-valued
2) we can add or subtract an integer multiple of 2*pi to any value of theta_k and not change the corresponding value of X[k]
3) X[k] must be real if theta_k = 0 or theta_k = pi, IF x[n] is real-valued
Let's look at the case with an even value of N, say
N = 8;
and define a sequence for evaluation and its DFT.
rng(100);
x = rand(N,1);
X = fft(x);
The angles that correspond to the DFT elements are
k = (0:N-1).';
theta_k = 2*sym(pi)*k/N;
For later use, we'll need the 1-based Matlab indices of the the arrays
index = (1:N).';
Put all of the information in a table for viewing
table(k,theta_k,X,index)
Use fact (2) from above and subtract 2*pi from the values of theta_k > pi
theta_k(theta_k > sym(pi)) = theta_k(theta_k > sym(pi)) -2*sym(pi);
table(k,theta_k,X,index)
Now we can see that fact(1) is verified by comparing the values of X for values of theta_k such that theta_k1 = -theta_k2.
For N even, the values of the DFT (because x[n] is real) base only on its first N/2 values in terms of Matlab indexing are (compared side-by-side with X just be sure)
[ [X(1); X(2:N/2); X(N/2+1); flip(conj(X(2:N/2)))] , X ]
Repeat the whole process for an odd value of N
N = 9;
x = rand(N,1);
X = fft(x);
k = (0:N-1).';
theta_k = 2*sym(pi)*k/N;
index = (1:N).';
theta_k(theta_k > sym(pi)) = theta_k(theta_k > sym(pi)) -2*sym(pi);
table(k,theta_k,X,index)
The odd case is the same as the even case, EXCEPT that we don't have a value of theta_k == pi in the middle. So for N odd, we have
[ [X(1); X(2:(N+1)/2); flip(conj(X(2:(N+1)/2)))] , X ]
In the even case, the key term is X(2:N/2) and in the odd case the key term is X(2:(N+1)/2). At first glance, it looks like a general solution would have to treat these cases separately. Fortunately, the indices in these expressions can both be computed using the same general expression
%subindex = 2:ceil(N/2)
Demonstrate for both cases
N = 8;
isequal(2:ceil(N/2),2:N/2)
N = 9;
isequal(2:ceil(N/2),2:(N+1)/2)
In summary, given the "lower" portion of the DFT, Xlower, and the length, N, of the real-valued time domain sequence, the full DFT of x(n) is simply (using row vectors here)
% X = [Xlower fliplr(conj(Xlower(2:ceil(N/2)))) ]
Whether N is even or odd must be determined through other means.
In this particular problem, @Angie specified that the final point in Xlower (which was called X in the Question) corresponded to the Nyquist frequency, which in turn corresponds to theta_k = pi. Hence, we can deduce for this question that N is even.
Let's apply to the original example to explore some practical effects.
Define the underlying continuous-time signal based on its Continuous Time Fourier Transform (CTFT)
m = 1;
k = 1;
c = 0.01;
syms X(w) x(t)
assume([w t],'real');
X(w) = 1/(-w.^2*m + 1i*c.*w + k);
x(t) = simplify(ifourier(X(w),w,t))
simplify(imag(x(t)))
x(t) = simplify(real(x(t)),100)
The lower portion of the frequency vector is:
f_lower = 0:0.001:5; % assuming f_lower(end) = Fs/2
f we assume that f_lower(end) is the Nyquist frequency, which corresponds to theta_k = pi, we have
Fs = f_lower(end)*2;
Ts = 1/Fs;I
N = 2*(numel(f_lower)-1)
omega = 2.*f_lower*pi;
Xlower = 1./(-omega.^2*m + 1i*c.*omega + k);
Note that the elements of Xlower are samples of the CTFT of x(t), which are very close to, but not exactly, the (scaled) DFT of samples of x(t).
Now we apply the recipe
Xrecon = [Xlower , fliplr(conj(Xlower(2:ceil(N/2)))) ];
xrecon = ifft(Xrecon);
Because we started by taking samples of the CTFT, we have to scale the ifft result by Ts if we want to compare xrecon to x(t)
xrecon = xrecon/Ts;
Compare the real part of xrecon to x(t)
tval = (0:N-1)/Fs;
figure
fplot(x(t),[0 tval(end)]);
hold on
plot(tval,real(xrecon));
ylim([-1 1])
xlim([0 200])
xlabel('Time (sec)');
ylabel('real(x)')
That looks like a good match.
Now check the imaginary part of xrecon.
figure
plot(tval,imag(xrecon))
xlabel('Time (sec)');
ylabel('Imag(x)')
It should be exactly zero because we applied the flip(conj()) recipe, but it isn't. The reason is fact (3) from the top of this answer. If Xlower is the lower part of the DFT of a real-valued, even-length sequence, then the last element of Xlower, which corresponds to theta_k = pi, must be real. But it isn't,
format short e
Xlower(end)
which is a consequence of how Xlower was formed by sampling the CTFT of x(t).
We can address this issue by forcing the imaginary part of xrecon to zero
%xrecon = real(xrecon)
or perhaps forcing that last element of Xlower to be real
% Xlower(end) = real(Xlower(end));
% etc.
or telling ifft that we intended Xrecon to satisfy the DFT criteria for a real-valued signal by using the symmetric flag
xrecon = ifft(Xrecon,'symmetric');
xrecon = xrecon/Ts;
Now ifft forces the result to be real
any(imag(xrecon))
and we still match the original signal
figure
fplot(x(t),[0 tval(end)]);
hold on
plot(tval,real(xrecon));
ylim([-1 1])
xlim([0 200])
xlabel('Time (sec)');
ylabel('real(x)')
Apparently ifft doesn't even use the upper portion of the input if the symmetric flag is set. So we don't need to fill out the upper portion of Xrecon with its actual values. Apparently any values will do, particularly zeros
Xrecon = [Xlower , zeros(1,ceil(N/2)-1)];
xrecon = ifft(Xrecon,'symmetric')/Ts;
any(imag(xrecon))
figure
fplot(x(t),[0 tval(end)]);
hold on
plot(tval,real(xrecon));
ylim([-1 1])
xlim([0 200])
xlabel('Time (sec)');
ylabel('real(x)')
More Answers (1)
Basit
on 5 Oct 2024
0 votes
I think we can not get the time domain data as good signal by using iFFT
Categories
Find more on Spectral Measurements in Help Center and File Exchange
