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.