FFT METHODS

"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.

pdfct()

For internal use. Don't use directly.

pdfst()

For internal use. Don't use directly.

spctrm()

See "Power Spectrum" below.

CLONING

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.

APPLICATIONS

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 = 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.

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 = 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

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.

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

STATISTICAL FUNCTIONS

    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.

Websites

• MetaCPAN

  A modern, open-source CPAN search engine, useful to view POD in HTML format.

  <https://metacpan.org/release/Math-FFT>

• RT: CPAN's Bug Tracker

  The RT ( Request Tracker ) website is the default bug/issue tracking system for CPAN.

  <https://rt.cpan.org/Public/Dist/Display.html?Name=Math-FFT>

• CPANTS

  The CPANTS is a website that analyzes the Kwalitee ( code metrics ) of a distribution.

  <http://cpants.cpanauthors.org/dist/Math-FFT>

• CPAN Testers

  The CPAN Testers is a network of smoke testers who run automated tests on uploaded CPAN distributions.

  <http://www.cpantesters.org/distro/M/Math-FFT>

• CPAN Testers Matrix

  The CPAN Testers Matrix is a website that provides a visual overview of the test results for a distribution on various Perls/platforms.

  <http://matrix.cpantesters.org/?dist=Math-FFT>

• CPAN Testers Dependencies

  The CPAN Testers Dependencies is a website that shows a chart of the test results of all dependencies for a distribution.

  <http://deps.cpantesters.org/?module=Math::FFT>

Bugs / Feature Requests

Source Code

<https://github.com/shlomif/perl-Math-FFT>

  git clone git://github.com/shlomif/perl-Math-FFT.git

Websites

• MetaCPAN

  A modern, open-source CPAN search engine, useful to view POD in HTML format.

  <https://metacpan.org/release/Math-FFT>

• RT: CPAN's Bug Tracker

  The RT ( Request Tracker ) website is the default bug/issue tracking system for CPAN.

  <https://rt.cpan.org/Public/Dist/Display.html?Name=Math-FFT>

• CPANTS

  The CPANTS is a website that analyzes the Kwalitee ( code metrics ) of a distribution.

  <http://cpants.cpanauthors.org/dist/Math-FFT>

• CPAN Testers

  The CPAN Testers is a network of smoke testers who run automated tests on uploaded CPAN distributions.

  <http://www.cpantesters.org/distro/M/Math-FFT>

• CPAN Testers Matrix

  The CPAN Testers Matrix is a website that provides a visual overview of the test results for a distribution on various Perls/platforms.

  <http://matrix.cpantesters.org/?dist=Math-FFT>

• CPAN Testers Dependencies

  The CPAN Testers Dependencies is a website that shows a chart of the test results of all dependencies for a distribution.

  <http://deps.cpantesters.org/?module=Math::FFT>

Bugs / Feature Requests

Source Code

<https://github.com/shlomif/perl-Math-FFT>

  git clone git://github.com/shlomif/perl-Math-FFT.git