FFT(3) User Contributed Perl Documentation FFT(3)NAMEMath::FFT - Perl module to calculate Fast Fourier Transforms
SYNOPSIS
use Math::FFT;
my $PI = 3.1415926539;
my $N = 64;
my ($series, $other_series);
for (my $k=0; $k<$N; $k++) {
$series->[$k] = sin(4*$k*$PI/$N) + cos(6*$k*$PI/$N);
}
my $fft = new Math::FFT($series);
my $coeff = $fft->rdft();
my $spectrum = $fft->spctrm;
my $original_data = $fft->invrdft($coeff);
for (my $k=0; $k<$N; $k++) {
$other_series->[$k] = sin(16*$k*$PI/$N) + cos(8*$k*$PI/$N);
}
my $other_fft = $fft->clone($other_series);
my $other_coeff = $other_fft->rdft();
my $correlation = $fft->correl($other_fft);
DESCRIPTION
This module implements some algorithms for calculating Fast Fourier
Transforms for one-dimensional data sets of size 2^n. The data,
assumed to arise from a constant sampling rate, is represented by an
array reference $data (as described in the methods below), which is
then used to create a "Math::FFT" object as
my $fft = new Math::FFT($data);
The methods available include the following.
FFT METHODS
"$coeff = $fft->cdft();"
This calculates the complex discrete Fourier transform for a data
set "x[j]". Here, $data is a reference to an array
"data[0...2*n-1]" holding the data
data[2*j] = Re(x[j]),
data[2*j+1] = Im(x[j]), 0<=j<n
An array reference $coeff is returned consisting of
coeff[2*k] = Re(X[k]),
coeff[2*k+1] = Im(X[k]), 0<=k<n
where
X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
"$orig_data = $fft->invcdft([$coeff]);"
Calculates the inverse complex discrete Fourier transform on a data
set "x[j]". If $coeff is not given, it will be set equal to an
earlier call to "$fft->cdft()". $coeff is a reference to an array
"coeff[0...2*n-1]" holding the data
coeff[2*j] = Re(x[j]),
coeff[2*j+1] = Im(x[j]), 0<=j<n
An array reference $orig_data is returned consisting of
orig_data[2*k] = Re(X[k]),
orig_data[2*k+1] = Im(X[k]), 0<=k<n
where, excluding the scale,
X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
A scaling "$orig_data->[$i] *= 2.0/$n" is then done so that
$orig_data coincides with the original $data.
"$coeff = $fft->rdft();"
This calculates the real discrete Fourier transform for a data set
"x[j]". On input, $data is a reference to an array "data[0...n-1]"
holding the data. An array reference $coeff is returned consisting
of
coeff[2*k] = R[k], 0<=k<n/2
coeff[2*k+1] = I[k], 0<k<n/2
coeff[1] = R[n/2]
where
R[k] = sum_j=0^n-1 data[j]*cos(2*pi*j*k/n), 0<=k<=n/2
I[k] = sum_j=0^n-1 data[j]*sin(2*pi*j*k/n), 0<k<n/2
"$orig_data = $fft->invrdft([$coeff]);"
Calculates the inverse real discrete Fourier transform on a data
set "coeff[j]". If $coeff is not given, it will be set equal to an
earlier call to "$fft->rdft()". $coeff is a reference to an array
"coeff[0...n-1]" holding the data
coeff[2*j] = R[j], 0<=j<n/2
coeff[2*j+1] = I[j], 0<j<n/2
coeff[1] = R[n/2]
An array reference $orig_data is returned where, excluding the
scale,
orig_data[k] = (R[0] + R[n/2]*cos(pi*k))/2 +
sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) +
sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
A scaling "$orig_data->[$i] *= 2.0/$n" is then done so that
$orig_data coincides with the original $data.
"$coeff = $fft->ddct();"
Computes the discrete cosine tranform on a data set "data[0...n-1]"
contained in an array reference $data. An array reference $coeff is
returned consisting of
coeff[k] = C[k], 0<=k<n
where
C[k] = sum_j=0^n-1 data[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
"$orig_data = $fft->invddct([$coeff]);"
Computes the inverse discrete cosine tranform on a data set
"coeff[0...n-1]" contained in an array reference $coeff. If $coeff
is not given, it will be set equal to an earlier call to
"$fft->ddct()". An array reference $orig_data is returned
consisting of
orig_data[k] = C[k], 0<=k<n
where, excluding the scale,
C[k] = sum_j=0^n-1 coeff[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
A scaling "$orig_data->[$i] *= 2.0/$n" is then done so that
$orig_data coincides with the original $data.
"$coeff = $fft->ddst();"
Computes the discrete sine transform of a data set "data[0...n-1]"
contained in an array reference $data. An array reference $coeff is
returned consisting of
coeff[k] = S[k], 0<k<n
coeff[0] = S[n]
where
S[k] = sum_j=0^n-1 data[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
"$orig_data = $fft->invddst($coeff);"
Computes the inverse discrete sine transform of a data set
"coeff[0...n-1]" contained in an array reference $coeff, arranged
as
coeff[j] = A[j], 0<j<n
coeff[0] = A[n]
If $coeff is not given, it will be set equal to an earlier call to
"$fft->ddst()". An array reference $orig_data is returned
consisting of
orig_data[k] = S[k], 0<=k<n
where, excluding a scale,
S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
The scaling "$a->[$i] *= 2.0/$n" is then done so that $orig_data
coincides with the original $data.
"$coeff = $fft->dfct();"
Computes the real symmetric discrete Fourier transform of a data
set "data[0...n]" contained in the array reference $data. An array
reference $coeff is returned consisting of
coeff[k] = C[k], 0<=k<=n
where
C[k] = sum_j=0^n data[j]*cos(pi*j*k/n), 0<=k<=n
"$orig_data = $fft->invdfct($coeff);"
Computes the inverse real symmetric discrete Fourier transform of a
data set "coeff[0...n]" contained in the array reference $coeff.
If $coeff is not given, it will be set equal to an earlier call to
"$fft->dfct()". An array reference $orig_data is returned
consisting of
orig_data[k] = C[k], 0<=k<=n
where, excluding the scale,
C[k] = sum_j=0^n coeff[j]*cos(pi*j*k/n), 0<=k<=n
A scaling "$coeff->[0] *= 0.5", "$coeff->[$n] *= 0.5", and
"$orig_data->[$i] *= 2.0/$n" is then done so that $orig_data
coincides with the original $data.
"$coeff = $fft->dfst();"
Computes the real anti-symmetric discrete Fourier transform of a
data set "data[0...n-1]" contained in the array reference $data. An
array reference $coeff is returned consisting of
coeff[k] = C[k], 0<k<n
where
C[k] = sum_j=0^n data[j]*sin(pi*j*k/n), 0<k<n
("coeff[0]" is used for a work area)
"$orig_data = $fft->invdfst($coeff);"
Computes the inverse real anti-symmetric discrete Fourier transform
of a data set "coeff[0...n-1]" contained in the array reference
$coeff. If $coeff is not given, it will be set equal to an earlier
call to "$fft->dfst()". An array reference $orig_data is returned
consisting of
orig_data[k] = C[k], 0<k<n
where, excluding the scale,
C[k] = sum_j=0^n coeff[j]*sin(pi*j*k/n), 0<k<n
A scaling "$orig_data->[$i] *= 2.0/$n" is then done so that
$orig_data coincides with the original $data.
CLONING
The algorithm used in the transforms makes use of arrays for a work
area and for a cos/sin lookup table dependent only on the size of the
data set. These arrays are initialized when the "Math::FFT" object is
created and then are populated when a transform method is first
invoked. After this, they persist for the lifetime of the object.
This aspect is exploited in a "cloning" method; if a "Math::FFT" object
is created for a data set $data1 of size "N":
$fft1 = new Math::FFT($data1);
then a new "Math::FFT" object can be created for a second data set
$data2 of the same size "N" by
$fft2 = $fft1->clone($data2);
The $fft2 object will copy the reuseable work area and lookup table
calculated from $fft1.
APPLICATIONS
This module includes some common applications - correlation,
convolution and deconvolution, and power spectrum - that arise with
real data sets. The conventions used here follow that of Numerical
Recipes in C, by Press, Teukolsky, Vetterling, and Flannery, in which
further details of the algorithms are given. Note in particular the
treatment of end effects by zero padding, which is assumed to be done
by the user, if required.
Correlation
The correlation between two functions is defined as
/
Corr(t) = | ds g(s+t) h(s)
/
This may be calculated, for two array references $data1 and $data2
of the same size $n, as either
$fft1 = new Math::FFT($data1);
$fft2 = new Math::FFT($data2);
$corr = $fft1->correl($fft2);
or as
$fft1 = new Math::FFT($data1);
$corr = $fft1->correl($data2);
The array reference $corr is returned in wrap-around order -
correlations at increasingly positive lags are in "$corr->[0]"
(zero lag) on up to "$corr->[$n/2-1]", while correlations at
increasingly negative lags are in "$corr->[$n-1]" on down to
"$corr->[$n/2]". The sign convention used is such that if $data1
lags $data2 (that is, is shifted to the right), then $corr will
show a peak at positive lags.
Convolution
The convolution of two functions is defined as
/
Convlv(t) = | ds g(s) h(t-s)
/
This is similar to calculating the correlation between the two
functions, but typically the functions here have a quite different
physical interpretation - one is a signal which persists
indefinitely in time, and the other is a response function of
limited duration. The convolution may be calculated, for two array
references $data and $respn, as
$fft = new Math::FFT($data);
$convlv = $fft->convlv($respn);
with the returned $convlv being an array reference. The method
assumes that the response function $respn has an odd number of
elements $m less than or equal to the number of elements $n of
$data. $respn is assumed to be stored in wrap-around order - the
first half contains the response at positive times, while the
second half, counting down from "$respn->[$m-1]", contains the
response at negative times.
Deconvolution
Deconvolution undoes the effects of convoluting a signal with a
known response function. In other words, in the relation
/
Convlv(t) = | ds g(s) h(t-s)
/
deconvolution reconstructs the original signal, given the
convolution and the response function. The method is implemented,
for two array references $data and $respn, as
$fft = new Math::FFT($data);
$deconvlv = $fft->deconvlv($respn);
As a result, if the convolution of a data set $data with a response
function $respn is calculated as
$fft1 = new Math::FFT($data);
$convlv = $fft1->convlv($respn);
then the deconvolution
$fft2 = new Math::FFT($convlv);
$deconvlv = $fft2->deconvlv($respn);
will give an array reference $deconvlv containing the same elements
as the original data $data.
Power Spectrum
If the FFT of a real function of "N" elements is calculated, the
"N/2+1" elements of the power spectrum are defined, in terms of the
(complex) Fourier coefficients "C[k]", as
P[0] = |C[0]|^2 / N^2
P[k] = 2 |C[k]|^2 / N^2 (k = 1, 2 ,..., N/2-1)
P[N/2] = |C[N/2]|^2 / N^2
Often for these purposes the data is partitioned into "K" segments,
each containing "2M" elements. The power spectrum for each segment
is calculated, and the net power spectrum is the average of all of
these segmented spectra.
Partitioning may be done in one of two ways: non-overlapping and
overlapping. Non-overlapping is useful when the data set is
gathered in real time, where the number of data points can be
varied at will. Overlapping is useful where there is a fixed number
of data points. In non-overlapping, the first <2M> elements
constitute segment 1, the next "2M" elements are segment 2, and so
on up to segment "K", for a total of "2KM" sampled points. In
overlapping, the first and second "M" elements are segment 1, the
second and third "M" elements are segment 2, and so on, for a total
of "(K+1)M" sampled points.
A problem that may arise in this procedure is leakage: the power
spectrum calculated for one bin contains contributions from nearby
bins. To lessen this effect data windowing is often used: multiply
the original data "d[j]" by a window function "w[j]", where j = 0,
1, ..., N-1. Some popular choices of such functions are
| j - N/2 |
w[j] = 1 - | ------- | ... Bartlett
| N/2 |
/ j - N/2 \ 2
w[j] = 1 - | ------- | ... Welch
\ N/2 /
1 / \
w[j] = --- |1 - cos(2 pi j / N) | ... Hann
2 \ /
The "spctrm" method, used as
$fft = Math::FFT->new($data);
$spectrum = $fft->spctrm(%options);
returns an array reference $spectrum representing the power
spectrum for a data set represented by an array reference $data.
The options available are
"window => window_name"
This specifies the window function; if not given, no such
function is used. Accepted values (see above) are "bartlett",
"welch", "hann", and "\&my_window", where "my_window" is a user
specified subroutine which must be of the form, for example,
sub my_window {
my ($j, $n) = @_;
return 1 - abs(2*($j-$n/2)/$n);
}
which implements the Bartlett window.
"overlap => 1"
This specifies whether overlapping should be done; if true (1),
overlapping will be used, whereas if false (0), or not
specified, no overlapping is used.
"segments => n"
This specifies that the data will be partitioned into "n"
segments. If not specified, no segmentation will be done.
"number => m"
This specifies that "2m" data points will be used for each
segment, and must be a power of 2. The power spectrum returned
will consist of "m+1" elements.
STATISTICAL FUNCTIONS
For convenience, a number of common statistical functions are included
for analyzing real data. After creating the object as
my $fft = new Math::FFT($data);
for a data set represented by the array reference $data of size "N",
these methods may be called as follows.
"$mean = $fft->mean([$data]);"
This returns the mean
1/N * sum_j=0^N-1 data[j]
If an array reference $data is not given, the data set used in
creating $fft will be used.
"$stdev = $fft->stdev([$data]);"
This returns the standard deviation
sqrt{ 1/(N-1) * sum_j=0^N-1 (data[j] - mean)**2 }
If an array reference $data is not given, the data set used in
creating $fft will be used.
"$rms = $fft->rms([$data]);"
This returns the root mean square
sqrt{ 1/N * sum_j=0^N-1 (data[j])**2 }
If an array reference $data is not given, the data set used in
creating $fft will be used.
"($min, $max) = $fft->range([$data]);"
This returns the minimum and maximum values of the data set. If an
array reference $data is not given, the data set used in creating
$fft will be used.
"$median = $fft->median([$data]);"
This returns the median of a data set. The median is defined, for
the 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 side of the middle, if the number of elements is even. If
an array reference $data is not given, the data set used in
creating $fft will be used.
BUGS
Please report any to Randy Kobes <randy@theoryx5.uwinnipeg.ca>
SEE ALSO
Math::Pari and PDL
COPYRIGHT
The algorithm used in this module to calculate the Fourier transforms
is based on the C routine of fft4g.c available at
http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html, which is copyrighted
1996-99 by Takuya OOURA. The file arrays.c included here to 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
copyright 2000,2005 by Randy Kobes <r.kobes@uwinnipeg.ca>, and is
distributed under the same terms as Perl itself.
perl v5.14.0 2005-09-03 FFT(3)
