# 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

**Category:**- Fourier transforms

### Note

Comments are closed.

## 9 CommentsOldest to Newest

**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

**2**of 9

Sumedh—I used Microsoft Visio for those.

**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

**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.

**5**of 9

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

**6**of 9

K—I used Visio.

**7**of 9

Thanks Steve!

K

**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

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

## Recent Comments