FFT of a real symmetric vector is not real and symmetric

Question

It’s hard for me to understand what should be a simple concept. I built a vector in MATLAB, which is real and symmetrical. When I take the FFT in MATLAB, the result has a significant imaginary component, although the rules of symmetry of the Fourier transform say that the FT of a real symmetric function must also be real and symmetric. My sample code is:

N = 1 + 2^8; k = linspace(-1,1,N); V = exp(-abs(k)); Vf1 = fft(fftshift(V)); Vf2 = fft(ifftshift(V)); Vf3 = ifft(fftshift(V)); Vf4 = ifft(ifftshift(V)); Vf5 = fft(V); Vf6 = ifft(V); disp([isreal(Vf1) isreal(Vf2) isreal(Vf3) isreal(Vf4) isreal(Vf5) isreal(Vf6)]) 

Result:

 0 0 0 0 0 0 

No combination of (i)fft or (i)fftshift results in a real symmetric vector. I tried with even and odd N ( N = 2^8 vs. N = 1+2^8 ).

I tried to look at k+flip(k) and there are some residuals of order eps(1) , but the residuals are also symmetrical, and the imaginary part of the FFT does not come out as fuzz of order eps(1) , but rather with a value comparable to the real part.

What dazzlingly obvious thing am I missing?

The dazzlingly obvious thing I was missing:

FFT is not an integral over the entire space; therefore, it receives a periodic signal. Above, I duplicate the last point in the period when I select even N , and therefore there is no way to shift it to put the zero frequency at the beginning without fractional indexing, which does not exist.

A word about my choice k . This is not arbitrary. The actual problem I'm trying to solve is to create a FTIR model interferogram, which I will FFT to get the spectrum. k is the distance by which the interferometer moves, which is converted into a frequency in wave numbers. In the real problem, there will be different scaling factors, so that the generating function V will give physically meaningful numbers.

Answer 1

@Yvon is absolutely right with his comments on symmetry. Your input signal looks symmetrical, but not because the symmetry is related to start 0. Using linspace in Matlab to build signals is in most cases a bad choice. Trying to restore results using fftshift is also a bad idea.

Use instead:

 k = 2*(0:N-1)/N - 1; 

and you get the expected result. However, the imaginary part of the converted values ​​will not be completely zero. There is some numerical noise.

 >> max(abs(imag(Vf5))) ans = 2.5535e-15 

The answer to Yvon’s question:

Why? → N = 1 + 2 ^ 4 N = 17 → x = linspace (-1,1, N) x = -1.0000 -0.8750 -0.7500 -0.6250 -0.5000 -0.3750 - 0.2500 -0.1250 0 0.1250 0.2500 0.3750 0.5000 0.6250 0.750 0.8750 1.0000 → y = 2 * (0: N-1) / N-1 y = -1 , 0000 -0.8824 -0.7647 -0.6471 -0.5294 -0.4118 -0.2941 -0.1765 -0.0588 0.0588 0.1765 0.291 0.4118 0.5294 0.641 0, 7647 0.824 - Yvon 1

Your example is not a symmetric (even) function, but an antisymmetric (odd) function. However, this does not matter.

For an antisymmetric function of length N, the following statement holds:

 f[i] == -f[-i] == -f[Ni] 

Index I works from 0 to N-1.

Let's see what happened with i = 2. Remember that count starts at 0 and ends at 16.

 x[2] = -0.75 -x[N-2] == -x[17-2] == -x[15] = (-1) 0.875 = -0.875 x[2] != -x[N-2] y[2] = -0.7647 -y[N-2] == -y[15] = (-1) 0.7647 y[2] == y[N-2] 

The problem is that the start of the Matlab vectors starts at 1. Modular (periodic) vectors start at the beginning 0. This difference leads to many misunderstandings.

Another way to explain why linspace (-1, + 1, N) is incorrect:

Imagine that you have a vector that contains one period of a periodic function, for example, the Cosine function. This single period is one of an infinite number of periods. The first value of your Cosinus vector must not be the same as the last value of your vector. However, this is exactly what linspace (-1, + 1, N) does. This leads to a sequence in which the last value of period 1 is the same value as the first pattern of the next period 2. This is not what you want. To avoid this error, use t = 2 * (0: N-1) / N - 1. The distance t [i + 1] -t [i] is 2 / N, and the last value should be t [N-1] = 1 - 2 / N, not 1.

Reply to Yvon's second comment

Whatever you insert the DFT / FFT into the input vector, in theory it is interpreted as a periodic function. But that is not the point.

DFT integrates.

 fft(m) = Sum_(k=0)^(N-1) (x(k) exp(-i 2 pi mk/N ) 

The first value x (k = 0) describes the amplitude of the first integration interval of length 1 / N. The second value x (k = 1) describes the amplitude of the second integration interval of length 1 / N. And so on.

The most recent integration interval of a symmetric function ends with the same value as the first sample. This means that the starting point of the last integration interval is k = N-1 = 1-1 / N. Your input vector contains the start points of the integration intervals.

Therefore, the last point of symmetry k = N is the point of the function, but it is not the starting point of the integration interval and therefore is not a member of the input vector.

Answer 2

it

 Vf = fftshift(fft(ifftshift(V))); 

That is, you need ifftshift in the time domain so that the patterns are interpreted as symbols of a symmetric function, and then fftshift in the frequency domain to make the symmetry again.

This only works for N odd. For N even concept of a symmetric function does not make sense: there is no way to shift the signal so that it is symmetrical relative to the origin (the origin must be “between two samples”, which is equally impossible).

In your example V , the code above gives Vf real and symmetrical. The following figure was created using semilogy(Vf) so that you can see both small and large values. (Of course, you can change the horizontal axis so that the graph is centered at 0 frequency, as you would expect, but in any case, the graph is considered symmetrical.)

Answer 3

You have a problem implementing the concept of "symmetry." A purely real, even (or "symmetric") function has a Fourier transform function, which is also real and even. "Even" is the symmetry about the y axis or t = 0.

However, when implementing a signal in Matlab, you always start with t = 0. That is, it is impossible to "determine" a signal starting from the beginning of time.

Searching the Internet for a while will lead me to this - Proper use of fftshift and ifftshift at the input to fft and ifft .

As Louis noted, you need to do ifftshift before giving the signal to fft . The reason has never been documented in Matlab, but only in this thread. For historical reasons, the outputs AND inputs of fft and ifft are reversed. That is, instead of ordered from -N/2 to N/2-1 (natural order), the signal in the time or frequency domain is ordered from 0 to N/2-1 , and then from -N/2 to -1 . This means that the correct coding method is fft( ifftshift(V) ) , but most people ignore this in most cases . Why he silently ignored, and did not raise huge problems, is that most of the problems were put on the amplitude of the signal, and not on the phase. Since the circular offset does not affect the amplitude spectrum, this is not a problem (even for the Matlab guys who wrote the documentation).

To check the equality of amplitude -

 Vf2 = fft(ifftshift(V)); Vf5 = fft(V); Va2 = abs(fftshift(Vf2)); Va5 = abs(fftshift(Vf5)); >> min(abs(Va2-Va5)<1e-10) ans = 1 

To see how badly the phase is wrong -

 Vp2 = angle(fftshift(Vf2)); Vp5 = angle(fftshift(Vf5)); 

In any case, as I wrote in the comment, after copying and pasting the code into a new pure Matlab, it gives 0 1 0 1 0 0 .

To your question about N = even and N = odd, I believe that when N = even, the signal is not symmetrical, since there is an unequal number of points on either side of the origin.

Answer 4

Just add the following line after "k = linspace (-1,1, N);"

 k(end)=[]; 

it will remove the last element of the array. This is defined as a symmetric array.

also believe that isreal (complex (1,0)) is false !!! The isreal function simply checks the memory format. so 1 + 0i is not real in the above example.

You defined your function to check for real numbers (like this)

 myisreal=@ (x) all((abs(imag(x))<1e-6*abs(real(x)))|(abs(x)<1e-8)); 

Finally, your source code should look something like this:

 N = 1 + 2^8; k = linspace(-1,1,N); k(end)=[]; V = exp(-abs(k)); Vf1 = fft(fftshift(V)); Vf2 = fft(ifftshift(V)); Vf3 = ifft(fftshift(V)); Vf4 = ifft(ifftshift(V)); Vf5 = fft(V); Vf6 = ifft(V); myisreal=@ (x) all((abs(imag(x))<1e-6*abs(real(x)))|(abs(x)<1e-8)); disp([myisreal(Vf1) myisreal(Vf2) myisreal(Vf3) myisreal(Vf4) myisreal(Vf5) myisreal(Vf6)]);