Steve on Image Processing

Concepts, algorithms & MATLAB

Complex surprises from fft 9

Posted by Steve Eddins,

One of the discrete-time Fourier transform properties that many people learn is that if a sequence is conjugate symmetric, , then the Fourier transform is real.

Therefore it surprises people sometimes when the output of fft is unexpectedly complex. For example:

x = [1 2 3 2 1]
x =

     1     2     3     2     1

The sequence above appears to be symmetric, but the output of fft(x) is complex:

fft(x)
ans =

   9.0000            -2.1180 - 1.5388i   0.1180 + 0.3633i   0.1180 - 0.3633i  -2.1180 + 1.5388i

So what's going on? Well, the fft function is computing the discrete Fourier transform of a sequence that is nonzero over the interval . Here's a plot:

You can see that in fact isn't actually symmetric about the origin: . is a shifted version of a symmetric sequence, and that shift causes the output of fft to be complex.

This is the symmetric sequence that has a real-valued Fourier transform:

So how do we compute the real-valued Fourier transform of that sequence using the fft function? If you'll allow me to be a bit hand-wavy here, I'll just say that the discrete Fourier transform has an implicit periodicity in both the time domain and in the frequency domain. Related to this periodicity, many of the discrete Fourier transform analogs to the more familiar Fourier transform properties involve circular shifts, or modulo indexing.

For example, the conjugate symmetry property for the discrete Fourier transform looks like this: if , then the DFT is real. What we want to do, then, is circularly shift x so that the center element moves to the left of the vector, like this:

xs = circshift(x,[0 -2])
xs =

     3     2     1     1     2

Now take the DFT of xs:

fft(xs)
ans =

    9.0000    2.6180    0.3820    0.3820    2.6180

Now the output is real, as expected.

If you're going to zero-pad, do that first and then apply the circular shift.

x128 = x;
x128(128) = 0;
x128s = circshift(x128,[0 -2]);
X128 = fft(x128s);
isreal(X128)
ans =

     0

Whoops. Was I wrong about using this procedure to produce a real result? Not really, it's just our old friend floating-point round-off error. Look how small the imaginary part is:

max(abs(imag(X128(:))))
ans =

  4.4409e-016

You can get rid of the negligible imaginary part by using real. Let's plot the real-valued Fourier transform. I'll use the frequency axis labeling technique I showed in my June 25 blog post.

w = unwrap(fftshift(2*pi * (0:(128-1)) / 128) - 2*pi);
plot(w/pi, fftshift(real(X128)))
xlabel('radians / \pi')


Get the MATLAB code

Published with MATLAB® 7.10

9 CommentsOldest to Newest

sumedh replied on : 1 of 9

Hi Steve,

Nice post; it’s something I’ve encountered (as recently as yesterday) many times.

I have an unrelated question though. I really like the plots you make for your posts (specifically the first two in this one). What do you use to make them?

Sumedh

Amanmeet Garg replied on : 3 of 9

Hi Steve,
I have read through this post a couple of times and have not been able to find out the solution to my problem.

I have a fft output of a time series, I add random phase to it, now I want to make it(random phased fft output) symmetric so that I get a real ifft output.

according to your post above, fftshift of the fft(signal) should make it symmetric and give me a real ifft output.

it is not so.

Please comment. Any help will be appreciated.

Regards
Aman

Steve replied on : 4 of 9

Amanmeet—I didn’t say anything at all like that in this post. If your time series is real, then the output of fft is conjugate symmetric modulo the FFT length. For an unspecified reason you then made it nonsymmetric by adding “random phase” to it. You’ll have to decide for yourself, based on your own application and reason for adding random phase, how to make it symmetric again. fftshift doesn’t have anything to do with it.

Hi Steve,

Your graphs on your post are really great. Could you tell me how to make first and second graph if you used MATLAB to make it?

K

Maeve replied on : 8 of 9

the loglog plot of the magnitude of the real and imaginary parts of my DFT [(real(DFT).^2+imag(DFT).^2).^0.5] versus frequency looks reasonable, can i trust that the energy is spread across frequency space as the plot indicates? what is the reason for forcing the data to be symmetric about zero (why not compute and plot the magnitude as i have done)? thanks for any advice you have! -maeve

Vittorio replied on : 9 of 9

Hi Steve!
Thanks for your posts, they have been really usefull!
I wanted to ask you a question about the phase of 3D rect pulse.
I have been able to create a 3D rectangular pulse and to evaluate the fft of it, but when it comes to the phase it looks like it’s wrong shifted.
I mean it looks like matlab is receiving the pulse not centered in the way it is wanted, so the phase is multiplied for a linear shifting.
Is there a way to fix it, to have a correct phase’s graph?
Thank you in advance! I hope it’s the right place to ask you questions…

npoints=256;
perc=1;
dt=6*1E-7/(npoints*perc); % Tempo di campionamento
df=1/(npoints*dt); % Frequenza di campionamento

t(1)=0;
f(1)=0;
for k=2:npoints/2
t(k)=(k-1)*dt;
t(npoints-k+2)=-t(k);
f(k)=(k-1)*df;
f(npoints-k+2)=-f(k);
end
t(npoints/2+1)=t(npoints/2+2)-dt;
f(npoints/2+1)=f(npoints/2+2)-df;

ts=ifftshift(t);
fs=fftshift(f);

figure
[X,Y]=meshgrid(ts,ts);

D = npoints/2; % to indicate origin at the center of the function
a = 100; % change it to enlarge or reduce the pulse
y = repmat(1:npoints,npoints,1);
x = y’;
rect = zeros(npoints);
rect(D-a:D+a-1,D-a:D+a-1) = ones(2*a);
rect=ifftshift(rect);
surf(X,Y,rect);
shading interp
axis tight
title (‘Rect 3D’);

R = fft2(ifftshift(rect));
R = fftshift(R);
[X,Y]=meshgrid(fs,fs);
figure;
surf(X,Y,abs(R));
shading interp
axis tight
title(‘Fourier Transform of Rectangular function’);
%plot real part
figure;
surf(X,Y,real(R));
shading interp
axis tight
title(‘Real part’);

Rm=abs(R);
imm=imag(R);
re=real(R);
re(abs(re) < 1e-12) = 0;
imm(abs(imm) < 1e-12) = 0;
phase=atan2(imm, re);

%plot phase
[X,Y]=meshgrid(fs,fs);
figure
surf(X,Y,phase);
shading flat
axis tight
title ('Spettro di Fase rect');

Add A Comment

What is 2 + 4?

Preview: hide