Math::FFT - Perl module to calculate Fast Fourier Transforms
version 1.34
use Math::FFT; my $PI = 3.1415926539; my $N = 64; my $series = [map { sin(4*$_*$PI/$N) + cos(6*$_*$PI/$N) } 0 .. $N-1]; my $fft = Math::FFT->new($series); my $coeff = $fft->rdft(); my $spectrum = $fft->spctrm; my $original_data = $fft->invrdft($coeff); my $other_series = [map { sin(16*$_*$PI/$N) + cos(8*$_*$PI/$N) } 0 .. $N-1]; my $other_fft = $fft->clone($other_series); my $other_coeff = $other_fft->rdft(); my $correlation = $fft->correl($other_fft);
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 = Math::FFT->new($data);
The methods available include the following.
my $fft = Math::FFT->new($series)
The constructor. Pass it an array of numbers, with a length that is a power of 2.
$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 transform 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 transform 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
.
my $other_fft = $fft->clone($other_series)
See "CLONING" below.
my $other_series = $fft->convlv($response_data)
See "Convolution" below.
my $corr = $fft->correl($other_fft)
See "Correlation" below.
my $deconvlv = $fft->deconvlv($respn)
See "Deconvolution" below.
For internal use. Don't use directly.
For internal use. Don't use directly.
See "Power Spectrum" below.
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 = Math::FFT->new($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
.
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.
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 = Math::FFT->new($data1); $fft2 = Math::FFT->new($data2); $corr = $fft1->correl($fft2);
or as
$fft1 = Math::FFT->new($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.
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 = Math::FFT->new($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 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 = Math::FFT->new($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 = Math::FFT->new($data); $convlv = $fft1->convlv($respn);
then the deconvolution
$fft2 = Math::FFT->new($convlv); $deconvlv = $fft2->deconvlv($respn);
will give an array reference $deconvlv
containing the same elements as the original data $data
.
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.
For convenience, a number of common statistical functions are included for analyzing real data. After creating the object as
my $fft = Math::FFT->new($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.
Please report any to Randy Kobes <randy@theoryx5.uwinnipeg.ca>
Math::Pari and PDL
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.
Shlomi Fish <shlomif@cpan.org>
This software is copyright (c) 2000 by Randy Kobes.
This is free software; you can redistribute it and/or modify it under the same terms as the Perl 5 programming language system itself.
Please report any bugs or feature requests on the bugtracker website https://rt.cpan.org/Public/Dist/Display.html?Name=Math-FFT or by email to bug-math-fft@rt.cpan.org.
When submitting a bug or request, please include a test-file or a patch to an existing test-file that illustrates the bug or desired feature.
You can find documentation for this module with the perldoc command.
perldoc Math::FFT
The following websites have more information about this module, and may be of help to you. As always, in addition to those websites please use your favorite search engine to discover more resources.
A modern, open-source CPAN search engine, useful to view POD in HTML format.
The default CPAN search engine, useful to view POD in HTML format.
The RT ( Request Tracker ) website is the default bug/issue tracking system for CPAN.
The AnnoCPAN is a website that allows community annotations of Perl module documentation.
The CPAN Ratings is a website that allows community ratings and reviews of Perl modules.
The CPAN Forum is a web forum for discussing Perl modules.
The CPANTS is a website that analyzes the Kwalitee ( code metrics ) of a distribution.
The CPAN Testers is a network of smokers who run automated tests on uploaded CPAN distributions.
The CPAN Testers Matrix is a website that provides a visual overview of the test results for a distribution on various Perls/platforms.
The CPAN Testers Dependencies is a website that shows a chart of the test results of all dependencies for a distribution.
Please report any bugs or feature requests by email to bug-math-fft at rt.cpan.org
, or through the web interface at https://rt.cpan.org/Public/Bug/Report.html?Queue=Math-FFT. You will be automatically notified of any progress on the request by the system.
The code is open to the world, and available for you to hack on. Please feel free to browse it and play with it, or whatever. If you want to contribute patches, please send me a diff or prod me to pull from your repository :)
https://github.com/shlomif/perl-Math-FFT
git clone https://github.com/shlomif/perl-Math-FFT.git