1Math::FFT(3)          User Contributed Perl Documentation         Math::FFT(3)
2
3
4

NAME

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

VERSION

9       version 1.36
10

SYNOPSIS

12           use Math::FFT;
13
14           my $PI = 3.1415926539;
15           my $N = 64;
16
17           my $series = [map { sin(4*$_*$PI/$N) + cos(6*$_*$PI/$N) } 0 .. $N-1];
18           my $fft = Math::FFT->new($series);
19           my $coeff = $fft->rdft();
20           my $spectrum = $fft->spctrm;
21           my $original_data = $fft->invrdft($coeff);
22
23           my $other_series =
24               [map { sin(16*$_*$PI/$N) + cos(8*$_*$PI/$N) } 0 .. $N-1];
25
26           my $other_fft = $fft->clone($other_series);
27           my $other_coeff = $other_fft->rdft();
28           my $correlation = $fft->correl($other_fft);
29

DESCRIPTION

31       This module implements some algorithms for calculating Fast Fourier
32       Transforms for one-dimensional data sets of size 2^n.  The data,
33       assumed to arise from a constant sampling rate, is represented by an
34       array reference $data (as described in the methods below), which is
35       then used to create a "Math::FFT" object as
36
37           my $fft = Math::FFT->new($data);
38
39       The methods available include the following.
40
41   FFT METHODS
42       "my $fft = Math::FFT->new($series)"
43           The constructor. Pass it an array of numbers, with a length that is
44           a power of 2.
45
46       "$coeff = $fft->cdft();"
47           This calculates the complex discrete Fourier transform for a data
48           set "x[j]". Here, $data is a reference to an array
49           "data[0...2*n-1]" holding the data
50
51               data[2*j] = Re(x[j]),
52               data[2*j+1] = Im(x[j]), 0<=j<n
53
54           An array reference $coeff is returned consisting of
55
56               coeff[2*k] = Re(X[k]),
57               coeff[2*k+1] = Im(X[k]), 0<=k<n
58
59           where
60
61               X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
62
63       "$orig_data = $fft->invcdft([$coeff]);"
64           Calculates the inverse complex discrete Fourier transform on a data
65           set "x[j]". If $coeff is not given, it will be set equal to an
66           earlier call to "$fft->cdft()". $coeff is a reference to an array
67           "coeff[0...2*n-1]" holding the data
68
69               coeff[2*j] = Re(x[j]),
70               coeff[2*j+1] = Im(x[j]), 0<=j<n
71
72           An array reference $orig_data is returned consisting of
73
74               orig_data[2*k] = Re(X[k]),
75               orig_data[2*k+1] = Im(X[k]), 0<=k<n
76
77           where, excluding the scale,
78
79               X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
80
81           A scaling "$orig_data->[$i] *= 2.0/$n" is then done so that
82           $orig_data coincides with the original $data.
83
84       "$coeff = $fft->rdft();"
85           This calculates the real discrete Fourier transform for a data set
86           "x[j]". On input, $data is a reference to an array "data[0...n-1]"
87           holding the data. An array reference $coeff is returned consisting
88           of
89
90               coeff[2*k] = R[k], 0<=k<n/2
91               coeff[2*k+1] = I[k], 0<k<n/2
92               coeff[1] = R[n/2]
93
94           where
95
96               R[k] = sum_j=0^n-1 data[j]*cos(2*pi*j*k/n), 0<=k<=n/2
97               I[k] = sum_j=0^n-1 data[j]*sin(2*pi*j*k/n), 0<k<n/2
98
99       "$orig_data = $fft->invrdft([$coeff]);"
100           Calculates the inverse real discrete Fourier transform on a data
101           set "coeff[j]". If $coeff is not given, it will be set equal to an
102           earlier call to "$fft->rdft()". $coeff is a reference to an array
103           "coeff[0...n-1]" holding the data
104
105               coeff[2*j] = R[j], 0<=j<n/2
106               coeff[2*j+1] = I[j], 0<j<n/2
107               coeff[1] = R[n/2]
108
109           An array reference $orig_data is returned where, excluding the
110           scale,
111
112               orig_data[k] = (R[0] + R[n/2]*cos(pi*k))/2 +
113                   sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) +
114                       sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
115
116           A scaling "$orig_data->[$i] *= 2.0/$n" is then done so that
117           $orig_data coincides with the original $data.
118
119       "$coeff = $fft->ddct();"
120           Computes the discrete cosine transform on a data set
121           "data[0...n-1]" contained in an array reference $data. An array
122           reference $coeff is returned consisting of
123
124               coeff[k] = C[k], 0<=k<n
125
126           where
127
128               C[k] = sum_j=0^n-1 data[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
129
130       "$orig_data = $fft->invddct([$coeff]);"
131           Computes the inverse discrete cosine transform on a data set
132           "coeff[0...n-1]" contained in an array reference $coeff.  If $coeff
133           is not given, it will be set equal to an earlier call to
134           "$fft->ddct()". An array reference $orig_data is returned
135           consisting of
136
137               orig_data[k] = C[k], 0<=k<n
138
139           where, excluding the scale,
140
141               C[k] = sum_j=0^n-1 coeff[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
142
143           A scaling "$orig_data->[$i] *= 2.0/$n" is then done so that
144           $orig_data coincides with the original $data.
145
146       "$coeff = $fft->ddst();"
147           Computes the discrete sine transform of a data set "data[0...n-1]"
148           contained in an array reference $data. An array reference $coeff is
149           returned consisting of
150
151               coeff[k] = S[k], 0<k<n
152               coeff[0] = S[n]
153
154           where
155
156               S[k] = sum_j=0^n-1 data[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
157
158       "$orig_data = $fft->invddst($coeff);"
159           Computes the inverse discrete sine transform of a data set
160           "coeff[0...n-1]" contained in an array reference $coeff, arranged
161           as
162
163               coeff[j] = A[j], 0<j<n
164               coeff[0] = A[n]
165
166           If $coeff is not given, it will be set equal to an earlier call to
167           "$fft->ddst()". An array reference $orig_data is returned
168           consisting of
169
170               orig_data[k] = S[k], 0<=k<n
171
172           where, excluding a scale,
173
174               S[k] =  sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
175
176           The scaling "$a->[$i] *= 2.0/$n" is then done so that $orig_data
177           coincides with the original $data.
178
179       "$coeff = $fft->dfct();"
180           Computes the real symmetric discrete Fourier transform of a data
181           set "data[0...n]" contained in the array reference $data. An array
182           reference $coeff is returned consisting of
183
184               coeff[k] = C[k], 0<=k<=n
185
186           where
187
188               C[k] = sum_j=0^n data[j]*cos(pi*j*k/n), 0<=k<=n
189
190       "$orig_data = $fft->invdfct($coeff);"
191           Computes the inverse real symmetric discrete Fourier transform of a
192           data set "coeff[0...n]" contained in the array reference $coeff.
193           If $coeff is not given, it will be set equal to an earlier call to
194           "$fft->dfct()". An array reference $orig_data is returned
195           consisting of
196
197               orig_data[k] = C[k], 0<=k<=n
198
199           where, excluding the scale,
200
201               C[k] = sum_j=0^n coeff[j]*cos(pi*j*k/n), 0<=k<=n
202
203           A scaling "$coeff->[0] *= 0.5", "$coeff->[$n] *= 0.5", and
204           "$orig_data->[$i] *= 2.0/$n" is then done so that $orig_data
205           coincides with the original $data.
206
207       "$coeff = $fft->dfst();"
208           Computes the real anti-symmetric discrete Fourier transform of a
209           data set "data[0...n-1]" contained in the array reference $data. An
210           array reference $coeff is returned consisting of
211
212               coeff[k] = C[k], 0<k<n
213
214           where
215
216               C[k] = sum_j=0^n data[j]*sin(pi*j*k/n), 0<k<n
217
218           ("coeff[0]" is used for a work area)
219
220       "$orig_data = $fft->invdfst($coeff);"
221           Computes the inverse real anti-symmetric discrete Fourier transform
222           of a data set "coeff[0...n-1]" contained in the array reference
223           $coeff.  If $coeff is not given, it will be set equal to an earlier
224           call to "$fft->dfst()". An array reference $orig_data is returned
225           consisting of
226
227               orig_data[k] = C[k], 0<k<n
228
229           where, excluding the scale,
230
231               C[k] = sum_j=0^n coeff[j]*sin(pi*j*k/n), 0<k<n
232
233           A scaling "$orig_data->[$i] *= 2.0/$n" is then done so that
234           $orig_data coincides with the original $data.
235
236       "my $other_fft = $fft->clone($other_series)"
237           See "CLONING" below.
238
239       "my $other_series = $fft->convlv($response_data)"
240           See "Convolution" below.
241
242       "my $corr = $fft->correl($other_fft)"
243           See "Correlation" below.
244
245       "my $deconvlv = $fft->deconvlv($respn)"
246           See "Deconvolution" below.
247
248       pdfct()
249           For internal use. Don't use directly.
250
251       pdfst()
252           For internal use. Don't use directly.
253
254       spctrm()
255           See "Power Spectrum" below.
256
257   CLONING
258       The algorithm used in the transforms makes use of arrays for a work
259       area and for a cos/sin lookup table dependent only on the size of the
260       data set. These arrays are initialized when the "Math::FFT" object is
261       created and then are populated when a transform method is first
262       invoked. After this, they persist for the lifetime of the object.
263
264       This aspect is exploited in a "cloning" method; if a "Math::FFT" object
265       is created for a data set $data1 of size "N":
266
267           $fft1 = Math::FFT->new($data1);
268
269       then a new "Math::FFT" object can be created for a second data set
270       $data2 of the same size "N" by
271
272           $fft2 = $fft1->clone($data2);
273
274       The $fft2 object will copy the reuseable work area and lookup table
275       calculated from $fft1.
276
277   APPLICATIONS
278       This module includes some common applications - correlation,
279       convolution and deconvolution, and power spectrum - that arise with
280       real data sets. The conventions used here follow that of Numerical
281       Recipes in C, by Press, Teukolsky, Vetterling, and Flannery, in which
282       further details of the algorithms are given. Note in particular the
283       treatment of end effects by zero padding, which is assumed to be done
284       by the user, if required.
285
286       Correlation
287           The correlation between two functions is defined as
288
289                         /
290               Corr(t) = | ds g(s+t) h(s)
291                         /
292
293           This may be calculated, for two array references $data1 and $data2
294           of the same size $n, as either
295
296               $fft1 = Math::FFT->new($data1);
297               $fft2 = Math::FFT->new($data2);
298               $corr = $fft1->correl($fft2);
299
300           or as
301
302               $fft1 = Math::FFT->new($data1);
303               $corr = $fft1->correl($data2);
304
305           The array reference $corr is returned in wrap-around order -
306           correlations at increasingly positive lags are in "$corr->[0]"
307           (zero lag) on up to "$corr->[$n/2-1]", while correlations at
308           increasingly negative lags are in "$corr->[$n-1]" on down to
309           "$corr->[$n/2]". The sign convention used is such that if $data1
310           lags $data2 (that is, is shifted to the right), then $corr will
311           show a peak at positive lags.
312
313       Convolution
314           The convolution of two functions is defined as
315
316                           /
317               Convlv(t) = | ds g(s) h(t-s)
318                           /
319
320           This is similar to calculating the correlation between the two
321           functions, but typically the functions here have a quite different
322           physical interpretation - one is a signal which persists
323           indefinitely in time, and the other is a response function of
324           limited duration. The convolution may be calculated, for two array
325           references $data and $respn, as
326
327               $fft = Math::FFT->new($data);
328               $convlv = $fft->convlv($respn);
329
330           with the returned $convlv being an array reference. The method
331           assumes that the response function $respn has an odd number of
332           elements $m less than or equal to the number of elements $n of
333           $data. $respn is assumed to be stored in wrap-around order - the
334           first half contains the response at positive times, while the
335           second half, counting down from "$respn->[$m-1]", contains the
336           response at negative times.
337
338       Deconvolution
339           Deconvolution undoes the effects of convoluting a signal with a
340           known response function. In other words, in the relation
341
342                           /
343               Convlv(t) = | ds g(s) h(t-s)
344                           /
345
346           deconvolution reconstructs the original signal, given the
347           convolution and the response function. The method is implemented,
348           for two array references $data and $respn, as
349
350               $fft = Math::FFT->new($data);
351               $deconvlv = $fft->deconvlv($respn);
352
353           As a result, if the convolution of a data set $data with a response
354           function $respn is calculated as
355
356               $fft1 = Math::FFT->new($data);
357               $convlv = $fft1->convlv($respn);
358
359           then the deconvolution
360
361               $fft2 = Math::FFT->new($convlv);
362               $deconvlv = $fft2->deconvlv($respn);
363
364           will give an array reference $deconvlv containing the same elements
365           as the original data $data.
366
367       Power Spectrum
368           If the FFT of a real function of "N" elements is calculated, the
369           "N/2+1" elements of the power spectrum are defined, in terms of the
370           (complex) Fourier coefficients "C[k]", as
371
372               P[0] = |C[0]|^2 / N^2
373               P[k] = 2 |C[k]|^2 / N^2   (k = 1, 2 ,..., N/2-1)
374               P[N/2] = |C[N/2]|^2 / N^2
375
376           Often for these purposes the data is partitioned into "K" segments,
377           each containing "2M" elements. The power spectrum for each segment
378           is calculated, and the net power spectrum is the average of all of
379           these segmented spectra.
380
381           Partitioning may be done in one of two ways: non-overlapping and
382           overlapping. Non-overlapping is useful when the data set is
383           gathered in real time, where the number of data points can be
384           varied at will. Overlapping is useful where there is a fixed number
385           of data points. In non-overlapping, the first <2M> elements
386           constitute segment 1, the next "2M" elements are segment 2, and so
387           on up to segment "K", for a total of "2KM" sampled points. In
388           overlapping, the first and second "M" elements are segment 1, the
389           second and third "M" elements are segment 2, and so on, for a total
390           of "(K+1)M" sampled points.
391
392           A problem that may arise in this procedure is leakage: the power
393           spectrum calculated for one bin contains contributions from nearby
394           bins. To lessen this effect data windowing is often used: multiply
395           the original data "d[j]" by a window function "w[j]", where j = 0,
396           1, ..., N-1. Some popular choices of such functions are
397
398                           | j - N/2 |
399               w[j] = 1 -  | ------- |     ... Bartlett
400                           |   N/2   |
401
402
403                           / j - N/2 \ 2
404               w[j] = 1 -  | ------- |     ... Welch
405                           \   N/2   /
406
407
408                        1   /                    \
409               w[j] =  ---  |1 - cos(2 pi j / N) |     ... Hann
410                        2   \                    /
411
412           The "spctrm" method, used as
413
414               $fft = Math::FFT->new($data);
415               $spectrum = $fft->spctrm(%options);
416
417           returns an array reference $spectrum representing the power
418           spectrum for a data set represented by an array reference $data.
419           The options available are
420
421           "window => window_name"
422               This specifies the window function; if not given, no such
423               function is used. Accepted values (see above) are "bartlett",
424               "welch", "hann", and "\&my_window", where "my_window" is a user
425               specified subroutine which must be of the form, for example,
426
427                   sub my_window {
428                       my ($j, $n) = @_;
429                       return 1 - abs(2*($j-$n/2)/$n);
430                   }
431
432               which implements the Bartlett window.
433
434           "overlap => 1"
435               This specifies whether overlapping should be done; if true (1),
436               overlapping will be used, whereas if false (0), or not
437               specified, no overlapping is used.
438
439           "segments => n"
440               This specifies that the data will be partitioned into "n"
441               segments. If not specified, no segmentation will be done.
442
443           "number => m"
444               This specifies that "2m" data points will be used for each
445               segment, and must be a power of 2. The power spectrum returned
446               will consist of "m+1" elements.
447
448   STATISTICAL FUNCTIONS
449       For convenience, a number of common statistical functions are included
450       for analyzing real data. After creating the object as
451
452           my $fft = Math::FFT->new($data);
453
454       for a data set represented by the array reference $data of size "N",
455       these methods may be called as follows.
456
457       "$mean = $fft->mean([$data]);"
458           This returns the mean
459
460               1/N * sum_j=0^N-1 data[j]
461
462           If an array reference $data is not given, the data set used in
463           creating $fft will be used.
464
465       "$stdev = $fft->stdev([$data]);"
466           This returns the standard deviation
467
468               sqrt{ 1/(N-1) * sum_j=0^N-1 (data[j] - mean)**2 }
469
470           If an array reference $data is not given, the data set used in
471           creating $fft will be used.
472
473       "$rms = $fft->rms([$data]);"
474           This returns the root mean square
475
476               sqrt{ 1/N * sum_j=0^N-1 (data[j])**2 }
477
478           If an array reference $data is not given, the data set used in
479           creating $fft will be used.
480
481       "($min, $max) = $fft->range([$data]);"
482           This returns the minimum and maximum values of the data set.  If an
483           array reference $data is not given, the data set used in creating
484           $fft will be used.
485
486       "$median = $fft->median([$data]);"
487           This returns the median of a data set. The median is defined, for
488           the sorted data set, as either the middle element, if the number of
489           elements is odd, or as the interpolated value of the the two values
490           on either side of the middle, if the number of elements is even. If
491           an array reference $data is not given, the data set used in
492           creating $fft will be used.
493

VERSION

495       version 1.36
496

BUGS

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

SEE ALSO

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

SUPPORT

514   Websites
515       The following websites have more information about this module, and may
516       be of help to you. As always, in addition to those websites please use
517       your favorite search engine to discover more resources.
518
519       •   MetaCPAN
520
521           A modern, open-source CPAN search engine, useful to view POD in
522           HTML format.
523
524           <https://metacpan.org/release/Math-FFT>
525
526       •   RT: CPAN's Bug Tracker
527
528           The RT ( Request Tracker ) website is the default bug/issue
529           tracking system for CPAN.
530
531           <https://rt.cpan.org/Public/Dist/Display.html?Name=Math-FFT>
532
533       •   CPANTS
534
535           The CPANTS is a website that analyzes the Kwalitee ( code metrics )
536           of a distribution.
537
538           <http://cpants.cpanauthors.org/dist/Math-FFT>
539
540       •   CPAN Testers
541
542           The CPAN Testers is a network of smoke testers who run automated
543           tests on uploaded CPAN distributions.
544
545           <http://www.cpantesters.org/distro/M/Math-FFT>
546
547       •   CPAN Testers Matrix
548
549           The CPAN Testers Matrix is a website that provides a visual
550           overview of the test results for a distribution on various
551           Perls/platforms.
552
553           <http://matrix.cpantesters.org/?dist=Math-FFT>
554
555       •   CPAN Testers Dependencies
556
557           The CPAN Testers Dependencies is a website that shows a chart of
558           the test results of all dependencies for a distribution.
559
560           <http://deps.cpantesters.org/?module=Math::FFT>
561
562   Bugs / Feature Requests
563       Please report any bugs or feature requests by email to "bug-math-fft at
564       rt.cpan.org", or through the web interface at
565       <https://rt.cpan.org/Public/Bug/Report.html?Queue=Math-FFT>. You will
566       be automatically notified of any progress on the request by the system.
567
568   Source Code
569       The code is open to the world, and available for you to hack on. Please
570       feel free to browse it and play with it, or whatever. If you want to
571       contribute patches, please send me a diff or prod me to pull from your
572       repository :)
573
574       <https://github.com/shlomif/perl-Math-FFT>
575
576         git clone git://github.com/shlomif/perl-Math-FFT.git
577

AUTHOR

579       Shlomi Fish <shlomif@cpan.org>
580

BUGS

582       Please report any bugs or feature requests on the bugtracker website
583       <https://github.com/shlomif/perl-Math-FFT/issues>
584
585       When submitting a bug or request, please include a test-file or a patch
586       to an existing test-file that illustrates the bug or desired feature.
587
589       This software is copyright (c) 2000 by Randy Kobes.
590
591       This is free software; you can redistribute it and/or modify it under
592       the same terms as the Perl 5 programming language system itself.
593

SUPPORT

595   Websites
596       The following websites have more information about this module, and may
597       be of help to you. As always, in addition to those websites please use
598       your favorite search engine to discover more resources.
599
600       •   MetaCPAN
601
602           A modern, open-source CPAN search engine, useful to view POD in
603           HTML format.
604
605           <https://metacpan.org/release/Math-FFT>
606
607       •   RT: CPAN's Bug Tracker
608
609           The RT ( Request Tracker ) website is the default bug/issue
610           tracking system for CPAN.
611
612           <https://rt.cpan.org/Public/Dist/Display.html?Name=Math-FFT>
613
614       •   CPANTS
615
616           The CPANTS is a website that analyzes the Kwalitee ( code metrics )
617           of a distribution.
618
619           <http://cpants.cpanauthors.org/dist/Math-FFT>
620
621       •   CPAN Testers
622
623           The CPAN Testers is a network of smoke testers who run automated
624           tests on uploaded CPAN distributions.
625
626           <http://www.cpantesters.org/distro/M/Math-FFT>
627
628       •   CPAN Testers Matrix
629
630           The CPAN Testers Matrix is a website that provides a visual
631           overview of the test results for a distribution on various
632           Perls/platforms.
633
634           <http://matrix.cpantesters.org/?dist=Math-FFT>
635
636       •   CPAN Testers Dependencies
637
638           The CPAN Testers Dependencies is a website that shows a chart of
639           the test results of all dependencies for a distribution.
640
641           <http://deps.cpantesters.org/?module=Math::FFT>
642
643   Bugs / Feature Requests
644       Please report any bugs or feature requests by email to "bug-math-fft at
645       rt.cpan.org", or through the web interface at
646       <https://rt.cpan.org/Public/Bug/Report.html?Queue=Math-FFT>. You will
647       be automatically notified of any progress on the request by the system.
648
649   Source Code
650       The code is open to the world, and available for you to hack on. Please
651       feel free to browse it and play with it, or whatever. If you want to
652       contribute patches, please send me a diff or prod me to pull from your
653       repository :)
654
655       <https://github.com/shlomif/perl-Math-FFT>
656
657         git clone git://github.com/shlomif/perl-Math-FFT.git
658

AUTHOR

660       Shlomi Fish <shlomif@cpan.org>
661

BUGS

663       Please report any bugs or feature requests on the bugtracker website
664       <https://github.com/shlomif/perl-Math-FFT/issues>
665
666       When submitting a bug or request, please include a test-file or a patch
667       to an existing test-file that illustrates the bug or desired feature.
668
670       This software is copyright (c) 2000 by Randy Kobes.
671
672       This is free software; you can redistribute it and/or modify it under
673       the same terms as the Perl 5 programming language system itself.
674
675
676
677perl v5.38.0                      2023-07-20                      Math::FFT(3)
Impressum