ifft convert frequency domain to time domain

Question

Angie 2024-2-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2081373-ifft-convert-frequency-domain-to-time-domain

评论： Paul 2024-2-18

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 个评论
显示无隐藏无

Star Strider 2024-2-13

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 2024-2-13

Thank you @Star Strider! it is clear now

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Paul 2024-2-12

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2081373-ifft-convert-frequency-domain-to-time-domain#answer_1407378

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 个评论
显示 3更早的评论隐藏 3更早的评论

Paul 2024-2-14

编辑：Paul 2024-2-14

在 MATLAB Online 中打开

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 = 
0000 + 0.0000i   4.4258 + 0.0000i
1000 + 0.0000i   0.2230 + 0.2742i
2000 + 0.0000i  -0.6682 - 0.5093i
3000 + 0.0000i   1.2644 + 1.0959i
4000 + 0.0000i   0.1177 + 0.0387i
5000 + 0.0000i  -0.8656 + 0.0000i
6000 + 0.0000i   0.1177 - 0.0387i
7000 + 0.0000i   1.2644 - 1.0959i
8000 + 0.0000i  -0.6682 + 0.5093i
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 2024-2-18

在 MATLAB Online 中打开

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)')

请先登录，再进行评论。

Answer 2

Star Strider 2024-2-14

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2081373-ifft-convert-frequency-domain-to-time-domain#answer_1408668

编辑：Star Strider 2024-2-14

在 MATLAB Online 中打开

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 个评论
显示 3更早的评论隐藏 3更早的评论

Star Strider 2024-2-14

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 2024-2-14

在 MATLAB Online 中打开

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