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

Question

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.

Answer 1

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.

Answer 2

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.

Answer 3

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

Answer 4

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 .