home *** CD-ROM | disk | FTP | other *** search
/ PC Professionell 2004 December / PCpro_2004_12.ISO / files / webserver / tsw / TSW_3.4.0.exe / Apache2 / perl / FFT.pm < prev    next >
Encoding:
Perl POD Document  |  2001-01-06  |  29.2 KB  |  1,079 lines
  1. package Math::FFT;
  2.  
  3. use strict;
  4. use vars qw($VERSION @ISA);
  5. require DynaLoader;
  6.  
  7. @ISA = qw(DynaLoader);
  8.  
  9. # Items to export into callers namespace by default. Note: do not export
  10. # names by default without a very good reason. Use EXPORT_OK instead.
  11. # Do not simply export all your public functions/methods/constants.
  12.  
  13. $VERSION = '0.25';
  14.  
  15. bootstrap Math::FFT $VERSION;
  16.  
  17. # Preloaded methods go here.
  18.  
  19. sub new {
  20.   my ($class, $data) = @_;
  21.   die 'Must call constructor with an array reference for the data'
  22.     unless ref($data) eq 'ARRAY';
  23.   $data->[0] ||= 0;    # keep warnings happy
  24.   my $n = @$data;
  25.   my $nip = int(3 + sqrt($n));
  26.   my $nw = int(2 + 5*$n/4);
  27.   my $ip = pack("i$nip", ());
  28.   my $w = pack("d$nw", ());
  29.   bless {
  30.      type => '',
  31.      mean => '',
  32.      coeff => '',
  33.      n   => $n,
  34.      data => $data,
  35.          ip => \$ip,
  36.          w => \$w,
  37.     }, $class;
  38. }
  39.  
  40. # clone method to copy the ip and w arrays for data of equal size
  41. sub clone {
  42.   my ($self, $data) = @_;
  43.   die 'Must call clone with an array reference for the data'
  44.     unless ref($data) eq 'ARRAY';
  45.   $data->[0] ||= 0;    # keep warnings happy
  46.   my $n = @$data;
  47.   die "Cannot clone data of unequal sizes" unless $n == $self->{n};
  48.   my $class = ref($self);
  49.   bless {
  50.      type => '',
  51.      coeff => '',
  52.      mean => '',
  53.      n   => $self->{n},
  54.      data => $data,
  55.      ip => $self->{ip},
  56.      w => $self->{w},
  57.     }, $class;
  58. }
  59.  
  60. # Complex Discrete Fourier Transform
  61. sub cdft {
  62.   my $self = shift;
  63.   my $n = $self->{n};
  64.   die 'size of data must be an integer power of 2' unless check_n($n);
  65.   my $data = [ @{$self->{data}} ];
  66.   _cdft($n, 1, $data, $self->{ip}, $self->{w});
  67.   $self->{type} = 'cdft';
  68.   $self->{coeff} = $data;
  69.   return $data;
  70. }
  71.  
  72. # Inverse Complex Discrete Fourier Transform
  73. sub invcdft {
  74.   my $self = shift;
  75.   my $data;
  76.   my $n = $self->{n};
  77.   if (my $arg = shift) {
  78.     die 'Must pass an array reference to invcdft'
  79.       unless ref($arg) eq 'ARRAY';
  80.     die "Size of data set must be $n"
  81.       unless $n == @$arg;
  82.     $data = [ @$arg ];
  83.   }
  84.   else {
  85.     die 'Must invert data created with cdft'
  86.       unless $self->{type} eq 'cdft';
  87.     $data =  [ @{$self->{coeff}} ];
  88.   }
  89.   _cdft($n, -1, $data, $self->{ip}, $self->{w});
  90.   $_ *= 2.0/$n for (@$data);
  91.   return $data;
  92. }
  93.  
  94. # Real Discrete Fourier Transform
  95. sub rdft {
  96.   my $self = shift;
  97.   my $n = $self->{n};
  98.   die 'size of data must be an integer power of 2' unless check_n($n);
  99.   my $data = [ @{$self->{data}} ];
  100.   _rdft($n, 1, $data, $self->{ip}, $self->{w});
  101.   $self->{type} = 'rdft';
  102.   $self->{coeff} = $data;
  103.   return $data;
  104. }
  105.  
  106. # Inverse Real Discrete Fourier Transform
  107. sub invrdft {
  108.   my $self = shift;
  109.   my $data;
  110.   my $n = $self->{n};
  111.   if (my $arg = shift) {
  112.     die 'Must pass an array reference to invrdft'
  113.       unless ref($arg) eq 'ARRAY';
  114.     die "Size of data set must be $n"
  115.       unless $n == @$arg;
  116.     $data = [ @$arg ];
  117.   }
  118.   else {
  119.     die 'Must invert data created with rdft'
  120.       unless $self->{type} eq 'rdft';
  121.     $data =  [ @{$self->{coeff}} ];
  122.   }
  123.   _rdft($n, -1, $data, $self->{ip}, $self->{w});
  124.   $_ *= 2.0/$n for (@$data);
  125.   return $data;
  126. }
  127.  
  128. # Discrete Cosine Transform
  129. sub ddct {
  130.   my $self = shift;
  131.   my $n = $self->{n};
  132.   die 'size of data must be an integer power of 2' unless check_n($n);
  133.   my $data = [ @{$self->{data}} ];
  134.   _ddct($n, -1, $data, $self->{ip}, $self->{w});
  135.   $self->{type} = 'ddct';
  136.   $self->{coeff} = $data;
  137.   return $data;
  138. }
  139.  
  140. # Inverse Discrete Cosine Transform
  141. sub invddct {
  142.   my $self = shift;
  143.   my $data;
  144.   my $n = $self->{n};
  145.   if (my $arg = shift) {
  146.     die 'Must pass an array reference to invddct'
  147.       unless ref($arg) eq 'ARRAY';
  148.     die "Size of data set must be $n"
  149.       unless $n == @$arg;
  150.     $data = [ @$arg ];
  151.   }
  152.   else {
  153.     die 'Must invert data created with ddct'
  154.       unless $self->{type} eq 'ddct';
  155.     $data =  [ @{$self->{coeff}} ];
  156.   }
  157.   $data->[0] *= 0.5;
  158.   _ddct($n, 1, $data, $self->{ip}, $self->{w});
  159.   $_ *= 2.0/$n for (@$data);
  160.   return $data;
  161. }
  162.  
  163. # Discrete Sine Transform
  164. sub ddst {
  165.   my $self = shift;
  166.   my $n = $self->{n};
  167.   die 'size of data must be an integer power of 2' unless check_n($n);
  168.   my $data = [ @{$self->{data}} ];
  169.   _ddst($n, -1, $data, $self->{ip}, $self->{w});
  170.   $self->{type} = 'ddst';
  171.   $self->{coeff} = $data;
  172.   return $data;
  173. }
  174.  
  175. # Inverse Discrete Sine Transform
  176. sub invddst {
  177.   my $self = shift;
  178.   my $data;
  179.   my $n = $self->{n};
  180.   if (my $arg = shift) {
  181.     die 'Must pass an array reference to invddst'
  182.       unless ref($arg) eq 'ARRAY';
  183.     die "Size of data set must be $n"
  184.       unless $n == @$arg;
  185.     $data = [ @$arg ];
  186.   }
  187.   else {
  188.     die 'Must invert data created with ddst'
  189.       unless $self->{type} eq 'ddst';
  190.     $data =  [ @{$self->{coeff}} ];
  191.   }
  192.   $data->[0] *= 0.5;
  193.   _ddst($n, 1, $data, $self->{ip}, $self->{w});
  194.   $_ *= 2.0/$n for (@$data);
  195.   return $data;
  196. }
  197.  
  198. # Cosine Transform of RDFT (Real Symmetric DFT)
  199. sub dfct {
  200.   my $self = shift;
  201.   my $np1 = $self->{n};
  202.   my $n = $np1 - 1;
  203.   die 'size of data must be an integer power of 2' unless check_n($n);
  204.   my $nt = int(2 + $n/2);
  205.   my $t = [];
  206.   my $data = [ @{$self->{data}} ];
  207.   pdfct($nt, $n, $data, $t, $self->{ip}, $self->{w});
  208.   $self->{type} = 'dfct';
  209.   $self->{coeff} = $data;
  210.   return $data;
  211. }
  212.  
  213. # Inverse Cosine Transform of RDFT (Real Symmetric DFT)
  214. sub invdfct {
  215.   my $self = shift;
  216.   my $data;
  217.   my $np1 = $self->{n};
  218.   my $n = $np1 - 1;
  219.   if (my $arg = shift) {
  220.     die 'Must pass an array reference to invdfct'
  221.       unless ref($arg) eq 'ARRAY';
  222.     die "Size of data set must be $n"
  223.       unless $np1 == @$data;
  224.     $data = [ @$arg ];
  225.   }
  226.   else {
  227.     die 'Must invert data created with dfct'
  228.       unless $self->{type} eq 'dfct';
  229.     $data =  [ @{$self->{coeff}} ];
  230.   }
  231.   my $nt = int(2 + $n/2);
  232.   my $t = [];
  233.   $data->[0] *= 0.5;
  234.   $data->[$n] *= 0.5;
  235.   pdfct($nt, $n, $data, $t, $self->{ip}, $self->{w});
  236.   $data->[0] *= 0.5;
  237.   $data->[$n] *= 0.5;
  238.   $_ *= 2.0/$n for (@$data);
  239.   return $data;
  240. }
  241.  
  242. # Sine Transform of RDFT (Real Anti-symmetric DFT)
  243. sub dfst {
  244.   my $self = shift;
  245.   my $n = $self->{n};
  246.   die 'size of data must be an integer power of 2' unless check_n($n);
  247.   my $data = [ @{$self->{data}} ];
  248.   my $nt = int(2 + $n/2);
  249.   my $t = [];
  250.   pdfst($nt, $n, $data, $t, $self->{ip}, $self->{w});
  251.   $self->{type} = 'dfst';
  252.   $self->{coeff} = $data;
  253.   return $data;
  254. }
  255.  
  256. # Inverse Sine Transform of RDFT (Real Anti-symmetric DFT)
  257. sub invdfst {
  258.   my $self = shift;
  259.   my $n = $self->{n};
  260.   my $data;
  261.   if (my $arg = shift) {
  262.     die 'Must pass an array reference to invdfst'
  263.       unless ref($arg) eq 'ARRAY';
  264.     die "Size of data set must be $n"
  265.       unless $n == @$arg;
  266.     $data = [ @$arg ];
  267.   }
  268.   else {
  269.     die 'Must invert data created with dfst'
  270.       unless $self->{type} eq 'dfst';
  271.     $data =  [ @{$self->{coeff}} ];
  272.   }
  273.   my $nt = int(2 + $n/2);
  274.   my $t = [];
  275.   pdfst($nt, $n, $data, $t, $self->{ip}, $self->{w});
  276.   $_ *= 2.0/$n for (@$data);
  277.   return $data;
  278. }
  279.  
  280. # check if $n is a power of 2
  281. sub check_n {
  282.   my $n = shift;
  283.   my $y = log($n) / 0.693147180559945309417;
  284.   return abs($y-int($y)) < 1e-6 ? 1 : 0;
  285. }
  286.  
  287. sub correl {
  288.   my ($self, $other) = @_;
  289.   my $n = $self->{n};
  290.   my $d1 = $self->{type} ? 
  291.     ($self->{type} eq 'rdft' ? [ @{$self->{coeff}} ] :
  292.     die 'correl must involve a real function' ) :
  293.       $self->rdft &&  [ @{$self->{coeff}} ];
  294.   my $d2 = [];
  295.   if (ref($other) eq 'Math::FFT') {
  296.     $d2 = $other->{type} ? 
  297.       ($other->{type} eq 'rdft' ? [ @{$other->{coeff}}] :
  298.        die 'correl must involve a real function' ) :
  299.      $other->rdft && [ @{$other->{coeff}}];
  300.   }
  301.   elsif (ref($other) eq 'ARRAY') {
  302.     $d2 = [ @$other ];
  303.     _rdft($n, 1, $d2, $self->{ip}, $self->{w});
  304.   }
  305.   else {
  306.     die 'Must call correl with either a Math::FFT object or an array ref';
  307.   }
  308.   my $corr = [];
  309.   _correl($n, $corr, $d1, $d2, $self->{ip}, $self->{w});
  310.   return $corr;
  311. }
  312.  
  313. sub convlv {
  314.   my ($self, $r) = @_;
  315.   die 'Must call convlv with an array reference for the response data'
  316.     unless ref($r) eq 'ARRAY';
  317.   my $respn = [ @$r ];
  318.   my $m = @$respn;
  319.   die 'size of response data must be an odd integer' unless $m % 2 == 1;
  320.   my $n = $self->{n};
  321.   my $d1 = $self->{type} ? 
  322.     ($self->{type} eq 'rdft' ? [ @{$self->{coeff}} ] :
  323.     die 'correl must involve a real function' ) :
  324.       $self->rdft &&  [ @{$self->{coeff}} ];
  325.   for (my $i=1; $i<=($m-1)/2; $i++) {
  326.     $respn->[$n-$i] = $respn->[$m-$i];
  327.   }
  328.   for (my $i=($m+3)/2; $i<=$n-($m-1)/2; $i++) {
  329.     $respn->[$i-1] = 0.0;
  330.   }
  331.   my $convlv = [];
  332.   _convlv($n, $convlv, $d1, $respn, $self->{ip}, $self->{w});
  333.   return $convlv;
  334. }
  335.  
  336. sub deconvlv {
  337.   my ($self, $r) = @_;
  338.   die 'Must call deconvlv with an array reference for the response data'
  339.     unless ref($r) eq 'ARRAY';
  340.   my $respn = [ @$r ];
  341.   my $m = @$respn;
  342.   die 'size of response data must be an odd integer' unless $m % 2 == 1;
  343.   my $n = $self->{n};
  344.   my $d1 = $self->{type} ? 
  345.     ($self->{type} eq 'rdft' ? [ @{$self->{coeff}} ] :
  346.     die 'correl must involve a real function' ) :
  347.       $self->rdft &&  [ @{$self->{coeff}} ];
  348.   for (my $i=1; $i<=($m-1)/2; $i++) {
  349.     $respn->[$n-$i] = $respn->[$m-$i];
  350.   }
  351.   for (my $i=($m+3)/2; $i<=$n-($m-1)/2; $i++) {
  352.     $respn->[$i-1] = 0.0;
  353.   }
  354.   my $convlv = [];
  355.   if (_deconvlv($n, $convlv, $d1, $respn, $self->{ip}, $self->{w}) != 0) {
  356.     die "Singularity encountered for response in deconvlv";
  357.   }
  358.   return $convlv;
  359. }
  360.  
  361. sub spctrm {
  362.   my ($self, %args) = @_;
  363.   my %accept = map {$_ => 1} qw(window segments number overlap);
  364.   for (keys %args) {
  365.     die "`$_' is not a valid argument to spctrm" if not $accept{$_};
  366.   }
  367.   my $win_fun = $args{window};
  368.   if ($win_fun and ref($win_fun) ne 'CODE') {
  369.     my %accept = map {$_ => 1} qw(hamm welch bartlett);
  370.     die "`$win_fun' is not a known window function in spctrm"
  371.       if not $accept{$win_fun};
  372.   }
  373.   die 'Please specify a value for "segments" in spctrm()'
  374.     if ($args{number} and ! $args{segments}); 
  375.   my $n = $self->{n};
  376.   my $d;
  377.   my $n2 = 0;
  378.   my $spctrm = [];
  379.   my $win_sub = {
  380.          'hamm' => sub {
  381.            my ($j, $n) = @_;
  382.            my $pi = 4.0*atan2(1,1);
  383.            return (1 - cos(2*$pi*$j/$n))/2;
  384.          },
  385.          'welch' => sub {
  386.            my ($j, $n) = @_;
  387.            return 1 - 4*($j-$n/2)*($j-$n/2)/$n/$n;
  388.          },
  389.          'bartlett' => sub {
  390.            my ($j, $n) = @_;
  391.            return 1 - abs(2*($j-$n/2)/$n);
  392.          },
  393.         };
  394.   if (not $args{segments} or $args{segments} == 1) {
  395.     if ($win_fun) {
  396.       $d = [ @{$self->{data}}];
  397.       $win_fun = $win_sub->{$win_fun} if ref($win_fun) ne 'CODE';
  398.       for (my $j=0; $j<$n; $j++) {
  399.     my $w = $win_fun->($j, $n);
  400.     $d->[$j] *= $w;
  401.     $n2 += $w * $w;
  402.       }
  403.       $n2 *= $n;
  404.       _spctrm($n, $spctrm, $d, $self->{ip}, $self->{w}, $n2, 1);
  405.     }
  406.     else {
  407.       $d = $self->{type} ? 
  408.     ($self->{type} eq 'rdft' ? $self->{coeff} :
  409.      die 'correl must involve a real function' ) :
  410.        $self->rdft && $self->{coeff};
  411.       $n2 = $n*$n;
  412.       _spctrm($n, $spctrm, $d, $self->{ip}, $self->{w}, $n2, 0);
  413.     }
  414.   }
  415.   else {
  416.     $d = [ @{$self->{data}}];
  417.     my ($data, @w);
  418.     my $k = $args{segments};
  419.     my $m = $args{number};
  420.     die 'Please specify a value for "number" in spctrm()'
  421.       if ($k and ! $m); 
  422.     die "number ($m) must an integer power of 2" unless check_n($m);
  423.     my $m2 = $m+$m;
  424.     my $overlap = $args{overlap};
  425.     my $N = $overlap ? ($k+1)*$m : 2*$k*$m;
  426.     die "Need $N data points (data only has $n)" if $N > $n;
  427.     if ($win_fun) {
  428.       $win_fun = $win_sub->{$win_fun} if ref($win_fun) ne 'CODE';
  429.       for (my $j=0; $j<$m2; $j++) {
  430.     $w[$j] = $win_fun->($j, $m2);
  431.     $n2 += $w[$j]*$w[$j];
  432.       }
  433.     }
  434.     else {
  435.       $n2 = $m2;      
  436.     }
  437.     if ($overlap) {
  438.       my @old =  splice(@$d, 0, $m);
  439.       for (0..$k-1) {
  440.     push @{$data->[$_]}, @old;
  441.     my @new = splice(@$d, 0, $m);
  442.     push @{$data->[$_]}, @new;
  443.     @old = @new;
  444.     if ($win_fun) {
  445.       my $j=0;
  446.       $data->[$_] = [ map {$w[$j++]*$_} @{$data->[$_]}];
  447.     }
  448.       }
  449.     }
  450.     else {
  451.       for (0..$k-1) {
  452.     push @{$data->[$_]}, splice(@$d, 0, $m2);
  453.     if ($win_fun) {
  454.       my $j=0;
  455.       $data->[$_] = [ map {$w[$j++]*$_} @{$data->[$_]}];
  456.     }
  457.       }
  458.     }
  459.     my $tmp = [];
  460.     my $nip = int(3 + sqrt($m2));
  461.     my $nw = int(2 + 5*$m2/4);
  462.     my $ip = pack("i$nip", ());
  463.     my $w = pack("d$nw", ());
  464.     _spctrm_bin($k, $m2, $spctrm, $data, \$ip, \$w, $n2, $tmp);
  465.   }
  466.   return $spctrm;
  467. }
  468.  
  469. sub mean {
  470.   my $self = shift;
  471.   my $sum = 0;
  472.   my ($n, $data);
  473.   my $flag = 0;
  474.   if ($data = shift) {
  475.     die 'Must call with an array reference'
  476.       unless ref($data) eq 'ARRAY';
  477.     $n = @$data;
  478.     $flag = 1;
  479.   }
  480.   else {
  481.     $data = $self->{data};
  482.     $n = $self->{n};
  483.   }
  484.   $sum += $_ for @$data;
  485.   my $mean = $sum / $n;
  486.   $self->{mean} = $mean unless $flag == 1;
  487.   return $mean;
  488. }
  489.  
  490. sub stdev {
  491.   my $self = shift;
  492.   my ($n, $data, $mean);
  493.   if ($data = shift) {
  494.     die 'Must call with an array reference'
  495.       unless ref($data) eq 'ARRAY';
  496.     $n = @$data;
  497.     $mean = $self->mean($data);
  498.   }
  499.   else {
  500.     $data = $self->{data};
  501.     $n = $self->{n};
  502.     $mean = $self->{mean} || $self->mean;
  503.   }
  504.   die 'Cannot find the standard deviation with n = 1'
  505.     if $n == 1;
  506.   my $sum = 0;
  507.   $sum += ($_ - $mean)*($_ - $mean) for @$data;
  508.   return sqrt($sum / ($n-1));
  509. }
  510.  
  511. sub range {
  512.   my $self = shift;
  513.   my ($n, $data);
  514.   if ($data = shift) {
  515.     die 'Must call with an array reference'
  516.       unless ref($data) eq 'ARRAY';
  517.     $n = @$data;
  518.   }
  519.   else {
  520.     $data = $self->{data};
  521.     $n = $self->{n};
  522.   }
  523.   my $min = $data->[0];
  524.   my $max = $data->[0];
  525.   for (@$data) {
  526.     $min = $_ if $_ < $min;
  527.     $max = $_ if $_ > $max;
  528.   }
  529.   return ($min, $max);
  530. }
  531.  
  532. sub median {
  533.   my $self = shift;
  534.   my ($n, $data);
  535.   if ($data = shift) {
  536.     die 'Must call with an array reference'
  537.       unless ref($data) eq 'ARRAY';
  538.     $n = @$data;
  539.   }
  540.   else {
  541.     $data = $self->{data};
  542.     $n = $self->{n};
  543.   }
  544.   my @sorted = sort {$a <=> $b} @$data;
  545.   return $n % 2 == 1 ?
  546.     $sorted[($n-1)/2] : ($sorted[$n/2] + $sorted[$n/2-1])/2;
  547. }
  548.  
  549. # Autoload methods go after =cut, and are processed by the autosplit program.
  550.  
  551.  
  552. 1;
  553.  
  554. __END__
  555.  
  556. =head1 NAME
  557.  
  558. Math::FFT - Perl module to calculate Fast Fourier Transforms
  559.  
  560. =head1 SYNOPSIS
  561.  
  562.   use Math::FFT;
  563.   my $PI = 3.1415926539;
  564.   my $N = 64;
  565.   my ($series, $other_series);
  566.   for (my $k=0; $k<$N; $k++) {
  567.       $series->[$k] = sin(4*$k*$PI/$N) + cos(6*$k*$PI/$N);
  568.   }
  569.   my $fft = new Math::FFT($series);
  570.   my $coeff = $fft->rdft();
  571.   my $spectrum = $fft->spctrm;
  572.   my $original_data = $fft->invrdft($coeff);
  573.  
  574.   for (my $k=0; $k<$N; $k++) {
  575.       $other_series->[$k] = sin(16*$k*$PI/$N) + cos(8*$k*$PI/$N);
  576.   }
  577.   my $other_fft = $fft->clone($other_series);
  578.   my $other_coeff = $other_fft->rdft();
  579.   my $correlation = $fft->correl($other_fft);
  580.  
  581. =head1 DESCRIPTION
  582.  
  583. This module implements some algorithms for calculating
  584. Fast Fourier Transforms for one-dimensional data sets of size 2^n.
  585. The data, assumed to arise from a constant sampling rate, is
  586. represented by an array reference C<$data> (as described in the
  587. methods below), which is then used to create a C<Math::FFT> object as
  588.  
  589.   my $fft = new Math::FFT($data);
  590.  
  591. The methods available include the following.
  592.  
  593. =head2 FFT METHODS
  594.  
  595. =over
  596.  
  597. =item C<$coeff = $fft-E<gt>cdft();>
  598.  
  599. This calculates the complex discrete Fourier transform
  600. for a data set C<x[j]>. Here, C<$data> is a reference to an 
  601. array C<data[0...2*n-1]> holding the data
  602.  
  603.   data[2*j] = Re(x[j]),
  604.   data[2*j+1] = Im(x[j]), 0<=j<n
  605.  
  606. An array reference C<$coeff> is returned consisting of
  607.  
  608.   coeff[2*k] = Re(X[k]),
  609.   coeff[2*k+1] = Im(X[k]), 0<=k<n
  610.  
  611. where
  612.  
  613.    X[k] = sum_j=0^n-1 x[j]*exp(2*pi*i*j*k/n), 0<=k<n
  614.  
  615. =item C<$orig_data = $fft-E<gt>invcdft([$coeff]);>
  616.  
  617. Calculates the inverse complex discrete Fourier transform
  618. on a data set C<x[j]>. If C<$coeff> is not given, it will be set 
  619. equal to an earlier call to C<$fft-E<gt>cdft()>. C<$coeff> is 
  620. a reference to an array C<coeff[0...2*n-1]> holding the data
  621.  
  622.   coeff[2*j] = Re(x[j]),
  623.   coeff[2*j+1] = Im(x[j]), 0<=j<n
  624.  
  625. An array reference C<$orig_data> is returned consisting of
  626.  
  627.   orig_data[2*k] = Re(X[k]),
  628.   orig_data[2*k+1] = Im(X[k]), 0<=k<n
  629.  
  630. where, excluding the scale,
  631.  
  632.    X[k] = sum_j=0^n-1 x[j]*exp(-2*pi*i*j*k/n), 0<=k<n
  633.  
  634. A scaling C<$orig_data-E<gt>[$i] *= 2.0/$n> is then done so that
  635. C<$orig_data> coincides with the original C<$data>.
  636.  
  637. =item C<$coeff = $fft-E<gt>rdft();>
  638.  
  639. This calculates the real discrete Fourier transform
  640. for a data set C<x[j]>. On input, $data is a reference to an
  641. array C<data[0...n-1]> holding the data. An array reference
  642. C<$coeff> is returned consisting of
  643.  
  644.   coeff[2*k] = R[k], 0<=k<n/2
  645.   coeff[2*k+1] = I[k], 0<k<n/2
  646.   coeff[1] = R[n/2]
  647.  
  648. where
  649.  
  650.   R[k] = sum_j=0^n-1 data[j]*cos(2*pi*j*k/n), 0<=k<=n/2
  651.   I[k] = sum_j=0^n-1 data[j]*sin(2*pi*j*k/n), 0<k<n/2
  652.  
  653. =item C<$orig_data = $fft-E<gt>invrdft([$coeff]);>
  654.  
  655. Calculates the inverse real discrete Fourier transform
  656. on a data set C<coeff[j]>. If C<$coeff> is not given, it will be set 
  657. equal to an earlier call to C<$fft-E<gt>rdft()>. C<$coeff> 
  658. is a reference to an array C<coeff[0...n-1]> holding the data
  659.  
  660.   coeff[2*j] = R[j], 0<=j<n/2
  661.   coeff[2*j+1] = I[j], 0<j<n/2
  662.   coeff[1] = R[n/2]
  663.  
  664. An array reference C<$orig_data> is returned where, excluding the scale,
  665.  
  666.   orig_data[k] = (R[0] + R[n/2]*cos(pi*k))/2 + 
  667.     sum_j=1^n/2-1 R[j]*cos(2*pi*j*k/n) + 
  668.       sum_j=1^n/2-1 I[j]*sin(2*pi*j*k/n), 0<=k<n
  669.  
  670. A scaling C<$orig_data-E<gt>[$i] *= 2.0/$n> is then done so that
  671. C<$orig_data> coincides with the original C<$data>.
  672.  
  673. =item C<$coeff = $fft-E<gt>ddct();>
  674.  
  675. Computes the discrete cosine tranform on a data set
  676. C<data[0...n-1]> contained in an array reference C<$data>. An
  677. array reference C<$coeff> is returned consisting of
  678.  
  679.   coeff[k] = C[k], 0<=k<n
  680.  
  681. where
  682.  
  683.   C[k] = sum_j=0^n-1 data[j]*cos(pi*(j+1/2)*k/n), 0<=k<n 
  684.  
  685. =item C<$orig_data = $fft-E<gt>invddct([$coeff]);>
  686.  
  687. Computes the inverse discrete cosine tranform on a data set
  688. C<coeff[0...n-1]> contained in an array reference C<$coeff>. 
  689. If C<$coeff> is not given, it will be set equal to an earlier 
  690. call to C<$fft-E<gt>ddct()>. An array reference C<$orig_data> 
  691. is returned consisting of
  692.  
  693.   orig_data[k] = C[k], 0<=k<n
  694.  
  695. where, excluding the scale,
  696.  
  697.   C[k] = sum_j=0^n-1 coeff[j]*cos(pi*j*(k+1/2)/n), 0<=k<n
  698.  
  699. A scaling C<$orig_data-E<gt>[$i] *= 2.0/$n> is then done so that
  700. C<$orig_data> coincides with the original C<$data>.
  701.  
  702. =item C<$coeff = $fft-E<gt>ddst();>
  703.  
  704. Computes the discrete sine transform of a data set 
  705. C<data[0...n-1]> contained in an array reference C<$data>. An
  706. array reference C<$coeff> is returned consisting of
  707.  
  708.  coeff[k] = S[k], 0<k<n
  709.  coeff[0] = S[n]
  710.  
  711. where
  712.  
  713.  S[k] = sum_j=0^n-1 data[j]*sin(pi*(j+1/2)*k/n), 0<k<=n
  714.  
  715. =item C<$orig_data = $fft-E<gt>invddst($coeff);>
  716.  
  717. Computes the inverse discrete sine transform of a data set 
  718. C<coeff[0...n-1]> contained in an array reference C<$coeff>, arranged as 
  719.  
  720.  coeff[j] = A[j], 0<j<n
  721.  coeff[0] = A[n]
  722.  
  723. If C<$coeff> is not given, it will be set equal to an earlier 
  724. call to C<$fft-E<gt>ddst()>. An array reference C<$orig_data> 
  725. is returned consisting of
  726.  
  727.  orig_data[k] = S[k], 0<=k<n
  728.  
  729. where, excluding a scale,
  730.  
  731.  S[k] =  sum_j=1^n A[j]*sin(pi*j*(k+1/2)/n), 0<=k<n
  732.  
  733. The scaling C<$a-E<gt>[$i] *= 2.0/$n> is then done so that
  734. C<$orig_data> coincides with the original C<$data>.
  735.  
  736. =item C<$coeff = $fft-E<gt>dfct();>
  737.  
  738. Computes the real symmetric discrete Fourier transform of a
  739. data set C<data[0...n]> contained in the array reference C<$data>. An
  740. array reference C<$coeff> is returned consisting of 
  741.  
  742.   coeff[k] = C[k], 0<=k<=n
  743.  
  744. where
  745.  
  746.   C[k] = sum_j=0^n data[j]*cos(pi*j*k/n), 0<=k<=n
  747.  
  748. =item C<$orig_data = $fft-E<gt>invdfct($coeff);>
  749.  
  750. Computes the inverse real symmetric discrete Fourier transform of a
  751. data set C<coeff[0...n]> contained in the array reference C<$coeff>. 
  752. If C<$coeff> is not given, it will be set equal to an earlier 
  753. call to C<$fft-E<gt>dfct()>. An array reference C<$orig_data> 
  754. is returned consisting of
  755.  
  756.   orig_data[k] = C[k], 0<=k<=n
  757.  
  758. where, excluding the scale,
  759.  
  760.   C[k] = sum_j=0^n coeff[j]*cos(pi*j*k/n), 0<=k<=n
  761.  
  762. A scaling C<$coeff-E<gt>[0] *= 0.5>, C<$coeff-E<gt>[$n] *= 0.5>, and 
  763. C<$orig_data-E<gt>[$i] *= 2.0/$n> is then done so that
  764. C<$orig_data> coincides with the original C<$data>.
  765.  
  766. =item C<$coeff = $fft-E<gt>dfst();>
  767.  
  768. Computes the real anti-symmetric discrete Fourier transform of a
  769. data set C<data[0...n-1]> contained in the array reference C<$data>. An
  770. array reference C<$coeff> is returned consisting of 
  771.  
  772.   coeff[k] = C[k], 0<k<n
  773.  
  774. where
  775.  
  776.   C[k] = sum_j=0^n data[j]*sin(pi*j*k/n), 0<k<n
  777.  
  778. (C<coeff[0]> is used for a work area)
  779.  
  780. =item C<$orig_data = $fft-E<gt>invdfst($coeff);>
  781.  
  782. Computes the inverse real anti-symmetric discrete Fourier transform of a
  783. data set C<coeff[0...n-1]> contained in the array reference C<$coeff>.
  784. If C<$coeff> is not given, it will be set equal to an earlier 
  785. call to C<$fft-E<gt>dfst()>. An array reference C<$orig_data> is 
  786. returned consisting of
  787.  
  788.   orig_data[k] = C[k], 0<k<n
  789.  
  790. where, excluding the scale,
  791.  
  792.   C[k] = sum_j=0^n coeff[j]*sin(pi*j*k/n), 0<k<n
  793.  
  794. A scaling C<$orig_data-E<gt>[$i] *= 2.0/$n> is then done so that
  795. C<$orig_data> coincides with the original C<$data>.
  796.  
  797. =back
  798.  
  799. =head2 CLONING
  800.  
  801. The algorithm used in the transforms makes use of arrays for a work 
  802. area and for a cos/sin lookup table dependent only on the size of 
  803. the data set. These arrays are initialized when the C<Math::FFT> object 
  804. is created and then are populated when a transform method is first 
  805. invoked. After this, they persist for the lifetime of the object.
  806.  
  807. This aspect is exploited in a C<cloning> method; if a C<Math::FFT>
  808. object is created for a data set C<$data1> of size C<N>:
  809.  
  810.   $fft1 = new Math::FFT($data1);
  811.  
  812. then a new C<Math::FFT> object can be created for a second data 
  813. set C<$data2> of the I<same> size C<N> by
  814.  
  815.    $fft2 = $fft1->clone($data2);
  816.  
  817. The C<$fft2> object will copy the reuseable work area and
  818. lookup table calculated from C<$fft1>.
  819.  
  820. =head2 APPLICATIONS
  821.  
  822. This module includes some common applications - correlation,
  823. convolution and deconvolution, and power spectrum - that
  824. arise with real data sets. The conventions used here
  825. follow that of I<Numerical Recipes in C>, by Press, Teukolsky,
  826. Vetterling, and Flannery, in which further details of the
  827. algorithms are given. Note in particular the treatment of end
  828. effects by zero padding, which is assumed to be done by the
  829. user, if required.
  830.  
  831. =over
  832.  
  833. =item Correlation
  834.  
  835. The correlation between two functions is defined as
  836.  
  837.              /
  838.    Corr(t) = | ds g(s+t) h(s) 
  839.              /
  840.  
  841. This may be calculated, for two array references C<$data1>
  842. and C<$data2> of the same size C<$n>, as either
  843.  
  844.    $fft1 = new Math::FFT($data1);
  845.    $fft2 = new Math::FFT($data2);
  846.    $corr = $fft1->correl($fft2);
  847.  
  848. or as
  849.  
  850.    $fft1 = new Math::FFT($data1);
  851.    $corr = $fft1->correl($data2);
  852.  
  853. The array reference C<$corr> is returned in wrap-around 
  854. order - correlations at increasingly positive lags are in 
  855. C<$corr-E<gt>[0]> (zero lag) on up to C<$corr-E<gt>[$n/2-1]>, 
  856. while correlations at increasingly negative lags are in 
  857. C<$corr-E<gt>[$n-1]> on down to C<$corr-E<gt>[$n/2]>. The sign 
  858. convention used is such that if C<$data1> lags C<$data2> (that 
  859. is, is shifted to the right), then C<$corr> will show a peak 
  860. at positive lags.
  861.  
  862. =item Convolution
  863.  
  864. The convolution of two functions is defined as
  865.  
  866.                /
  867.    Convlv(t) = | ds g(s) h(t-s) 
  868.                /
  869.  
  870. This is similar to calculating the correlation between the
  871. two functions, but typically the functions here have a quite
  872. different physical interpretation - one is a signal which 
  873. persists indefinitely in time, and the other is a response 
  874. function of limited duration. The convolution may be calculated, 
  875. for two array references C<$data> and C<$respn>, as
  876.  
  877.    $fft = new Math::FFT($data);
  878.    $convlv = $fft->convlv($respn);
  879.  
  880. with the returned C<$convlv> being an array reference. The method 
  881. assumes that the response function C<$respn> has an I<odd> number 
  882. of elements C<$m> less than or equal to the number of elements C<$n> 
  883. of C<$data>. C<$respn> is assumed to be stored in wrap-around order - 
  884. the first half contains the response at positive times, while the 
  885. second half, counting down from C<$respn-E<gt>[$m-1]>, contains the
  886. response at negative times.
  887.  
  888. =item Deconvolution
  889.  
  890. Deconvolution undoes the effects of convoluting a signal
  891. with a known response function. In other words, in the relation
  892.  
  893.                /
  894.    Convlv(t) = | ds g(s) h(t-s) 
  895.                /
  896.  
  897. deconvolution reconstructs the original signal, given the convolution
  898. and the response function. The method is implemented, for two array 
  899. references C<$data> and C<$respn>, as
  900.  
  901.    $fft = new Math::FFT($data);
  902.    $deconvlv = $fft->deconvlv($respn);
  903.  
  904. As a result, if the convolution of a data set C<$data> with
  905. a response function C<$respn> is calculated as
  906.  
  907.    $fft1 = new Math::FFT($data);
  908.    $convlv = $fft1->convlv($respn);
  909.  
  910. then the deconvolution
  911.  
  912.    $fft2 = new Math::FFT($convlv);
  913.    $deconvlv = $fft2->deconvlv($respn);
  914.  
  915. will give an array reference C<$deconvlv> containing the
  916. same elements as the original data C<$data>.
  917.  
  918. =item Power Spectrum
  919.  
  920. If the FFT of a real function of C<N> elements is calculated, 
  921. the C<N/2+1> elements of the power spectrum are defined, in terms 
  922. of the (complex) Fourier coefficients C<C[k]>, as
  923.  
  924.    P[0] = |C[0]|^2 / N^2
  925.    P[k] = 2 |C[k]|^2 / N^2   (k = 1, 2 ,..., N/2-1)
  926.    P[N/2] = |C[N/2]|^2 / N^2
  927.  
  928. Often for these purposes the data is partitioned into C<K>
  929. segments, each containing C<2M> elements. The power spectrum
  930. for each segment is calculated, and the net power spectrum
  931. is the average of all of these segmented spectra.
  932.  
  933. Partitioning may be done in one of two ways: I<non-overlapping> and
  934. I<overlapping>. Non-overlapping is useful when the data set
  935. is gathered in real time, where the number of data points
  936. can be varied at will. Overlapping is useful where there is
  937. a fixed number of data points. In non-overlapping, the first
  938. <2M> elements constitute segment 1, the next C<2M> elements
  939. are segment 2, and so on up to segment C<K>, for a total of
  940. C<2KM> sampled points. In overlapping, the first and second
  941. C<M> elements are segment 1, the second and third C<M> elements
  942. are segment 2, and so on, for a total of C<(K+1)M> sampled points.
  943.  
  944. A problem that may arise in this procedure is I<leakage>: the
  945. power spectrum calculated for one bin contains contributions
  946. from nearby bins. To lessen this effect I<data windowing> is
  947. often used: multiply the original data C<d[j]> by a window
  948. function C<w[j]>, where j = 0, 1, ..., N-1. Some popular choices 
  949. of such functions are
  950.  
  951.               | j - N/2 |
  952.   w[j] = 1 -  | ------- |     ... Bartlett   
  953.               |   N/2   |
  954.  
  955.  
  956.               / j - N/2 \ 2
  957.   w[j] = 1 -  | ------- |     ... Welch  
  958.               \   N/2   /
  959.  
  960.  
  961.            1   /                    \
  962.   w[j] =  ---  |1 - cos(2 pi j / N) |     ... Hamm  
  963.            2   \                    /
  964.  
  965.  
  966. The C<spctrm> method, used as
  967.  
  968.     $fft = Math::FFT->new($data);
  969.     $spectrum = $fft->spctrm([$key => $value, ...]);
  970.  
  971. returns an array reference C<$spectrum> representing the power 
  972. spectrum for a data set represented by an array reference C<$data>.
  973. The options available are
  974.  
  975. =over
  976.  
  977. =item C<window =E<gt> window_name>
  978.  
  979. This specifies the window function; if not given, no such
  980. function is used. Accepted values (see above) are C<"bartlett">, 
  981. C<"welch">, C<"hamm">, and C<\&my_window>, where C<my_window> is a 
  982. user specified subroutine which must be of the form, for example,
  983.  
  984.    sub my_window {
  985.       my ($j, $n) = @_;
  986.       return 1 - abs(2*($j-$n/2)/$n);
  987.    }
  988.  
  989. which implements the Bartlett window.
  990.  
  991. =item C<overlap =E<gt> 1>
  992.  
  993. This specifies whether overlapping should be done; if true (1),
  994. overlapping will be used, whereas if false (0), or not
  995. specified, no overlapping is used.
  996.  
  997. =item C<segments =E<gt> n>
  998.  
  999. This specifies that the data will be partitioned into C<n>
  1000. segments. If not specified, no segmentation will be done.
  1001.  
  1002. =item C<number =E<gt> m>
  1003.  
  1004. This specifies that C<2m> data points will be used for 
  1005. each segment, and must be a power of 2. The power 
  1006. spectrum returned will consist of C<m+1> elements.
  1007.  
  1008. =back
  1009.  
  1010. =back
  1011.  
  1012. =head2 STATISTICAL FUNCTIONS
  1013.  
  1014. For convenience, a number of common statistical functions are 
  1015. included for analyzing real data. After creating the object as
  1016.  
  1017.   my $fft = new Math::FFT($data);
  1018.  
  1019. for a data set represented by the array reference C<$data>
  1020. of size C<N>, these methods may be called as follows.
  1021.  
  1022. =over
  1023.  
  1024. =item C<$mean = $fft-E<gt>mean([$data]);>
  1025.  
  1026. This returns the mean
  1027.  
  1028.   1/N * sum_j=0^N-1 data[j]
  1029.  
  1030. If an array reference C<$data> is not given, the data set used 
  1031. in creating C<$fft> will be used.
  1032.  
  1033. =item C<$stdev = $fft-E<gt>stdev([$data]);>
  1034.  
  1035. This returns the standard deviation
  1036.  
  1037.   sqrt{ 1/(N-1) * sum_j=0^N-1 (data[j] - mean)**2 }
  1038.  
  1039. If an array reference C<$data> is not given, the data set used 
  1040. in creating C<$fft> will be used.
  1041.  
  1042. =item C<($min, $max) = $fft-E<gt>range([$data]);>
  1043.  
  1044. This returns the minimum and maximum values of the data set.
  1045. If an array reference C<$data> is not given, the data set used 
  1046. in creating C<$fft> will be used.
  1047.  
  1048. =item C<$median = $fft-E<gt>median([$data]);>
  1049.  
  1050. This returns the median of a data set. The median is defined,
  1051. for the I<sorted> data set, as either the middle element, if the
  1052. number of elements is odd, or as the interpolated value of
  1053. the the two values on either side of the middle, if the number
  1054. of elements is even. If an array reference C<$data> is not given, 
  1055. the data set used in creating C<$fft> will be used.
  1056.  
  1057. =back
  1058.  
  1059. =head1 BUGS
  1060.  
  1061. Please report any to Randy Kobes <randy@theoryx5.uwinnipeg.ca>
  1062.  
  1063. =head1 SEE ALSO
  1064.  
  1065. L<Math::Pari> and L<PDL>
  1066.  
  1067. =head1 COPYRIGHT
  1068.  
  1069. The algorithm used in this module to calculate the Fourier
  1070. transforms is based on the C routine of fft4g.c available
  1071. at http://momonga.t.u-tokyo.ac.jp/~ooura/fft.html, which is
  1072. copyrighted 1996-99 by Takuya OOURA. The file arrays.c included 
  1073. here to handle passing arrays to and from C comes from the PGPLOT 
  1074. module of Karl Glazebrook <kgb@aaoepp.aao.gov.au>. The perl code 
  1075. of Math::FFT is copyright 2000 by Randy Kobes, and is distributed 
  1076. under the same terms as Perl itself.
  1077.  
  1078. =cut
  1079.