{"id":336,"date":"2010-07-16T18:20:22","date_gmt":"2010-07-16T22:20:22","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/2010\/07\/16\/complex-surprises-from-fft\/"},"modified":"2019-10-29T13:31:18","modified_gmt":"2019-10-29T17:31:18","slug":"complex-surprises-from-fft","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2010\/07\/16\/complex-surprises-from-fft\/","title":{"rendered":"Complex surprises from fft"},"content":{"rendered":"
\r\n

One of the discrete-time Fourier transform properties that many people learn is that if a sequence is conjugate symmetric,\r\n , then the Fourier transform is real.\r\n <\/p>\r\n

Therefore it surprises people sometimes when the output of fft<\/tt> is unexpectedly complex. For example:\r\n <\/p>

x = [1 2 3 2 1]<\/pre>
\r\nx =\r\n\r\n     1     2     3     2     1\r\n\r\n<\/pre>
The sequence above appears to be symmetric, but the output of fft(x)<\/tt> is complex:\r\n   <\/p>
fft(x)<\/pre>
\r\nans =\r\n\r\n   9.0000            -2.1180 - 1.5388i   0.1180 + 0.3633i   0.1180 - 0.3633i  -2.1180 + 1.5388i\r\n\r\n<\/pre>
So what's going on? Well, the fft<\/tt> function is computing the discrete Fourier transform of a sequence   that is nonzero over the interval  . Here's a plot:\r\n   <\/p>\r\n   
 <\/p>\r\n   
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<\/tt> to be complex.\r\n   <\/p>\r\n   
This<\/b> is the symmetric sequence that has a real-valued Fourier transform:\r\n   <\/p>\r\n   
 <\/p>\r\n   
So how do we compute the real-valued Fourier transform of that sequence using the fft<\/tt> 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\r\n      periodicity in both the time domain and in the frequency domain.  Related to this periodicity, many of the discrete Fourier\r\n      transform analogs to the more familiar Fourier transform properties involve circular shifts, or modulo indexing.\r\n   <\/p>\r\n   
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<\/tt> so that the center element moves to the left of the vector, like this:\r\n   <\/p>
xs = circshift(x,[0 -2])<\/pre>
\r\nxs =\r\n\r\n     3     2     1     1     2\r\n\r\n<\/pre>
Now take the DFT of xs<\/tt>:\r\n   <\/p>
fft(xs)<\/pre>
\r\nans =\r\n\r\n    9.0000    2.6180    0.3820    0.3820    2.6180\r\n\r\n<\/pre>
Now the output is real, as expected.<\/p>\r\n   
If you're going to zero-pad, do that first and then apply the circular shift.<\/p>
x128 = x;\r\nx128(128) = 0;\r\nx128s = circshift(x128,[0 -2]);\r\nX128 = fft(x128s);\r\nisreal(X128)<\/pre>
\r\nans =\r\n\r\n     0\r\n\r\n<\/pre>
Whoops. Was I wrong about using this procedure to produce a real result? Not really, it's just our old friend floating-point\r\n      round-off error. Look how small the imaginary part is:\r\n   <\/p>
max(abs(imag(X128(:))))<\/pre>
\r\nans =\r\n\r\n  4.4409e-016\r\n\r\n<\/pre>
You can get rid of the negligible imaginary part by using real<\/tt>. Let's plot the real-valued Fourier transform. I'll use the frequency axis labeling technique I showed in my June 25 blog post<\/a>.\r\n   <\/p>
w = unwrap(fftshift(2*pi * (0:(128-1)) \/ 128) - 2*pi);\r\nplot(w\/pi, fftshift(real(X128)))\r\nxlabel('radians \/ \\pi'<\/span>)<\/pre>