1FFT(3) User Contributed Perl Documentation FFT(3)
2
6 Math::FFT - Perl module to calculate Fast Fourier Transforms
7
9 use Math::FFT;
10 my $PI = 3.1415926539;
11 my $N = 64;
12 my ($series, $other_series);
13 for (my $k=0; $k<$N; $k++) {
14 $series->[$k] = sin(4*$k*$PI/$N) + cos(6*$k*$PI/$N);
15 }
16 my $fft = new Math::FFT($series);
17 my $coeff = $fft->rdft();
18 my $spectrum = $fft->spctrm;
19 my $original_data = $fft->invrdft($coeff);
20 for (my $k=0; $k<$N; $k++) {
22 $other_series->[$k] = sin(16*$k*$PI/$N) + cos(8*$k*$PI/$N);
23 }
24 my $other_fft = $fft->clone($other_series);
25 my $other_coeff = $other_fft->rdft();
26 my $correlation = $fft->correl($other_fft);
27
29 This module implements some algorithms for calculating Fast Fourier
30 Transforms for one-dimensional data sets of size 2^n. The data,
31 assumed to arise from a constant sampling rate, is represented by an
32 array reference $data (as described in the methods below), which is
33 then used to create a "Math::FFT" object as
34 my $fft = new Math::FFT($data);
36 The methods available include the following.
38 FFT METHODS
40 "$coeff = $fft->cdft();"
41 This calculates the complex discrete Fourier transform for a data
42 set "x[j]". Here, $data is a reference to an array
43 "data[0...2*n-1]" holding the data
44 data[2*j] = Re(x[j]),
46 data[2*j+1] = Im(x[j]), 0<=j<n
47 An array reference $coeff is returned consisting of
49 coeff[2*k] = Re(X[k]),
51 coeff[2*k+1] = Im(X[k]), 0<=k<n
52 where
54 X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
56 "$orig_data = $fft->invcdft([$coeff]);"
58 Calculates the inverse complex discrete Fourier transform on a data
59 set "x[j]". If $coeff is not given, it will be set equal to an
60 earlier call to "$fft->cdft()". $coeff is a reference to an array
61 "coeff[0...2*n-1]" holding the data
62 coeff[2*j] = Re(x[j]),
64 coeff[2*j+1] = Im(x[j]), 0<=j<n
65 An array reference $orig_data is returned consisting of
67 orig_data[2*k] = Re(X[k]),
69 orig_data[2*k+1] = Im(X[k]), 0<=k<n
70 where, excluding the scale,
72 X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
74 A scaling "$orig_data->[$i] *= 2.0/$n" is then done so that
76 $orig_data coincides with the original $data.
77 "$coeff = $fft->rdft();"
79 This calculates the real discrete Fourier transform for a data set
80 "x[j]". On input, $data is a reference to an array "data[0...n-1]"
81 holding the data. An array reference $coeff is returned consisting
82 of
83 coeff[2*k] = R[k], 0<=k<n/2
85 coeff[2*k+1] = I[k], 0<k<n/2
86 coeff[1] = R[n/2]
87 where
89 R[k] = sum_j=0^n-1 data[j]*cos(2*pi*j*k/n), 0<=k<=n/2
91 I[k] = sum_j=0^n-1 data[j]*sin(2*pi*j*k/n), 0<k<n/2
92 "$orig_data = $fft->invrdft([$coeff]);"
94 Calculates the inverse real discrete Fourier transform on a data
95 set "coeff[j]". If $coeff is not given, it will be set equal to an
96 earlier call to "$fft->rdft()". $coeff is a reference to an array
97 "coeff[0...n-1]" holding the data
98 coeff[2*j] = R[j], 0<=j<n/2
100 coeff[2*j+1] = I[j], 0<j<n/2
101 coeff[1] = R[n/2]
102 An array reference $orig_data is returned where, excluding the
104 scale,
105 orig_data[k] = (R[0] + R[n/2]*cos(pi*k))/2 +
107 sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) +
108 sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
109 A scaling "$orig_data->[$i] *= 2.0/$n" is then done so that
111 $orig_data coincides with the original $data.
112 "$coeff = $fft->ddct();"
114 Computes the discrete cosine tranform on a data set "data[0...n-1]"
115 contained in an array reference $data. An array reference $coeff is
116 returned consisting of
117 coeff[k] = C[k], 0<=k<n
119 where
121 C[k] = sum_j=0^n-1 data[j]*cos(pi*(j+1/2)*k/n), 0<=k<n
123 "$orig_data = $fft->invddct([$coeff]);"
125 Computes the inverse discrete cosine tranform on a data set
126 "coeff[0...n-1]" contained in an array reference $coeff. If $coeff
127 is not given, it will be set equal to an earlier call to
128 "$fft->ddct()". An array reference $orig_data is returned
129 consisting of
130 orig_data[k] = C[k], 0<=k<n
132 where, excluding the scale,
134 C[k] = sum_j=0^n-1 coeff[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
136 A scaling "$orig_data->[$i] *= 2.0/$n" is then done so that
138 $orig_data coincides with the original $data.
139 "$coeff = $fft->ddst();"
141 Computes the discrete sine transform of a data set "data[0...n-1]"
142 contained in an array reference $data. An array reference $coeff is
143 returned consisting of
144 coeff[k] = S[k], 0<k<n
146 coeff[0] = S[n]
147 where
149 S[k] = sum_j=0^n-1 data[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
151 "$orig_data = $fft->invddst($coeff);"
153 Computes the inverse discrete sine transform of a data set
154 "coeff[0...n-1]" contained in an array reference $coeff, arranged
155 as
156 coeff[j] = A[j], 0<j<n
158 coeff[0] = A[n]
159 If $coeff is not given, it will be set equal to an earlier call to
161 "$fft->ddst()". An array reference $orig_data is returned
162 consisting of
163 orig_data[k] = S[k], 0<=k<n
165 where, excluding a scale,
167 S[k] = sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
169 The scaling "$a->[$i] *= 2.0/$n" is then done so that $orig_data
171 coincides with the original $data.
172 "$coeff = $fft->dfct();"
174 Computes the real symmetric discrete Fourier transform of a data
175 set "data[0...n]" contained in the array reference $data. An array
176 reference $coeff is returned consisting of
177 coeff[k] = C[k], 0<=k<=n
179 where
181 C[k] = sum_j=0^n data[j]*cos(pi*j*k/n), 0<=k<=n
183 "$orig_data = $fft->invdfct($coeff);"
185 Computes the inverse real symmetric discrete Fourier transform of a
186 data set "coeff[0...n]" contained in the array reference $coeff.
187 If $coeff is not given, it will be set equal to an earlier call to
188 "$fft->dfct()". An array reference $orig_data is returned
189 consisting of
190 orig_data[k] = C[k], 0<=k<=n
192 where, excluding the scale,
194 C[k] = sum_j=0^n coeff[j]*cos(pi*j*k/n), 0<=k<=n
196 A scaling "$coeff->[0] *= 0.5", "$coeff->[$n] *= 0.5", and
198 "$orig_data->[$i] *= 2.0/$n" is then done so that $orig_data
199 coincides with the original $data.
200 "$coeff = $fft->dfst();"
202 Computes the real anti-symmetric discrete Fourier transform of a
203 data set "data[0...n-1]" contained in the array reference $data. An
204 array reference $coeff is returned consisting of
205 coeff[k] = C[k], 0<k<n
207 where
209 C[k] = sum_j=0^n data[j]*sin(pi*j*k/n), 0<k<n
211 ("coeff[0]" is used for a work area)
213 "$orig_data = $fft->invdfst($coeff);"
215 Computes the inverse real anti-symmetric discrete Fourier transform
216 of a data set "coeff[0...n-1]" contained in the array reference
217 $coeff. If $coeff is not given, it will be set equal to an earlier
218 call to "$fft->dfst()". An array reference $orig_data is returned
219 consisting of
220 orig_data[k] = C[k], 0<k<n
222 where, excluding the scale,
224 C[k] = sum_j=0^n coeff[j]*sin(pi*j*k/n), 0<k<n
226 A scaling "$orig_data->[$i] *= 2.0/$n" is then done so that
228 $orig_data coincides with the original $data.
229 CLONING
231 The algorithm used in the transforms makes use of arrays for a work
232 area and for a cos/sin lookup table dependent only on the size of the
233 data set. These arrays are initialized when the "Math::FFT" object is
234 created and then are populated when a transform method is first
235 invoked. After this, they persist for the lifetime of the object.
236 This aspect is exploited in a "cloning" method; if a "Math::FFT" object
238 is created for a data set $data1 of size "N":
239 $fft1 = new Math::FFT($data1);
241 then a new "Math::FFT" object can be created for a second data set
243 $data2 of the same size "N" by
244 $fft2 = $fft1->clone($data2);
246 The $fft2 object will copy the reuseable work area and lookup table
248 calculated from $fft1.
249 APPLICATIONS
251 This module includes some common applications - correlation,
252 convolution and deconvolution, and power spectrum - that arise with
253 real data sets. The conventions used here follow that of Numerical
254 Recipes in C, by Press, Teukolsky, Vetterling, and Flannery, in which
255 further details of the algorithms are given. Note in particular the
256 treatment of end effects by zero padding, which is assumed to be done
257 by the user, if required.
258 Correlation
260 The correlation between two functions is defined as
261 /
263 Corr(t) = | ds g(s+t) h(s)
264 /
265 This may be calculated, for two array references $data1 and $data2
267 of the same size $n, as either
268 $fft1 = new Math::FFT($data1);
270 $fft2 = new Math::FFT($data2);
271 $corr = $fft1->correl($fft2);
272 or as
274 $fft1 = new Math::FFT($data1);
276 $corr = $fft1->correl($data2);
277 The array reference $corr is returned in wrap-around order -
279 correlations at increasingly positive lags are in "$corr->[0]"
280 (zero lag) on up to "$corr->[$n/2-1]", while correlations at
281 increasingly negative lags are in "$corr->[$n-1]" on down to
282 "$corr->[$n/2]". The sign convention used is such that if $data1
283 lags $data2 (that is, is shifted to the right), then $corr will
284 show a peak at positive lags.
285 Convolution
287 The convolution of two functions is defined as
288 /
290 Convlv(t) = | ds g(s) h(t-s)
291 /
292 This is similar to calculating the correlation between the two
294 functions, but typically the functions here have a quite different
295 physical interpretation - one is a signal which persists
296 indefinitely in time, and the other is a response function of
297 limited duration. The convolution may be calculated, for two array
298 references $data and $respn, as
299 $fft = new Math::FFT($data);
301 $convlv = $fft->convlv($respn);
302 with the returned $convlv being an array reference. The method
304 assumes that the response function $respn has an odd number of
305 elements $m less than or equal to the number of elements $n of
306 $data. $respn is assumed to be stored in wrap-around order - the
307 first half contains the response at positive times, while the
308 second half, counting down from "$respn->[$m-1]", contains the
309 response at negative times.
310 Deconvolution
312 Deconvolution undoes the effects of convoluting a signal with a
313 known response function. In other words, in the relation
314 /
316 Convlv(t) = | ds g(s) h(t-s)
317 /
318 deconvolution reconstructs the original signal, given the
320 convolution and the response function. The method is implemented,
321 for two array references $data and $respn, as
322 $fft = new Math::FFT($data);
324 $deconvlv = $fft->deconvlv($respn);
325 As a result, if the convolution of a data set $data with a response
327 function $respn is calculated as
328 $fft1 = new Math::FFT($data);
330 $convlv = $fft1->convlv($respn);
331 then the deconvolution
333 $fft2 = new Math::FFT($convlv);
335 $deconvlv = $fft2->deconvlv($respn);
336 will give an array reference $deconvlv containing the same elements
338 as the original data $data.
339 Power Spectrum
341 If the FFT of a real function of "N" elements is calculated, the
342 "N/2+1" elements of the power spectrum are defined, in terms of the
343 (complex) Fourier coefficients "C[k]", as
344 P[0] = |C[0]|^2 / N^2
346 P[k] = 2 |C[k]|^2 / N^2 (k = 1, 2 ,..., N/2-1)
347 P[N/2] = |C[N/2]|^2 / N^2
348 Often for these purposes the data is partitioned into "K" segments,
350 each containing "2M" elements. The power spectrum for each segment
351 is calculated, and the net power spectrum is the average of all of
352 these segmented spectra.
353 Partitioning may be done in one of two ways: non-overlapping and
355 overlapping. Non-overlapping is useful when the data set is
356 gathered in real time, where the number of data points can be
357 varied at will. Overlapping is useful where there is a fixed number
358 of data points. In non-overlapping, the first <2M> elements
359 constitute segment 1, the next "2M" elements are segment 2, and so
360 on up to segment "K", for a total of "2KM" sampled points. In
361 overlapping, the first and second "M" elements are segment 1, the
362 second and third "M" elements are segment 2, and so on, for a total
363 of "(K+1)M" sampled points.
364 A problem that may arise in this procedure is leakage: the power
366 spectrum calculated for one bin contains contributions from nearby
367 bins. To lessen this effect data windowing is often used: multiply
368 the original data "d[j]" by a window function "w[j]", where j = 0,
369 1, ..., N-1. Some popular choices of such functions are
370 | j - N/2 |
372 w[j] = 1 - | ------- | ... Bartlett
373 | N/2 |
374 / j - N/2 \ 2
377 w[j] = 1 - | ------- | ... Welch
378 \ N/2 /
379 1 / \
382 w[j] = --- |1 - cos(2 pi j / N) | ... Hann
383 2 \ /
384 The "spctrm" method, used as
386 $fft = Math::FFT->new($data);
388 $spectrum = $fft->spctrm(%options);
389 returns an array reference $spectrum representing the power
391 spectrum for a data set represented by an array reference $data.
392 The options available are
393 "window => window_name"
395 This specifies the window function; if not given, no such
396 function is used. Accepted values (see above) are "bartlett",
397 "welch", "hann", and "\&my_window", where "my_window" is a user
398 specified subroutine which must be of the form, for example,
399 sub my_window {
401 my ($j, $n) = @_;
402 return 1 - abs(2*($j-$n/2)/$n);
403 }
404 which implements the Bartlett window.
406 "overlap => 1"
408 This specifies whether overlapping should be done; if true (1),
409 overlapping will be used, whereas if false (0), or not
410 specified, no overlapping is used.
411 "segments => n"
413 This specifies that the data will be partitioned into "n"
414 segments. If not specified, no segmentation will be done.
415 "number => m"
417 This specifies that "2m" data points will be used for each
418 segment, and must be a power of 2. The power spectrum returned
419 will consist of "m+1" elements.
420 STATISTICAL FUNCTIONS
422 For convenience, a number of common statistical functions are included
423 for analyzing real data. After creating the object as
424 my $fft = new Math::FFT($data);
426 for a data set represented by the array reference $data of size "N",
428 these methods may be called as follows.
429 "$mean = $fft->mean([$data]);"
431 This returns the mean
432 1/N * sum_j=0^N-1 data[j]
434 If an array reference $data is not given, the data set used in
436 creating $fft will be used.
437 "$stdev = $fft->stdev([$data]);"
439 This returns the standard deviation
440 sqrt{ 1/(N-1) * sum_j=0^N-1 (data[j] - mean)**2 }
442 If an array reference $data is not given, the data set used in
444 creating $fft will be used.
445 "$rms = $fft->rms([$data]);"
447 This returns the root mean square
448 sqrt{ 1/N * sum_j=0^N-1 (data[j])**2 }
450 If an array reference $data is not given, the data set used in
452 creating $fft will be used.
453 "($min, $max) = $fft->range([$data]);"
455 This returns the minimum and maximum values of the data set. If an
456 array reference $data is not given, the data set used in creating
457 $fft will be used.
458 "$median = $fft->median([$data]);"
460 This returns the median of a data set. The median is defined, for
461 the sorted data set, as either the middle element, if the number of
462 elements is odd, or as the interpolated value of the the two values
463 on either side of the middle, if the number of elements is even. If
464 an array reference $data is not given, the data set used in
465 creating $fft will be used.
466
468 Please report any to Randy Kobes <randy@theoryx5.uwinnipeg.ca>
469
471 Math::Pari and PDL
472
474 The algorithm used in this module to calculate the Fourier transforms
475 is based on the C routine of fft4g.c available at
476 http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html, which is copyrighted
477 1996-99 by Takuya OOURA. The file arrays.c included here to handle
478 passing arrays to and from C comes from the PGPLOT module of Karl
479 Glazebrook <kgb@aaoepp.aao.gov.au>. The perl code of Math::FFT is
480 copyright 2000,2005 by Randy Kobes <r.kobes@uwinnipeg.ca>, and is
481 distributed under the same terms as Perl itself.
482perl v5.12.0 2005-09-03 FFT(3)