NAME=20
SYN= OPSIS=20
= DESCRIPTION=20
- = FFT=20 METHODS=20
- CLON= ING=20
- APPLICATIONS=20
- STATISTICAL=20 FUNCTIONS
BUGS=20
SEE= =20 ALSO=20
CO= PYRIGHT=20

NAME=20

Math::FFT - Perl module to calculate Fast Fourier Transforms

SYNOPSIS

  use =
Math::FFT;
  my $PI =3D 3.1415926539;
  my $N =3D 64;
  my ($series, $other_series);
  for (my $k=3D0; $k<$N; $k++) {
      $series->[$k] =3D sin(4*$k*$PI/$N) + cos(6*$k*$PI/$N);
  }
  my $fft =3D new Math::FFT($series);
  my $coeff =3D $fft->rdft();
  my $spectrum =3D $fft->spctrm;
  my $original_data =3D $fft->invrdft($coeff);

  for (my $k=3D0; $k<$N; $k++) {
      $other_series->[$k] =3D sin(16*$k*$PI/$N) + cos(8*$k*$PI/$N);
  }
  my $other_fft =3D $fft->clone($other_series);
  my $other_coeff =3D $other_fft->rdft();
  my $correlation =3D $fft->correl($other_fft);

DESCRIPTION

This module implements some algorithms for calculating Fast Fourier=20 Transforms for one-dimensional data sets of size 2^n. The data, assumed = to arise=20 from a constant sampling rate, is represented by an array reference=20 $data (as described in the methods below), which is then = used to=20 create a Math::FFT object as

  my $fft =3D new =
Math::FFT($data);

The methods available include the following.

FFT METHODS

cdft();">$coeff =3D=20 $fft->cdft();=20

This calculates the complex discrete Fourier transform for a data = set=20 x[j]. Here, $data is a reference to an array = data[0...2*n-1] holding the data

  data[2*j] =3D =
Re(x[j]),
  data[2*j+1] =3D Im(x[j]), 0<=3Dj<n

An array reference $coeff is returned consisting = of

  coeff[2*k] =3D Re(X[k]),
  coeff[2*k+1] =3D Im(X[k]), 0<=3Dk<n

where

   X[k] =3D sum_j=3D0^n-1 x[j]*exp(2*pi*i*j*k/n), =
0<=3Dk<n

invcdft([$coeff]);">$orig_data =3D=20 $fft->invcdft([$coeff]);=20

Calculates the inverse complex discrete Fourier transform on a data = set=20 x[j]. If $coeff is not given, it will be set = equal=20 to an earlier call to $fft->cdft(). = $coeff is a=20 reference to an array coeff[0...2*n-1] holding the = data

  coeff[2*j] =3D Re(x[j]),
  coeff[2*j+1] =3D Im(x[j]), 0<=3Dj<n

An array reference $orig_data is returned consisting = of

  orig_data[2*k] =3D Re(X[k]),
  orig_data[2*k+1] =3D Im(X[k]), 0<=3Dk<n

where, excluding the scale,

   X[k] =3D sum_j=3D0^n-1 =
x[j]*exp(-2*pi*i*j*k/n), 0<=3Dk<n

A scaling $orig_data->[$i] *=3D 2.0/$n is then done = so that=20 $orig_data coincides with the original = $data.

rdft();">$coeff =3D=20 $fft->rdft();=20

This calculates the real discrete Fourier transform for a data set=20 x[j]. On input, $data is a reference to an array=20 data[0...n-1] holding the data. An array reference=20 $coeff is returned consisting of

  coeff[2*k] =3D =
R[k], 0<=3Dk<n/2
  coeff[2*k+1] =3D I[k], 0<k<n/2
  coeff[1] =3D R[n/2]

where

  R[k] =3D sum_j=3D0^n-1 data[j]*cos(2*pi*j*k/n), =
0<=3Dk<=3Dn/2
  I[k] =3D sum_j=3D0^n-1 data[j]*sin(2*pi*j*k/n), 0<k<n/2

invrdft([$coeff]);">$orig_data =3D=20 $fft->invrdft([$coeff]);=20

Calculates the inverse real discrete Fourier transform on a data = set=20 coeff[j]. If $coeff is not given, it will be = set=20 equal to an earlier call to $fft->rdft(). = $coeff=20 is a reference to an array coeff[0...n-1] holding the = data

  coeff[2*j] =3D R[j], 0<=3Dj<n/2
  coeff[2*j+1] =3D I[j], 0<j<n/2
  coeff[1] =3D R[n/2]

An array reference $orig_data is returned where, = excluding the=20 scale,

  orig_data[k] =3D (R[0] + R[n/2]*cos(pi*k))/2 +=20
    sum_j=3D1^n/2-1 R[j]*cos(2*pi*j*k/n) +=20
      sum_j=3D1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=3Dk<n

A scaling $orig_data->[$i] *=3D 2.0/$n is then done = so that=20 $orig_data coincides with the original = $data.

ddct();">$coeff =3D=20 $fft->ddct();=20

Computes the discrete cosine tranform on a data set=20 data[0...n-1] contained in an array reference = $data.=20 An array reference $coeff is returned consisting = of

  coeff[k] =3D C[k], 0<=3Dk<n

where

  C[k] =3D sum_j=3D0^n-1 data[j]*cos(pi*(j+1/2)*k/n), =
0<=3Dk<n

invddct([$coeff]);">$orig_data =3D=20 $fft->invddct([$coeff]);=20

Computes the inverse discrete cosine tranform on a data set=20 coeff[0...n-1] contained in an array reference=20 $coeff. If $coeff is not given, it will be = set equal=20 to an earlier call to $fft->ddct(). An array reference = $orig_data is returned consisting of

  =
orig_data[k] =3D C[k], 0<=3Dk<n

where, excluding the scale,

  C[k] =3D sum_j=3D0^n-1 =
coeff[j]*cos(pi*j*(k+1/2)/n), 0<=3Dk<n

A scaling $orig_data->[$i] *=3D 2.0/$n is then done = so that=20 $orig_data coincides with the original = $data.

ddst();">$coeff =3D=20 $fft->ddst();=20

Computes the discrete sine transform of a data set=20 data[0...n-1] contained in an array reference = $data.=20 An array reference $coeff is returned consisting = of

 coeff[k] =3D S[k], 0<k<n
 coeff[0] =3D S[n]

where

 S[k] =3D sum_j=3D0^n-1 data[j]*sin(pi*(j+1/2)*k/n), =
0<k<=3Dn

invddst($coeff);">$orig_data =3D=20 $fft->invddst($coeff);=20

Computes the inverse discrete sine transform of a data set=20 coeff[0...n-1] contained in an array reference=20 $coeff, arranged as

 coeff[j] =3D A[j], =
0<j<n
 coeff[0] =3D A[n]

If $coeff is not given, it will be set equal to an = earlier=20 call to $fft->ddst(). An array reference=20 $orig_data is returned consisting of

 =
orig_data[k] =3D S[k], 0<=3Dk<n

where, excluding a scale,

 S[k] =3D  sum_j=3D1^n =
A[j]*sin(pi*j*(k+1/2)/n), 0<=3Dk<n

The scaling $a->[$i] *=3D 2.0/$n is then done so = that=20 $orig_data coincides with the original = $data.

dfct();">$coeff =3D=20 $fft->dfct();=20

Computes the real symmetric discrete Fourier transform of a data = set=20 data[0...n] contained in the array reference = $data.=20 An array reference $coeff is returned consisting = of

  coeff[k] =3D C[k], 0<=3Dk<=3Dn

where

  C[k] =3D sum_j=3D0^n data[j]*cos(pi*j*k/n), =
0<=3Dk<=3Dn

invdfct($coeff);">$orig_data =3D=20 $fft->invdfct($coeff);=20

Computes the inverse real symmetric discrete Fourier transform of a = data=20 set coeff[0...n] contained in the array reference=20 $coeff. If $coeff is not given, it will be = set equal=20 to an earlier call to $fft->dfct(). An array reference = $orig_data is returned consisting of

  =
orig_data[k] =3D C[k], 0<=3Dk<=3Dn

where, excluding the scale,

  C[k] =3D sum_j=3D0^n =
coeff[j]*cos(pi*j*k/n), 0<=3Dk<=3Dn

A scaling $coeff->[0] *=3D 0.5, = $coeff->[$n] *=3D=20 0.5, and $orig_data->[$i] *=3D 2.0/$n is then = done so=20 that $orig_data coincides with the original=20 $data.

dfst();">$coeff =3D=20 $fft->dfst();=20

Computes the real anti-symmetric discrete Fourier transform of a = data set=20 data[0...n-1] contained in the array reference=20 $data. An array reference $coeff is returned = consisting of

  coeff[k] =3D C[k], 0<k<n

where

  C[k] =3D sum_j=3D0^n data[j]*sin(pi*j*k/n), =
0<k<n

(coeff[0] is used for a work area)

invdfst($coeff);">$orig_data =3D=20 $fft->invdfst($coeff);=20

Computes the inverse real anti-symmetric discrete Fourier transform = of a=20 data set coeff[0...n-1] contained in the array reference=20 $coeff. If $coeff is not given, it will be = set equal=20 to an earlier call to $fft->dfst(). An array reference = $orig_data is returned consisting of

  =
orig_data[k] =3D C[k], 0<k<n

where, excluding the scale,

  C[k] =3D sum_j=3D0^n =
coeff[j]*sin(pi*j*k/n), 0<k<n

A scaling $orig_data->[$i] *=3D 2.0/$n is then done = so that=20 $orig_data coincides with the original=20 $data.

CLONING

The algorithm used in the transforms makes use of arrays for a work = area and=20 for a cos/sin lookup table dependent only on the size of the data set. = These=20 arrays are initialized when the Math::FFT object is created = and=20 then are populated when a transform method is first invoked. After this, = they=20 persist for the lifetime of the object.

This aspect is exploited in a cloning method; if a=20 Math::FFT object is created for a data set = $data1 of=20 size N:

  $fft1 =3D new Math::FFT($data1);

then a new Math::FFT object can be created for a second = data set=20 $data2 of the same size N by

   =
$fft2 =3D $fft1->clone($data2);

The $fft2 object will copy the reuseable work area and = lookup=20 table calculated from $fft1.

APPLICATIONS

This module includes some common applications - correlation, = convolution and=20 deconvolution, and power spectrum - that arise with real data sets. The=20 conventions used here follow that of Numerical Recipes in C, by = Press,=20 Teukolsky, Vetterling, and Flannery, in which further details of the = algorithms=20 are given. Note in particular the treatment of end effects by zero = padding,=20 which is assumed to be done by the user, if required.

Correlation=20

The correlation between two functions is defined as

        =
     /
   Corr(t) =3D | ds g(s+t) h(s)=20
             /

This may be calculated, for two array references = $data1 and=20 $data2 of the same size $n, as = either

   $fft1 =3D new Math::FFT($data1);
   $fft2 =3D new Math::FFT($data2);
   $corr =3D $fft1->correl($fft2);

or as

   $fft1 =3D new Math::FFT($data1);
   $corr =3D $fft1->correl($data2);

The array reference $corr is returned in wrap-around = order -=20 correlations at increasingly positive lags are in = $corr->[0]=20 (zero lag) on up to $corr->[$n/2-1], while = correlations at=20 increasingly negative lags are in $corr->[$n-1] on = down to=20 $corr->[$n/2]. The sign convention used is such that = if=20 $data1 lags $data2 (that is, is shifted to = the=20 right), then $corr will show a peak at positive lags.

Convolution=20

The convolution of two functions is defined as

             =
  /
   Convlv(t) =3D | ds g(s) h(t-s)=20
               /

This is similar to calculating the correlation between the two = functions,=20 but typically the functions here have a quite different physical=20 interpretation - one is a signal which persists indefinitely in time, = and the=20 other is a response function of limited duration. The convolution may = be=20 calculated, for two array references $data and=20 $respn, as

   $fft =3D new Math::FFT($data);
   $convlv =3D $fft->convlv($respn);

with the returned $convlv being an array reference. = The method=20 assumes that the response function $respn has an = odd=20 number of elements $m less than or equal to the number of = elements $n of $data. $respn is = assumed=20 to be stored in wrap-around order - the first half contains the = response at=20 positive times, while the second half, counting down from=20 $respn->[$m-1], contains the response at negative = times.

Deconvolution=20

Deconvolution undoes the effects of convoluting a signal with a = known=20 response function. In other words, in the relation

            =
   /
   Convlv(t) =3D | ds g(s) h(t-s)=20
               /

deconvolution reconstructs the original signal, given the = convolution and=20 the response function. The method is implemented, for two array = references=20 $data and $respn, as

   $fft =3D new =
Math::FFT($data);
   $deconvlv =3D $fft->deconvlv($respn);

As a result, if the convolution of a data set $data = with a=20 response function $respn is calculated as

   =
$fft1 =3D new Math::FFT($data);
   $convlv =3D $fft1->convlv($respn);

then the deconvolution

   $fft2 =3D new Math::FFT($convlv);
   $deconvlv =3D $fft2->deconvlv($respn);

will give an array reference $deconvlv containing the = same=20 elements as the original data $data.

Power Spectrum=20

If the FFT of a real function of N elements is = calculated, the=20 N/2+1 elements of the power spectrum are defined, in = terms of the=20 (complex) Fourier coefficients C[k], as

   P[0] =
=3D |C[0]|^2 / N^2
   P[k] =3D 2 |C[k]|^2 / N^2   (k =3D 1, 2 ,..., N/2-1)
   P[N/2] =3D |C[N/2]|^2 / N^2

Often for these purposes the data is partitioned into = K=20 segments, each containing 2M elements. The power spectrum = for=20 each segment is calculated, and the net power spectrum is the average = of all=20 of these segmented spectra.

Partitioning may be done in one of two ways: non-overlapping = and=20 overlapping. Non-overlapping is useful when the data set is = gathered in=20 real time, where the number of data points can be varied at will. = Overlapping=20 is useful where there is a fixed number of data points. In = non-overlapping,=20 the first <2M> elements constitute segment 1, the next = 2M=20 elements are segment 2, and so on up to segment K, for a = total of=20 2KM sampled points. In overlapping, the first and second=20 M elements are segment 1, the second and third = M=20 elements are segment 2, and so on, for a total of (K+1)M = sampled=20 points.

A problem that may arise in this procedure is leakage: the = power=20 spectrum calculated for one bin contains contributions from nearby = bins. To=20 lessen this effect data windowing is often used: multiply the = original=20 data d[j] by a window function w[j], where j = =3D 0, 1,=20 ..., N-1. Some popular choices of such functions are

          =
    | j - N/2 |
  w[j] =3D 1 -  | ------- |     ... Bartlett  =20
              |   N/2   |


              / j - N/2 \ 2
  w[j] =3D 1 -  | ------- |     ... Welch =20
              \   N/2   /


           1   /                    \
  w[j] =3D  ---  |1 - cos(2 pi j / N) |     ... Hann =20
           2   \                    /

The spctrm method, used as

    $fft =3D =
Math::FFT->new($data);
    $spectrum =3D $fft->spctrm(%options);

returns an array reference $spectrum representing the = power=20 spectrum for a data set represented by an array reference = $data.=20 The options available are

_window_name">window =3D>=20 window_name=20

This specifies the window function; if not given, no such = function is=20 used. Accepted values (see above) are "bartlett",=20 "welch", "hann", and = \&my_window,=20 where my_window is a user specified subroutine which = must be of=20 the form, for example,

   sub my_window {
      my ($j, $n) =3D @_;
      return 1 - abs(2*($j-$n/2)/$n);
   }

which implements the Bartlett window.

_1">overlap = =3D>=20 1=20

This specifies whether overlapping should be done; if true (1),=20 overlapping will be used, whereas if false (0), or not specified, no = overlapping is used.

_n">segments =3D>=20 n=20

This specifies that the data will be partitioned into = n=20 segments. If not specified, no segmentation will be done.

_m">number = =3D> m=20

This specifies that 2m data points will be used for = each=20 segment, and must be a power of 2. The power spectrum returned will = consist=20 of m+1 elements.

STATISTICAL FUNCTIONS

For convenience, a number of common statistical functions are = included for=20 analyzing real data. After creating the object as

  my $fft =3D =
new Math::FFT($data);

for a data set represented by the array reference $data = of size=20 N, these methods may be called as follows.

mean([$data]);">$mean =3D=20 $fft->mean([$data]);=20

This returns the mean

  1/N * sum_j=3D0^N-1 data[j]

If an array reference $data is not given, the data set = used in=20 creating $fft will be used.

stdev([$data]);">$stdev =3D=20 $fft->stdev([$data]);=20

This returns the standard deviation

  sqrt{ 1/(N-1) * =
sum_j=3D0^N-1 (data[j] - mean)**2 }

If an array reference $data is not given, the data set = used in=20 creating $fft will be used.

rms([$data]);">$rms =3D=20 $fft->rms([$data]);=20

This returns the root mean square

  sqrt{ 1/N * =
sum_j=3D0^N-1 (data[j])**2 }

If an array reference $data is not given, the data set = used in=20 creating $fft will be used.

range([$data]);">($min, $max) =3D = $fft->range([$data]);=20

This returns the minimum and maximum values of the data set. If an = array=20 reference $data is not given, the data set used in = creating=20 $fft will be used.

median([$data]);">$median =3D=20 $fft->median([$data]);=20

This returns the median of a data set. The median is defined, for = the=20 sorted data set, as either the middle element, if the number of = elements is odd, or as the interpolated value of the the two values on = either=20 side of the middle, if the number of elements is even. If an array = reference=20 $data is not given, the data set used in creating=20 $fft will be used.

BUGS=20

Please report any to Randy Kobes = <randy@theoryx5.uwinnipeg.ca>

COPYRIGHT

The algorithm used in this module to calculate the Fourier transforms = is=20 based on the C routine of fft4g.c available at http://momonga.t.= u-tokyo.ac.jp/~ooura/fft.html,=20 which is copyrighted 1996-99 by Takuya OOURA. The file arrays.c included = here to=20 handle passing arrays to and from C comes from the PGPLOT module of Karl = Glazebrook <kgb@aaoepp.aao.gov.au>. The perl code of Math::FFT is=20 copyright 2000,2005 by Randy Kobes <r.kobes@uwinnipeg.ca>, and is=20 distributed under the same terms as Perl itself.