# DFT (FFT) of a Real Even Function Doesn't Yield Real Only DFT Signal

There are a lot of similar questions already, but I can't find one answering this problem.

I used the following MATLAB code to perform the FFT on a Lorentzian input signal

x = linspace(-10,10,2^10);
y = 1./(x.^2 + 1);
F = fft(y);
IF = imag(F).^2;
power = trapz(IF); %integrating the imaginary part squared


Now, from everything I've read, the FFT of a real, even function is supposed to give $0$ imaginary part. But in this case, I get

power = 3.86


which doesn't seem to be just a small numerical inaccuracy. What's going on here?

Attached is the imaginary spectrum after fftshift. • you have to understand that the points $x$ and $x[N/2]$ (where $N$ is the size of the DFT) have no reflected points. (and remember that silly MATLAB adds 1 to each index.) Jun 4 '18 at 6:01

The issue is that the DFT matrix is not actually symmetric (hence the FFT is not symmetric). This can be adjusted if you need the symmetry. Modifying your code a bit, here is my session where I use a diagonal matrix D to shift all rows of the DFT matrix to make them explicitly symmetric:

x = linspace(-10,10,2^10);
y = 1./(x.^2 + 1);
D = diag(exp(-1j*2*pi*(0:1023)/1024/2));
F = D*fft(y).';
imag(F)'*imag(F) / (F'*F)
ans =    2.8350e-32
IF = imag(F).^2;
power = trapz(IF)
power =    2.3315e-27


To visualize this, consider the DFT matrix that has a first column of all ones. The last column is $not$ all ones; hence it is not symmetric. By shifting each row by a half step, we get symmetry. The resulting symmetric DFT matrix is still unitary (I think I got this right but please double-check me):

V = fliplr(vander(exp(1j*2*pi*(0:(N-1))/N)))';
Vsym = 1/sqrt(N)*diag(sqrt(V(:,2)))*V;
norm(Vsym*Vsym'-eye(1024))
ans =    2.5503e-13
norm(Vsym'*Vsym-eye(1024))
ans =    2.5508e-13


An alternate solution is to shift your data rather than add a linear phase shift in frequency (i.e. a time shift), meaning that your data needs to not be strictly symmetric into the FFT.

• Ah, that is something peculiar I didn't ​realize. Thanks! Jun 4 '18 at 5:50
• @HiddenBabel, there's nothing really peculiar. remember that the FFT is a fast method of evaluating the DFT. so mathematically, you're dealing with the DFT and all of the properties and theorems of the DFT apply here. the DFT maps a periodic sequence of length $N$ to another periodic sequence of the same length. so $$x[n+N] = x[n] \quad \forall n \in \mathbb{Z}$$ and same for $$X[k+N] = X[k] \quad \forall k \in \mathbb{Z}$$ if $N$ is odd, everything is symmetrical $X[-k]=X[k]$ for $|k|<N/2$. but if $N$ is even, then $X[-N/2]=X[N/2]$, but they are not two reflected points. Jun 4 '18 at 16:07
• @robert bristow-johnson I was thinking about the DFT wrongly. I thought that simply distributing points symmetrically around 0 was good enough, but that was wrong. Thanks! Jun 4 '18 at 16:46
• it's not wrong. if you come away with this thinking that's wrong, your understanding of the DFT is diminished. Jun 4 '18 at 17:00
• But that is exactly what I did. The points linspace(-10,10,1024) are symmetrically placed about x=0. However, that doesn't work as we've seen. Jun 4 '18 at 20:34

When you work with DFT you need to remember one big assumption of the DFT - The samples are part of periodic signal.

For example, assume your signal in on the grid vX = [-3:3].
Let's say it is the simplest symmetric function: vY = vX .^ 2.
If you apply fft on it the result won't be real.

Why is that?
Because the DFT assumes the input is periodic.
Hence the input signal is something like [vY, vY, vY, ..., vY].

Moreover, for the DFT the indexing is [0, 1, ..., N - 1] which in this case is [0:6]. Now if we do periodic continuation, what will you get on the grid of [-7:6]?

figure();
plot([-7:6], [vY, vY]);


It won't be a be symmetric around 0 (Of course look at [-6, 6]).

Let's define vYY = vY(1:6). If you apply fft(vYY) the result is pure real.
Why is that?
Because for the DFT it sits on the grid [0:5].
If you plot a periodic continuation of it:

figure();
plot([-6:5], [vYY, vYY]);


Look on [-5, 5], now it is symmetric.
Hence its DFT is pure real.

## MATLAB Code

% The Grid
vX = [-3:3];
vY = vX .^ 2; %<! Even Function on the Grid

figure();
plot(vX, vY);
title('Figure 001');

vYDft = fft(vY);
norm(imag(vYDft), 'inf')

norm(imag(fft(vY(1:6))), 'inf')

figure();
plot([-7:6], [vY, vY])
title('Figure 002');

figure();
plot([-6:5], [vY(1:6), vY(1:6)]);
title('Figure 003');


Output:

ans = 7.2678
ans = 0

Figure 001 Figure 002 Figure 003 Your x axis goes from -10 to 10 with an even number of samples. Thus your FFT is probably not symmetric around the 1st sample or N/2 sample in the input vector, which is required for the imaginary part of an FFT result to be (within rounding error of) zero.

Symmetric around the middle only works for odd length FFTs. (try 2^10 + 1)

Try shifting your function half a sample to the left.

• power = 3.8642 (even 2^10) power = 3.8642 (odd 2^10+1) first thing I tried
– user28715
Jun 4 '18 at 4:24
• hotpaw hits it straight up. Jun 4 '18 at 4:39

try

x = linspace(10*(-1024/1024), 10*(1023/1024), 2048);
y = 1./(x.^2 + 1);
y = fftshift(y);
Y = fft(y);
max(abs(imag(Y)))/mean(abs(Y))


i didn't run this, but i bet the number that appears at bottom is very very tiny. like around 1e-15 .

• Let's not change out metric, power = 0.9666, better , but still higher than I expected
– user28715
Jun 4 '18 at 5:05
• i don't understand what you mean, Stan. Jun 4 '18 at 5:09
• whoops, have to modify this a teeny bit Jun 4 '18 at 5:12
• we should stick to the same error metric, power = trapz(IF); %integrating the imaginary part squared
– user28715
Jun 4 '18 at 5:25
• I'm not sure why this was tagged as not helpful. The point here is that the input waveform is slightly asymmetric to match the asymmetry of the fft. All of the points other than index 0 and N/2 are paired up symmetrically so that the imaginary part goes to zero (or very close to it). I believe this is the lowest cost method to get to zero imaginary.
– user35336
Jun 4 '18 at 7:22