General Data Analysis Facilities

Fourier Analysis


The STSDAS package
fourier contains tasks to compute forward and reverse Fourier transforms, auto-correlation and cross-correlation, power-spectra, and convolution. In addition, the package provides some test and diagnostic functions to help in understanding and visualizing Fourier transforms. The task mkfunc allows you to create your own data sets for experimentation, shift lets you shift them in phase, and ftplot provides a convenient graphics display of both input and output data (real and imaginary parts). bracewell demonstrates Fourier transforms from the catalog of standard transform pairs in Bracewell's book on Fourier transforms. carith does complex multiplication or division of images.

Tasks in this package can transform data arrays of arbitrary size, i.e., there is no restriction for the arrays to be a power of 2 in length. The tasks are layered on code that does a prime factor decomposition of the Fourier transform. Arrays that cannot be decomposed into prime factors will not be transformed very quickly (i.e., if your data file is 509 x 509 pixels, you should probably copy it to a 512 x 512 array; if it is 510 x 510, with prime factors 2, 3, 5, 17, it's okay). You can use the tasks listprimes and factor to check the size of your data set.

The Fourier transform tasks make a reasonable attempt to keep track of the coordinates attached to your input data, and a number of Fourier transform coordinate pairs are recognized: time/freq, lambda/wavenumb, ll-mm/uu-vv.

The forward task performs a Fourier transform, defined by:

The inverse task performs an inverse Fourier transform, defined by:

The autocorr task performs auto-correlation, defined by:

The crosscorr task performs cross-correlation defined by:

The fconvolve task performs convolution, as defined by:

The powerspec task computes a power spectrum, defined by:

A frequently asked question about tasks in the fourier package is, "Which pixel receives the DC signal when taking a forward transform or power spectrum? That is, which pixel is zero frequency?" Zero frequency is at the reference pixel, which is given by the value of the image header keyword CRPIX1 (or CRPIX1 and CRPIX2 for a two-dimensional image). This will be pixel 1 (or [1,1] for two-dimensional) if the center parameter was set to "no".

The coordinate parameters of the input image are used to compute appropriate coordinate parameters for the output image. This makes it straightforward to associate a "world coordinate" value with any pixel in the output image. The coordinates can be seen using igi or splot for a one-dimensional image. For two-dimensional images, the image can be displayed using SAOimage, and then coordinates of features can be found using the rimcursor task with the parameter setting wcs="world". For an image produced by the forward or powerspec task, the world coordinate at a pixel is the frequency in physical units, such as cycles per second. As an example, consider an image N pixels long with pixel spacing of D seconds per pixel. The pixel spacing in the Fourier domain is 1 / (N D) Hz per pixel. Suppose the input image contains a periodic signal with period P pixels, which is P D seconds. The periodic signal results in a peak in the Fourier domain at N / P pixels from the reference pixel. The world coordinate (frequency in Hertz, in this example) at a pixel is the offset from the reference pixel multiplied by the pixel spacing, or (N / P) 1 / (N D) = 1 / (P D) Hz, as you would expect for a period of P D seconds.

For a two-dimensional image the pixel spacing changes independently in each dimension. If the length and pixel spacing are N1 and D1 for the first image axis and N2 and D2 for the second image axis, then the pixel spacing in the Fourier domain is 1 / (N1 D1) for the first axis and 1 / (N2 D2) for the second axis. In general, therefore, the coordinate axes will not be orthogonal in the Fourier domain. Tasks such as forward, inverse, and powerspec compute the output CD matrix by multiplying the input CD matrix on the right by the following matrix:



The interpretation of a power spectrum deserves some discussion. The power spectrum shows the variance of the input data. While the total variance is just the average of the squares of the data values, the power spectrum tells whether the variance is due primarily to low frequencies or high frequencies or what combination. The integral of the power spectrum over a given interval is the contribution to the variance attributed to frequencies within that interval. For discrete data, we can interpret that principle as follows. The sum of the squares of the pixel values of the input data is equal to the sum of all the pixel values of the power spectrum. (This includes negative frequencies.) So the contribution of a range of frequencies to the variance of the input data is the sum of the power spectrum over the range of pixel numbers corresponding to those frequencies, divided by the total number of pixels. We can translate between pixel number P and frequency F in the usual way using the coordinate parameters in the power spectrum, e.g., for one-dimensional data, F = CD1_1 (P - CRPIX1). The value of the power spectrum at zero frequency (i.e., at the reference pixel) is the square of the average of the input data multiplied by the number of pixels.

When the input data contain a periodic term of sufficient amplitude, it will show up as a feature in the power spectrum. The amplitude of that periodic term can be obtained by adding up the counts in the power spectrum over the feature, as follows:



where N is the number of pixels in the input data. The factor of two accounts for the mirror image feature at negative frequency, so the sum of counts is taken only over the feature at positive frequency. Of course, when taking the sum of counts over the feature, the background level of the power spectrum in the neighborhood of the feature should first be subtracted.

The units of a power spectrum can be determined from the units of the input data. An integral over a portion of a power spectrum, as described above, will have the same units as the variance of the input data, e.g., (counts/pixel)2. In the particular case of a one-dimensional time series with 0.2 second per pixel, if one of the input data values is 20 (counts/pixel), that value could also be regarded as 100 counts/sec. So if the variance has been computed in units of (counts/pixel)2 by integrating over the power spectrum, the variance may be converted to (counts/sec)2 by dividing by (0.2)2.


Generated with WebMaker