when do I need use fftshift?
108 ビュー (過去 30 日間)
古いコメントを表示
Bx = 10;
A = sqrt(log(2))/(2*pi*Bx);
fs = 500; %sampling frequency
dt = 1/fs; %time step
T=1; %total time window
t = -T/2:dt:T/2-dt; %time grids
df = 1/T; %freq step
Fmax = 1/2/dt; %freq window
f=-Fmax:df:Fmax-df; %freq grids, not used in our examples, could be used by plot(f, X)
%-------------------------------------------------------------------
% Numerical evaluations
x = exp(-t.^2/2/A^2);
Xan = A*sqrt(2*pi)*exp(-2*pi^2*f.^2*A^2); %X(f), analytical Fourier transform of x(t), real
Xfft = dt * fft(x); %directly using fft()
Xfftshift = dt * fft(fftshift(x)); %using fftshift() before fft()
Xfinal = dt * fftshift(fft(fftshift(x)));
figure;
plot(f, imag(Xfinal), 'LineWidth', 3, 'DisplayName', 'Imaginary (Numerical)');
hold on;
plot(f, real(Xfinal), 'LineWidth', 3, 'DisplayName', 'Real (Numerical)');
plot(f, Xan, '--', 'LineWidth', 2, 'DisplayName', 'Analytical');
legend('show');
xlabel('Frequency (Hz)');
ylabel('Amplitude');
title('Fourier Transform of Gaussian Pulse');
grid on;
%shift property of Fourier transform, gaussian function move right T/2
%should be a complex number. after test, the x should not use fftshift can
%get the right result.
%different from whyusefftshift.m is t and x
clear
Bx = 10;
A = sqrt(log(2))/(2*pi*Bx);
fs = 500; %sampling frequency
dt = 1/fs; %time step
T = 1; %total time window
t = 0:dt:T-dt; %time grids starting from t=0
df = 1/T; %freq step
Fmax = 1/2/dt; %freq window
f = -Fmax:df:Fmax-df; %freq grids, not used in our examples, could be used by plot(f, X)
% Numerical evaluations
x = exp(-(t - T/2).^2 / (2*A^2)); % Adjusted the time grid
Xan_real = A * sqrt(2*pi) * exp(-2*pi^2*f.^2*A^2).*cos(pi*T*f); %X(f), analytical Fourier transform of x(t), real
Xfft = dt * fft(x); %directly using fft()
%Xfinal = dt * fft(fftshift(x)); %using fftshift() before fft()
%Xfinal = dt * fftshift(fft(fftshift(x)));
%Xfinal = dt * fftshift(fft(fftshift(x)));
Xfinal = dt * fftshift(fft(x));
figure;
plot(f, imag(Xfinal), 'LineWidth', 10, 'DisplayName', 'Imaginary (Numerical)');
legend('show');
figure;
plot(f, imag(Xfinal), 'LineWidth', 10, 'DisplayName', 'Imaginary (Numerical)');
hold on;
plot(f, real(Xfinal), 'LineWidth', 3, 'DisplayName', 'Real (Numerical)');
plot(f, Xan_real, '--', 'LineWidth', 2, 'DisplayName', 'Analytical-real');
legend('show');
xlabel('Frequency (Hz)');
ylabel('Amplitude');
title('Fourier Transform of Gaussian Pulse');
grid on;
In the above two code snippet, I can figure out when I need use the fftshift.
In the first code, I must using fftshift before using fft function to get the right result. But In the second code,
I must not using fftshift before x and must using fft directly to get the right result.
Your help would be highly appreciated!
0 件のコメント
採用された回答
Paul
2023 年 12 月 18 日
編集済み: Paul
2023 年 12 月 20 日
Hi Daniel,
Let's plot the first function
Bx = 10;
A = sqrt(log(2))/(2*pi*Bx);
fs = 500; %sampling frequency
dt = 1/fs; %time step
T=1; %total time window
t = -T/2:dt:T/2-dt; %time grids
x = exp(-t.^2/2/A^2);
plot(t,x)
Two things that are very important to remember about the Discrete Fourier Transfrom (DFT), as is computed by fft.
One: there is an underlying assumption that the time-domain sequence is periodic.
Two: the fft implemenation of the DFT operates on one period of the sequence starting at t = 0.
The periodic sequence in item one is the so-called periodic extension of x, with period equal to the time duration of x. We can graph a couple of periods of this periodic extension by
figure
plot([t, t+T , t + 2*T],[x x x])
Now, in order for fft to work the way we want based on point 2, we need to provide it the elements of this periodic extension between 0 and 1. Because of the symmetry in the time vector around t = 0 by construction, we can use ifftshift (ifftshift should be used because numel(t) is even and the first point is at -T/2 and the last at T/2-dt; for aperiodic symmetric signals, it's better (in principle) to define the time vector with an odd number of points symmetric around the point of symmetry, which is t = 0 in this case)
figure
plot((0:numel(t)-1)*dt,ifftshift(x)) % note the last point is at T-dt
We see that ifftshift produces a single period of the periodic extension of the signal over the interval starting at t = 0
Let's check the second case
clear
Bx = 10;
A = sqrt(log(2))/(2*pi*Bx);
fs = 500; %sampling frequency
dt = 1/fs; %time step
T = 1; %total time window
t = 0:dt:T-dt; %time grids starting from t=0
% Numerical evaluations
x = exp(-(t - T/2).^2 / (2*A^2)); % Adjusted the time grid
figure
plot(t,x)
Here, x is defined over the desired interval, so there is no need for additional shifting.
Finally, it's important to remember that bolded clause above. ifftshift (and fftshift) assumes a particular symmetry (depending on if the sequence is even or odd length). If that symmetry is not present in the original time vector, then ifftshift (and fftshift) won't give the correct result. However, in the general case, we can use circshift to map an input sequence defined over any interval to the interval expected by fft.
4 件のコメント
Paul
2023 年 12 月 21 日
編集済み: Paul
2023 年 12 月 21 日
I wouldn't say that's the default input as that implies that there's some way to specify something other than a default, but there isn't. And I think it would be from 0 <= t <= T - dt.
If the first element in the input to fft does not correspond to t = 0, then the result will have the correct magntiude but the phase will be offset by a linear phase factor.
For example, suppose we have an 11 point sequence that defined for n = -5:5;
rng(100)
x = rand(1,11);
n = -5:5;
Assuming we want to use fft to get frequency domain samples of one period of the Discrete Time Fourier Transform of x, we'd do
X1 = fft(ifftshift(x));
w = (0:10)/11*2*pi; % rad
If we just take the fft
X2 = fft(x);
The magnitudes are the same
figure
hold on
stem(w,abs(X1))
stem(w,abs(X2))
But the phases are different
figure
hold on
stem(w,angle(X1))
stem(w,angle(X2))
The difference in phase is a lnear phase term determined by the offset of the first point in x from the origin
figure
hold on
stem(w,angle(X1))
stem(w,angle(X2.*exp(-1j*w*n(1))))
If magnitude is all that's needed, then no need to worry about the shifting.
その他の回答 (0 件)
参考
カテゴリ
Help Center および File Exchange で Spectral Measurements についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!