1
$\begingroup$

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. imaginary spectrum

$\endgroup$
  • $\begingroup$ you have to understand that the points $x[0]$ 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.) $\endgroup$ – robert bristow-johnson Jun 4 '18 at 6:01
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Ah, that is something peculiar I didn't ​realize. Thanks! $\endgroup$ – HiddenBabel Jun 4 '18 at 5:50
  • $\begingroup$ @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. $\endgroup$ – robert bristow-johnson Jun 4 '18 at 16:07
  • $\begingroup$ @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! $\endgroup$ – HiddenBabel Jun 4 '18 at 16:46
  • $\begingroup$ it's not wrong. if you come away with this thinking that's wrong, your understanding of the DFT is diminished. $\endgroup$ – robert bristow-johnson Jun 4 '18 at 17:00
  • $\begingroup$ 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. $\endgroup$ – HiddenBabel Jun 4 '18 at 20:34
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ power = 3.8642 (even 2^10) power = 3.8642 (odd 2^10+1) first thing I tried $\endgroup$ – Stanley Pawlukiewicz Jun 4 '18 at 4:24
  • $\begingroup$ hotpaw hits it straight up. $\endgroup$ – robert bristow-johnson Jun 4 '18 at 4:39
1
$\begingroup$

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

enter image description here

Figure 002

enter image description here

Figure 003

enter image description here

$\endgroup$
-1
$\begingroup$

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 .

$\endgroup$
  • $\begingroup$ Let's not change out metric, power = 0.9666, better , but still higher than I expected $\endgroup$ – Stanley Pawlukiewicz Jun 4 '18 at 5:05
  • $\begingroup$ i don't understand what you mean, Stan. $\endgroup$ – robert bristow-johnson Jun 4 '18 at 5:09
  • $\begingroup$ whoops, have to modify this a teeny bit $\endgroup$ – robert bristow-johnson Jun 4 '18 at 5:12
  • $\begingroup$ we should stick to the same error metric, power = trapz(IF); %integrating the imaginary part squared $\endgroup$ – Stanley Pawlukiewicz Jun 4 '18 at 5:25
  • 1
    $\begingroup$ 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. $\endgroup$ – user35336 Jun 4 '18 at 7:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.