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!

Answer 1

[Paul]

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.

[Daniel Niu]

Dear Paul,

Thank you for your elaborate answers. I also want to know what the default input and output of ifft?

thank you!

[Paul]

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

ans = 

0

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 sampled, discrete-time signal is periodic, 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 low period (as the period gets lower, i.e., more samples, 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]

Dear @Paul,

would you please see my another question?

why the two split-step fourier methods using fft and ifft is equivalent ? - MATLAB Answers - MATLAB Central (mathworks.com)

according to your answers, the fft 's default input is t=[0 T], what happen if I give the [-T/2 T/2] to fft function?

thank you

[Paul]

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.