Program: FFTN Purpose: Apply fast Fourier transform in N dimensions, i.e. for a profile (N=1), a map or a series of maps (N=2) or a data cube (N=3). Category: MANIPULATION File: fftn.src Author: M.G.R. Vogelaar Keywords: INVERSE= FFT or Inverse FFT: ... [F]/I Set modus of operation: FFT or the inverse transformation. PHASEIN= Input data in amplitude and phase? ... Y/[N] If True, then the phases must be entered in radians. RINSET= Set, subset(s) real part of input: Note that if your input is not contiguous (e.g.: Aurora FREQ 10 14:16 20) then the output will be from FREQ 10 to 20 and the subsets that are not included in the input are set to blank in the output. The transform is only applied to the given subsets. IINSET= Set, subset imaginary part of input: BOX= Enter box in .... PHASEOUT= Output data in amplitude and phase? ... Y/[N] The output phases are in radians. The set names get the extension _amp and _phase Else they get extensions _real and _imag If INVERSE=n SHIFTIN= Set first element in sample to the center of the input data? ... Y/[N] SHIFTOUT= Shift the zero-freq. component of OUTPUT to the center of the spectrum? ... Y/[N] If INVERSE=y SHIFTIN= Shift the zero-freq. component of INPUT to the left of the spectrum? ... Y/[N] SHIFTOUT= Set first element of transform to be the center of the output data? ... Y/[N]" OUTSET= Name of output set: Description: The program calculates an N dimensional Fourier transform of the input data. The number of axes in the selected subset(s) sets the value of N. With this program it is possible to do a 3-dimensional FFT on a data cube or one or more 1-dimensional FFT's in an image or 2-dimensional FFT's on channel maps in a data cube etc. Gaussians --------- As a testcase, one often uses a Gaussian function as a real input. The FFT of this symmetrical real function should be also gaussian and the imaginary part is zero. An inverse FFT should return the original image. With numpy's fft() function (and therefore in this task), we need to do something extra to achieve the same result. First of all, the values in the result of an FFT follow so-called standard order: If A = fft(a, n), then A[0] contains the zero-frequency term (the mean of the signal), which is always purely real for real inputs. Then A[1:n/2] contains the positive-frequency terms, and A[n/2+1:] contains the negative- frequency terms, in order of decreasingly negative frequency. The routine fftshift(A) shifts transforms and their frequencies to put the zero-frequency components in the middle. This action is set with keyword SHIFTOUT=Y But if your input is a Gaussian with the peak in the center, your output wil not be Gaussian and the imaginary part is non-zero. The reason is that the input is not a real-valued, symmetrical time series because the peak is not at t=0. We have to shift the peak of the Gaussian to time = 0 and wrap around the shifted data for t < 0. The function that does this is called ifftshift(A) and is executed if we set SHIFTIN=Y From Gaussian with the peak in the center to Gaussian with the peak in the center, use: BOX= IINSET= INVERSE= n OUTSET= gaussout PHASEIN= PHASEOUT= RINSET= gaussin SHIFTIN= y SHIFTOUT= y To get the original Gaussian back: BOX= IINSET= gaussout_imag INVERSE= y OUTSET= igauss PHASEIN= PHASEOUT= RINSET= gaussout_real SHIFTIN= y SHIFTOUT= y Notes: If you shifted a transform and want to apply the inverse transform, you should also set SHIFTIN=Y because the input is not the FFT but the shifted FFT. Updates: Jul 30, 2014: VOG, Document created.