ifft convert frequency domain to time domain

49 visualizzazioni (ultimi 30 giorni)
Angie
Angie il 12 Feb 2024
Commentato: Paul il 18 Feb 2024
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 Commenti
Star Strider
Star Strider il 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.

Accedi per commentare.

Risposta accettata

Paul
Paul il 12 Feb 2024
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 Commenti
Paul
Paul il 14 Feb 2024
Modificato: Paul il 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)
z = 1×10
0.5434 0.2784 0.4245 0.8448 0.0047 0.1216 0.6707 0.8259 0.1367 0.5751
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
ans =
0.0000 + 0.0000i 4.4258 + 0.0000i 0.1000 + 0.0000i 0.2230 + 0.2742i 0.2000 + 0.0000i -0.6682 - 0.5093i 0.3000 + 0.0000i 1.2644 + 1.0959i 0.4000 + 0.0000i 0.1177 + 0.0387i 0.5000 + 0.0000i -0.8656 + 0.0000i 0.6000 + 0.0000i 0.1177 - 0.0387i 0.7000 + 0.0000i 1.2644 - 1.0959i 0.8000 + 0.0000i -0.6682 + 0.5093i 0.9000 + 0.0000i 0.2230 - 0.2742i
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)
fangie = 1×6
0 0.1000 0.2000 0.3000 0.4000 0.5000
and you have the corresponding elements of Z
Zangie = Z(1:6)
Zangie =
4.4258 + 0.0000i 0.2230 + 0.2742i -0.6682 - 0.5093i 1.2644 + 1.0959i 0.1177 + 0.0387i -0.8656 + 0.0000i
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.
Paul
Paul il 18 Feb 2024
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.
The DFT is defined as a sum of complex exponentials weighted by the elements of the input sequence
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)
ans = 8×4 table
k theta_k X index _ ________ _________________ _____ 0 0 3.714+0i 1 1 pi/4 0.63618+0.12198i 2 2 pi/2 -0.54714+1.2707i 3 3 (3*pi)/4 0.44119-0.37049i 4 4 pi -0.42718+0i 5 5 (5*pi)/4 0.44119+0.37049i 6 6 (3*pi)/2 -0.54714-1.2707i 7 7 (7*pi)/4 0.63618-0.12198i 8
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)
ans = 8×4 table
k theta_k X index _ _________ _________________ _____ 0 0 3.714+0i 1 1 pi/4 0.63618+0.12198i 2 2 pi/2 -0.54714+1.2707i 3 3 (3*pi)/4 0.44119-0.37049i 4 4 pi -0.42718+0i 5 5 -(3*pi)/4 0.44119+0.37049i 6 6 -pi/2 -0.54714-1.2707i 7 7 -pi/4 0.63618-0.12198i 8
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 ]
ans =
3.7140 + 0.0000i 3.7140 + 0.0000i 0.6362 + 0.1220i 0.6362 + 0.1220i -0.5471 + 1.2707i -0.5471 + 1.2707i 0.4412 - 0.3705i 0.4412 - 0.3705i -0.4272 + 0.0000i -0.4272 + 0.0000i 0.4412 + 0.3705i 0.4412 + 0.3705i -0.5471 - 1.2707i -0.5471 - 1.2707i 0.6362 - 0.1220i 0.6362 - 0.1220i
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)
ans = 9×4 table
k theta_k X index _ _________ _________________ _____ 0 0 4.116+0i 1 1 (2*pi)/9 1.0333+0.22082i 2 2 (4*pi)/9 -1.3691+0.30323i 3 3 (2*pi)/3 -1.2096+0.062645i 4 4 (8*pi)/9 0.10258+0.10967i 5 5 -(8*pi)/9 0.10258-0.10967i 6 6 -(2*pi)/3 -1.2096-0.062645i 7 7 -(4*pi)/9 -1.3691-0.30323i 8 8 -(2*pi)/9 1.0333-0.22082i 9
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 ]
ans =
4.1160 + 0.0000i 4.1160 + 0.0000i 1.0333 + 0.2208i 1.0333 + 0.2208i -1.3691 + 0.3032i -1.3691 + 0.3032i -1.2096 + 0.0626i -1.2096 + 0.0626i 0.1026 + 0.1097i 0.1026 + 0.1097i 0.1026 - 0.1097i 0.1026 - 0.1097i -1.2096 - 0.0626i -1.2096 - 0.0626i -1.3691 - 0.3032i -1.3691 - 0.3032i 1.0333 - 0.2208i 1.0333 - 0.2208i
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)
ans = logical
1
N = 9;
isequal(2:ceil(N/2),2:(N+1)/2)
ans = logical
1
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))
x(t) = 
simplify(imag(x(t)))
ans = 
0
x(t) = simplify(real(x(t)),100)
x(t) = 
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)
N = 10000
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)
ans =
-1.0142e-03 - 3.2317e-07i
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))
ans = logical
0
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))
ans = logical
0
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)')

Accedi per commentare.

Più risposte (1)

Star Strider
Star Strider il 14 Feb 2024
Modificato: Star Strider il 14 Feb 2024
The ‘reconstructed_fft’ vector is symmetric because of the way you constructed it (the default assumption that fft makes is that it is asymmetric), so you need to inform fft of that fact. Then the imaginary parts of the inversion disappear.
Try this —
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);
SzX = size(X)
SzX = 1×2
1 5001
reconstructed_fft = [X flip(conj(X))];
time_domain_signal = ifft(reconstructed_fft,'symmetric');
[a,b] = bounds(imag(time_domain_signal))
a = 0
b = 0
Fn = 2 * max(f) ;
N = length(time_domain_signal);
t=0:1/(Fn): (N-1)/(Fn) ;
figure
plot(t,time_domain_signal)
xlim([0 25])
.
  5 Commenti
Star Strider
Star Strider il 14 Feb 2024
I copied Angie’s posted code in this Comment, then ran it to understand the problem.
The only changes I made were to return the size of ‘X’ and to add 'symmetric' to the fft call.
Keeping the value of ‘c’ as in the original produces a pure sine curve. (Changing it temporarily was an experiment. I changed it back just now.)
As for ‘reconstructed_fft’ ask Angie! That approach was in my original Comment.
Paul
Paul il 14 Feb 2024
Yes, I know that approach was in your original comment. As I've said, I don't believe it's correct. Again, let's try a simple example. Because reconstructed_fft has an even length, so must the original time domain sequence. Here's an original time domain signal with an even length and its FFT
rng(100);
xoriginal = rand(1,8);
Xoriginal = fft(xoriginal);
The X in the problem would be the first half
X = Xoriginal(1:4);
Now follow the recipe (w/o symmetric):
reconstructed_fft = [X flip(conj(X))];
time_domain_signal = ifft(reconstructed_fft);
[xoriginal; time_domain_signal].'
ans =
0.5434 + 0.0000i 1.0610 + 0.0000i 0.2784 + 0.0000i 0.5385 - 0.2231i 0.4245 + 0.0000i 0.6055 - 0.6055i 0.8448 + 0.0000i 0.0990 - 0.2389i 0.0047 + 0.0000i 0.0000 + 0.3798i 0.1216 + 0.0000i -0.0685 - 0.1654i 0.6707 + 0.0000i 0.4110 + 0.4110i 0.8259 + 0.0000i 1.0674 + 0.4421i
These are not the same.
Adding the symmetric directive gets rid of the imaginary part, but does not recover the original signal
reconstructed_fft = [X flip(conj(X))];
time_domain_signal = ifft(reconstructed_fft,'symmetric');
[xoriginal; time_domain_signal].'
ans = 8×2
0.5434 0.6520 0.2784 0.1698 0.4245 0.5331 0.8448 0.7362 0.0047 0.1133 0.1216 0.0130 0.6707 0.7793 0.8259 0.7173

Accedi per commentare.

Prodotti


Release

R2020b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by