Dicrete Fourier Transform: how to use fftshift with fft properly

Question

I want to calculate the FFT numerically on the set Y. For testing, I use the Gauss function Y = exp (-x ^ 2). The symbolic Fourier transform is Y '= constant * exp (-k ^ 2/4).

import numpy X = numpy.arange(-100,100) Y = numpy.exp(-(X/5.0)**2) 

The naive approach fails:

 from numpy.fft import * from matplotlib import pyplot def plotReIm(x,y): f = pyplot.figure() ax = f.add_subplot(111) ax.plot(x, numpy.real(y), 'b', label='R()') ax.plot(x, numpy.imag(y), 'r:', label='I()') ax.plot(x, numpy.abs(y), 'k--', label='abs()') ax.legend() Y_k = fftshift(fft(Y)) k = fftshift(fftfreq(len(Y))) plotReIm(k,Y_k) 

real (Y_k) jumps between positive and negative values, which correspond to the phase of the jump, which is absent in the symbolic result. This, of course, is undesirable. (The result is technically correct in the sense that abs (Y_k) gives the amplitudes as expected if if (Y_k) equals Y.)

Here the function fftshift () displays the array k monotonically increasing and accordingly changes Y_k. The pair zip (k, Y_k) does not change, applying this operation to both vectors.

These changes fix the problem:

 Y_k = fftshift(fft(ifftshift(Y))) k = fftshift(fftfreq(len(Y))) plotReIm(k,Y_k) 

Is it correct to use the fft () function if monotone Y and Y_k are required?

Reverse operation above:

 Yx = fftshift(ifft(ifftshift(Y_k))) x = fftshift(fftfreq(len(Y_k), k[1] - k[0])) plotReIm(x,Yx) 

In this case, the document clearly indicates that Y_k must be sorted compatible with the outputs fft () and fftfreq (), which we can achieve by using ifftshift ().

These questions have bothered me for a long time: Are the output and input arrays of both fft () and ifft () always such that a[0] should contain the zero frequency term, a[1:n/2+1] should contain the positive-frequency terms, and a[n/2+1:] should contain the negative-frequency terms, in order of decreasingly negative frequency [numpy reference], where "frequency" is an independent variable?

The answer to the Fourier transform of a Gaussian is not Gaussian does not answer my question.

Answer 1

FFT can be considered as the creation of set vectors, each of which has an amplitude and a phase. The fft_shift operation changes the reference point for the phase angle of zero, from the edge of the FFT aperture, to the center of the original input data vector.

The phase (and, therefore, the real component of the complex vector) of the result sometimes turns out to be less “acne” when this is done, especially if any input function is completed in such a way that it is discontinuous along the edges of the FFT aperture. Or, if the input is symmetrical around the center of the FFT aperture, the phase of the FFT result will always be zero after fft_shift.

Fft_shift can be accomplished by rotating the N / 2 vector or simply flipping alternating sign bits as a result of the FFT, which can be more dcache friendly to the CPU.

Answer 2

The definition of the fft (and ifft ) ifft is here: http://docs.scipy.org/doc/numpy/reference/routines.fft.html#background-information

This is what routines compute, no more and no less. Note that the discrete Fourier transform is quite different from the continuous Fourier transform. For a dense-sampling function, there is a relationship between the two, but the ratio also includes phase factors and scaling in addition to fftshift . This is the reason for the fluctuations that you see in your plot. You can solve the necessary phase factor from the above mathematical formula for the DFT.