when do I need use fftshift?
显示 更早的评论
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!
采纳的回答
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 个评论
Daniel Niu
2023-12-20
I'm not 100% sure what you're asking. I think you're generally asking about how to use fft and ifft to analyze continuous-time signals after sampling.
Maybe it's easiest to illustrate with an example:
Assume we have a continuous-time signal that is a triangular pulse that is symmetric around the origin
syms t real
x(t) = triangularPulse(t);
figure
fplot(x(t)),axis padded
Its Continuous-Time Fourier Transform (CTFT) is
syms w real
X(w) = fourier(x(t),t,w);
X(w) is always real and positive, so in plots below we won't bother with phase.
angle(simplify(expand(rewrite(X(w),'sincos'))))
And we need to deal with a situation when w = 0
X(w) = piecewise(w == 0, limit(X(w),w,0),X(w)); % deal with special case at w = 0
Now we sample the signal so that we can use fft. Recall that the DFT (implemented by fft) always assumes that the input is one period of a periodic, discrete-time signal, even though our continuous-time signal is not periodic. So we have to "trick" the DFT by taking enough zero-valued samples off the ends so that the DFT thinks the periodic extension of the sampled signal has a very large period (as the period gets larger we get closer and closer to the continuous-time signal). I'll assume that taking samples over -10 < t < 10 (20 second period) is sufficient. And we need to select a small sampling period that also provides equal spacing over that interval
Ts = 0.1;
tval = -10:Ts:10; N = numel(tval);
Sample the signal
xval = double(x(tval));
As shown previously, because we've taken samples symmetrically around t =0, we use ifftshift before calling fft (also have to multiply by Ts for later comparison with the CTFT):
Xdft = Ts*fft(ifftshift(xval));
The frequency points corresponding to the DFT samples are (rad/sec)
wdft = (0:N-1)/N*2*pi/Ts;
Plot the DFT
figure
stem(wdft,abs(Xdft)),title('DFT')
In discrete-time, the Fourier transform is 2pi-periodic in angular frequency (i.e., rad/sec*Ts) and Matlab's convention with fft is to output the period from 0 to 2*pi.
To compare back to the CTFT we want the "central" period, which can be obtained with fftshift
Xdftshift = fftshift(Xdft);
The corresponding frequency vector is now
wdftshift = ceil(-N/2:(N-1)/2)/N*2*pi/Ts;
Compare to the CTFT
figure
hold on
fplot(abs(X(w)),[wdftshift(1) wdftshift(end)])
plot(wdftshift,abs(Xdftshift),'.')
axis padded
Now, going the other way, suppose we obtain frequency-domain samples from the CTFT
wctft = wdftshift;
Xctft = double(X(wctft));
figure
stem(wctft,abs(Xctft))
Before we can use ifft, we have to shift so that we get samples of what the DFT would be over 0-2pi (normalized), i.e., the plot above titled 'DFT'. So we use ifftshift and then call ifft (again, we need to scale by Ts)
xval = ifft(ifftshift(Xctft),'symmetric')/Ts;
figure
stem((0:numel(xval)-1)*Ts,xval)
Just like fft, the output of ifft can be viewed as one period of a periodic signal and ifft returns the period from 0 < t < T (T=20). So we use fftshift to get back the "central" period
xvalshift = fftshift(xval);
Compare to the original signal
figure
hold on
fplot(x(t),[-10 10])
plot(tval,xvalshift,'.')
axis padded
Very close. I'm sure if we zoom in we'd see some differences due to the sampling and windowing that yielded the frequency domain samples.
Keep in mind that this example worked because the time and frequency domain signals have a particular kind of symmetry around the origin and the sampling of those signals was also symmetric around the origin. If, for example, the triangular pulse was shifted to the right or left, the analysis would have to be done a bit differently.
Daniel Niu
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 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Multirate Signal Processing 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
