1Math::FFT(3) User Contributed Perl Documentation Math::FFT(3)
2
6 Math::FFT - Perl module to calculate Fast Fourier Transforms
7
9 version 1.34
10
12 use Math::FFT;
13 my $PI = 3.1415926539;
15 my $N = 64;
16 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 my $other_series =
24 [map { sin(16*$_*$PI/$N) + cos(8*$_*$PI/$N) } 0 .. $N-1];
25 my $other_fft = $fft->clone($other_series);
27 my $other_coeff = $other_fft->rdft();
28 my $correlation = $fft->correl($other_fft);
29
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 my $fft = Math::FFT->new($data);
38 The methods available include the following.
40 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 "$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 data[2*j] = Re(x[j]),
52 data[2*j+1] = Im(x[j]), 0<=j<n
53 An array reference $coeff is returned consisting of
55 coeff[2*k] = Re(X[k]),
57 coeff[2*k+1] = Im(X[k]), 0<=k<n
58 where
60 X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
62 "$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 coeff[2*j] = Re(x[j]),
70 coeff[2*j+1] = Im(x[j]), 0<=j<n
71 An array reference $orig_data is returned consisting of
73 orig_data[2*k] = Re(X[k]),
75 orig_data[2*k+1] = Im(X[k]), 0<=k<n
76 where, excluding the scale,
78 X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
80 A scaling "$orig_data->[$i] *= 2.0/$n" is then done so that
82 $orig_data coincides with the original $data.
83 "$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 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 where
95 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 "$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 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 An array reference $orig_data is returned where, excluding the
110 scale,
111 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 A scaling "$orig_data->[$i] *= 2.0/$n" is then done so that
117 $orig_data coincides with the original $data.
118 "$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 coeff[k] = C[k], 0<=k<n
125 where
127 C[k] = sum_j=0^n-1 data[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
129 "$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 orig_data[k] = C[k], 0<=k<n
138 where, excluding the scale,
140 C[k] = sum_j=0^n-1 coeff[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
142 A scaling "$orig_data->[$i] *= 2.0/$n" is then done so that
144 $orig_data coincides with the original $data.
145 "$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 coeff[k] = S[k], 0<k<n
152 coeff[0] = S[n]
153 where
155 S[k] = sum_j=0^n-1 data[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
157 "$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 coeff[j] = A[j], 0<j<n
164 coeff[0] = A[n]
165 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 orig_data[k] = S[k], 0<=k<n
171 where, excluding a scale,
173 S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
175 The scaling "$a->[$i] *= 2.0/$n" is then done so that $orig_data
177 coincides with the original $data.
178 "$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 coeff[k] = C[k], 0<=k<=n
185 where
187 C[k] = sum_j=0^n data[j]*cos(pi*j*k/n), 0<=k<=n
189 "$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 orig_data[k] = C[k], 0<=k<=n
198 where, excluding the scale,
200 C[k] = sum_j=0^n coeff[j]*cos(pi*j*k/n), 0<=k<=n
202 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 "$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 coeff[k] = C[k], 0<k<n
213 where
215 C[k] = sum_j=0^n data[j]*sin(pi*j*k/n), 0<k<n
217 ("coeff[0]" is used for a work area)
219 "$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 orig_data[k] = C[k], 0<k<n
228 where, excluding the scale,
230 C[k] = sum_j=0^n coeff[j]*sin(pi*j*k/n), 0<k<n
232 A scaling "$orig_data->[$i] *= 2.0/$n" is then done so that
234 $orig_data coincides with the original $data.
235 "my $other_fft = $fft->clone($other_series)"
237 See "CLONING" below.
238 "my $other_series = $fft->convlv($response_data)"
240 See "Convolution" below.
241 "my $corr = $fft->correl($other_fft)"
243 See "Correlation" below.
244 "my $deconvlv = $fft->deconvlv($respn)"
246 See "Deconvolution" below.
247 pdfct()
249 For internal use. Don't use directly.
250 pdfst()
252 For internal use. Don't use directly.
253 spctrm()
255 See "Power Spectrum" below.
256 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 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 $fft1 = Math::FFT->new($data1);
268 then a new "Math::FFT" object can be created for a second data set
270 $data2 of the same size "N" by
271 $fft2 = $fft1->clone($data2);
273 The $fft2 object will copy the reuseable work area and lookup table
275 calculated from $fft1.
276 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 Correlation
287 The correlation between two functions is defined as
288 /
290 Corr(t) = | ds g(s+t) h(s)
291 /
292 This may be calculated, for two array references $data1 and $data2
294 of the same size $n, as either
295 $fft1 = Math::FFT->new($data1);
297 $fft2 = Math::FFT->new($data2);
298 $corr = $fft1->correl($fft2);
299 or as
301 $fft1 = Math::FFT->new($data1);
303 $corr = $fft1->correl($data2);
304 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 Convolution
314 The convolution of two functions is defined as
315 /
317 Convlv(t) = | ds g(s) h(t-s)
318 /
319 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 $fft = Math::FFT->new($data);
328 $convlv = $fft->convlv($respn);
329 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 Deconvolution
339 Deconvolution undoes the effects of convoluting a signal with a
340 known response function. In other words, in the relation
341 /
343 Convlv(t) = | ds g(s) h(t-s)
344 /
345 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 $fft = Math::FFT->new($data);
351 $deconvlv = $fft->deconvlv($respn);
352 As a result, if the convolution of a data set $data with a response
354 function $respn is calculated as
355 $fft1 = Math::FFT->new($data);
357 $convlv = $fft1->convlv($respn);
358 then the deconvolution
360 $fft2 = Math::FFT->new($convlv);
362 $deconvlv = $fft2->deconvlv($respn);
363 will give an array reference $deconvlv containing the same elements
365 as the original data $data.
366 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 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 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 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 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 | j - N/2 |
399 w[j] = 1 - | ------- | ... Bartlett
400 | N/2 |
401 / j - N/2 \ 2
404 w[j] = 1 - | ------- | ... Welch
405 \ N/2 /
406 1 / \
409 w[j] = --- |1 - cos(2 pi j / N) | ... Hann
410 2 \ /
411 The "spctrm" method, used as
413 $fft = Math::FFT->new($data);
415 $spectrum = $fft->spctrm(%options);
416 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 "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 sub my_window {
428 my ($j, $n) = @_;
429 return 1 - abs(2*($j-$n/2)/$n);
430 }
431 which implements the Bartlett window.
433 "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 "segments => n"
440 This specifies that the data will be partitioned into "n"
441 segments. If not specified, no segmentation will be done.
442 "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 STATISTICAL FUNCTIONS
449 For convenience, a number of common statistical functions are included
450 for analyzing real data. After creating the object as
451 my $fft = Math::FFT->new($data);
453 for a data set represented by the array reference $data of size "N",
455 these methods may be called as follows.
456 "$mean = $fft->mean([$data]);"
458 This returns the mean
459 1/N * sum_j=0^N-1 data[j]
461 If an array reference $data is not given, the data set used in
463 creating $fft will be used.
464 "$stdev = $fft->stdev([$data]);"
466 This returns the standard deviation
467 sqrt{ 1/(N-1) * sum_j=0^N-1 (data[j] - mean)**2 }
469 If an array reference $data is not given, the data set used in
471 creating $fft will be used.
472 "$rms = $fft->rms([$data]);"
474 This returns the root mean square
475 sqrt{ 1/N * sum_j=0^N-1 (data[j])**2 }
477 If an array reference $data is not given, the data set used in
479 creating $fft will be used.
480 "($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 "$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
495 Please report any to Randy Kobes <randy@theoryx5.uwinnipeg.ca>
496
498 Math::Pari and PDL
499
501 The algorithm used in this module to calculate the Fourier transforms
502 is based on the C routine of fft4g.c available at
503 http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html, which is copyrighted
504 1996-99 by Takuya OOURA. The file arrays.c included here to handle
505 passing arrays to and from C comes from the PGPLOT module of Karl
506 Glazebrook <kgb@aaoepp.aao.gov.au>. The perl code of Math::FFT is
507 copyright 2000,2005 by Randy Kobes <r.kobes@uwinnipeg.ca>, and is
508 distributed under the same terms as Perl itself.
509
511 Shlomi Fish <shlomif@cpan.org>
512
514 This software is copyright (c) 2000 by Randy Kobes.
515 This is free software; you can redistribute it and/or modify it under
517 the same terms as the Perl 5 programming language system itself.
518
520 Please report any bugs or feature requests on the bugtracker website
521 <https://rt.cpan.org/Public/Dist/Display.html?Name=Math-FFT> or by
522 email to bug-math-fft@rt.cpan.org <mailto:bug-math-fft@rt.cpan.org>.
523 When submitting a bug or request, please include a test-file or a patch
525 to an existing test-file that illustrates the bug or desired feature.
526
528 Perldoc
529 You can find documentation for this module with the perldoc command.
530 perldoc Math::FFT
532 Websites
534 The following websites have more information about this module, and may
535 be of help to you. As always, in addition to those websites please use
536 your favorite search engine to discover more resources.
537 · MetaCPAN
539 A modern, open-source CPAN search engine, useful to view POD in
541 HTML format.
542 <http://metacpan.org/release/Math-FFT>
544 · Search CPAN
546 The default CPAN search engine, useful to view POD in HTML format.
548 <http://search.cpan.org/dist/Math-FFT>
550 · RT: CPAN's Bug Tracker
552 The RT ( Request Tracker ) website is the default bug/issue
554 tracking system for CPAN.
555 <https://rt.cpan.org/Public/Dist/Display.html?Name=Math-FFT>
557 · AnnoCPAN
559 The AnnoCPAN is a website that allows community annotations of Perl
561 module documentation.
562 <http://annocpan.org/dist/Math-FFT>
564 · CPAN Ratings
566 The CPAN Ratings is a website that allows community ratings and
568 reviews of Perl modules.
569 <http://cpanratings.perl.org/d/Math-FFT>
571 · CPAN Forum
573 The CPAN Forum is a web forum for discussing Perl modules.
575 <http://cpanforum.com/dist/Math-FFT>
577 · CPANTS
579 The CPANTS is a website that analyzes the Kwalitee ( code metrics )
581 of a distribution.
582 <http://cpants.cpanauthors.org/dist/Math-FFT>
584 · CPAN Testers
586 The CPAN Testers is a network of smokers who run automated tests on
588 uploaded CPAN distributions.
589 <http://www.cpantesters.org/distro/M/Math-FFT>
591 · CPAN Testers Matrix
593 The CPAN Testers Matrix is a website that provides a visual
595 overview of the test results for a distribution on various
596 Perls/platforms.
597 <http://matrix.cpantesters.org/?dist=Math-FFT>
599 · CPAN Testers Dependencies
601 The CPAN Testers Dependencies is a website that shows a chart of
603 the test results of all dependencies for a distribution.
604 <http://deps.cpantesters.org/?module=Math::FFT>
606 Bugs / Feature Requests
608 Please report any bugs or feature requests by email to "bug-math-fft at
609 rt.cpan.org", or through the web interface at
610 <https://rt.cpan.org/Public/Bug/Report.html?Queue=Math-FFT>. You will
611 be automatically notified of any progress on the request by the system.
612 Source Code
614 The code is open to the world, and available for you to hack on. Please
615 feel free to browse it and play with it, or whatever. If you want to
616 contribute patches, please send me a diff or prod me to pull from your
617 repository :)
618 <https://github.com/shlomif/perl-Math-FFT>
620 git clone https://github.com/shlomif/perl-Math-FFT.git
622perl v5.28.0 2017-04-07 Math::FFT(3)