JMSLTM Numerical Library 4.0

com.imsl.math
Class FFT

java.lang.Object
  extended bycom.imsl.math.FFT
All Implemented Interfaces:
Cloneable, Serializable

public class FFT
extends Object
implements Serializable, Cloneable

FFT functions.

Class FFT computes the discrete Fourier transform of a real vector of size n. The method used is a variant of the Cooley-Tukey algorithm, which is most efficient when n is a product of small prime factors. If n satisfies this condition, then the computational effort is proportional to n log n.

The forward method computes the forward transform. If n is even, then the forward transform is

q_{2m - 1} = sumlimits_{k = 0}^{n - 1} 
  {p_k } cos frac{{2pi km}}{n} ,,,, m = 1,; ldots ,;n/2

q_{2m - 2}   =  - sumlimits_{k = 0}^{n - 1} 
  {p_k } sin frac{{2pi km}}{n} ,,,, m = 1,; ldots ,;n/2 - 1

q_0   = sumlimits_{k = 0}^{n - 1} {p_k }

If n is odd, q_m is defined as above for m from 1 to (n - 1)/2.

Let f be a real valued function of time. Suppose we sample f at n equally spaced time intervals of length delta seconds starting at time t_0. That is, we have

p_i : = fleft( {t_0  + iDelta } right),i = 0,,1, 
  ldots ,,n - 1

We will assume that n is odd for the remainder of this discussion. The class FFT treats this sequence as if it were periodic of period n. In particular, it assumes that fleft( {t_0 } right) = fleft( {t_0  + nDelta } 
  right). Hence, the period of the function is assumed to be T = nDelta. We can invert the above transform for p as follows:

p_m  = {1 over n}left[ {q_0  + 
  2sumlimits_{k = 0}^{left( {n - 3} right)/2} {quad q_{2k + 1} } cos 
  {{2pi (k+1)m} over n} - 2sumlimits_{k = 0}^{left( {n - 3} right)/2} 
  {quad q_{2k + 2} } sin {{2pi (k+1)m} over n}} right]

This formula is very revealing. It can be interpreted in the following manner. The coefficients q produced by FFT determine an interpolating trigonometric polynomial to the data. That is, if we define

gleft( t right) = {1 over n}left[ 
  {q_0  + 2sumlimits_{k = 0}^{left( {n - 3} right)/2} {quad q_{2k + 1} } 
  cos {{2pi (k+1)left( {t - t_0 } right)} over {nDelta }} - 
  2sumlimits_{k = 0}^{left( {n - 3} right)/2} {quad q_{2k + 2} } 
  sin {{2pi (k+1)left( {t - t_0 } right)} over {nDelta }}} right]

= {1 over n}left[ {q_0  + 2sumlimits_{k = 0}^{left( 
  {n - 3} right)/2} {quad q_{2k + 1} } cos {{2pi (k+1)left( {t - t_0 } 
  right)} over T} - 2sumlimits_{k = 0}^{left( {n - 3} right)/2} {quad 
  q_{2k + 2} } sin {{2pi (k+1)left( {t - t_0 } right)} over T}} right]

then we have

fleft( {{rm{t}}_{rm{0}}  + left( {i - 1} 
  right)Delta } right) = gleft( {{rm{t}}_{rm{0}}  + left( {i - 1} 
  right)} right)Delta

Now suppose we want to discover the dominant frequencies, forming the vector P of length (n + 1)/2 as follows:

P_0: = left| {q_0 } right|

P_k: = sqrt {q_{2k - 2}^2  + q_{2k - 1}^2 } 
  ,,,, k = 1,;2,; ldots ,;left( {n - 1} right)/2

These numbers correspond to the energy in the spectrum of the signal. In particular, P_k corresponds to the energy level at frequency

{k over T} = {k over {nDelta }} 
  ,,,,,,,, k = 0,;1,; ldots ,;{{n - 1} over 2}

Furthermore, note that there are only (n + 1)/2 approx T/(2Delta) resolvable frequencies when n observations are taken. This is related to the Nyquist phenomenon, which is induced by discrete sampling of a continuous signal. Similar relations hold for the case when n is even.

If the backward method is used, then the backward transform is computed. If n is even, then the backward transform is

q_m  = p_0  + left( { - 1} right)^m p_{n - 1} 
  + 2sumlimits_{k = 0}^{n/2 - 1} {p_{2k + 1} } cos frac{{2pi (k+1)m}}{n} - 
  2sumlimits_{k = 0}^{n/2 - 2} {p_{2k + 2} } sin frac{{2pi (k+1)m}}{n}

If n is odd,

q_m  = p_0  + 2sumlimits_{k = 0}^{left( 
  {n - 3} right)/2} {;p_{2k + 1} } cos {{2pi (k+1)m} over n} - 
  2sumlimits_{k = 0}^{left( {n - 3} right)/2} {;p_{2k + 2} } sin 
  {{2pi (k+1)m} over n}

The backward Fourier transform is the unnormalized inverse of the forward Fourier transform.

FFT is based on the real FFT in FFTPACK, which was developed by Paul Swarztrauber at the National Center for Atmospheric Research.

See Also:
Example, Serialized Form

Constructor Summary
FFT(int n)
          Constructs an FFT object.
 
Method Summary
 double[] backward(double[] coef)
          Compute the real periodic sequence from its Fourier coefficients.
 double[] forward(double[] seq)
          Compute the Fourier coefficients of a real periodic sequence.
 
Methods inherited from class java.lang.Object
clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait
 

Constructor Detail

FFT

public FFT(int n)
Constructs an FFT object.

Parameters:
n - is the length of the sequence to be transformed
Method Detail

backward

public double[] backward(double[] coef)
Compute the real periodic sequence from its Fourier coefficients.

Parameters:
coef - a double array containing the Fourier coefficients
Returns:
a double array containing the periodic sequence

forward

public double[] forward(double[] seq)
Compute the Fourier coefficients of a real periodic sequence.

Parameters:
seq - a double array containing the sequence to be transformed
Returns:
a double array containing the transformed sequence

JMSLTM Numerical Library 4.0

Copyright 1970-2006 Visual Numerics, Inc.
Built June 1 2006.