FFT.i

   1 %module "Math::GSL::FFT"
   2 %include "typemaps.i"
   3 %include "gsl_typemaps.i"
   4
   5
   6 %typemap(argout) (double data[], const size_t stride, const size_t n) {
   7     int i=0;
   8     AV* tempav = newAV();
   9
  10     while( i < $3 ) {
  11         av_push(tempav, newSVnv((double) $1[i]));
  12         i++;
  13     }
  14
  15     $result = sv_2mortal( newRV_noinc( (SV*) tempav) );
  16     argvi++;
  17 }
  18
  19 %typemap(argout) (gsl_complex_packed_array data[], const size_t stride, const size_t n) {
  20     int i=0;
  21     AV* tempav = newAV();
  22
  23     while( i < 2*$3 ) {
  24         av_push(tempav, newSVnv((double) $1[i]));
  25         i++;
  26     }
  27
  28     $result = sv_2mortal( newRV_noinc( (SV*) tempav) );
  29     argvi++;
  30 }
  31
  32 %{
  33     #include "gsl/gsl_fft.h"
  34     #include "gsl/gsl_fft_complex.h"
  35     #include "gsl/gsl_fft_halfcomplex.h"
  36     #include "gsl/gsl_fft_real.h"
  37     #include "gsl/gsl_complex.h"
  38 %}
  39
  40 %include "gsl/gsl_complex.h"
  41 %include "gsl/gsl_fft.h"
  42 %include "gsl/gsl_fft_complex.h"
  43 %include "gsl/gsl_fft_halfcomplex.h"
  44 %include "gsl/gsl_fft_real.h"
  45
  46 %perlcode %{
  47 @EXPORT_complex = qw/
  48                gsl_fft_complex_radix2_forward
  49                gsl_fft_complex_radix2_backward
  50                gsl_fft_complex_radix2_inverse
  51                gsl_fft_complex_radix2_transform
  52                gsl_fft_complex_radix2_dif_forward
  53                gsl_fft_complex_radix2_dif_backward
  54                gsl_fft_complex_radix2_dif_inverse
  55                gsl_fft_complex_radix2_dif_transform
  56                gsl_fft_complex_wavetable_alloc
  57                gsl_fft_complex_wavetable_free
  58                gsl_fft_complex_workspace_alloc
  59                gsl_fft_complex_workspace_free
  60                gsl_fft_complex_memcpy
  61                gsl_fft_complex_forward
  62                gsl_fft_complex_backward
  63                gsl_fft_complex_inverse
  64                gsl_fft_complex_transform
  65                /;
  66 @EXPORT_halfcomplex = qw/
  67                gsl_fft_halfcomplex_radix2_backward
  68                gsl_fft_halfcomplex_radix2_inverse
  69                gsl_fft_halfcomplex_radix2_transform
  70                gsl_fft_halfcomplex_wavetable_alloc
  71                gsl_fft_halfcomplex_wavetable_free
  72                gsl_fft_halfcomplex_backward
  73                gsl_fft_halfcomplex_inverse
  74                gsl_fft_halfcomplex_transform
  75                gsl_fft_halfcomplex_unpack
  76                gsl_fft_halfcomplex_radix2_unpack
  77                /;
  78 @EXPORT_real = qw/
  79                gsl_fft_real_radix2_transform
  80                gsl_fft_real_wavetable_alloc
  81                gsl_fft_real_wavetable_free
  82                gsl_fft_real_workspace_alloc
  83                gsl_fft_real_workspace_free
  84                gsl_fft_real_transform
  85                gsl_fft_real_unpack
  86              /;
  87 @EXPORT_vars = qw/
  88                 $gsl_fft_forward
  89                 $gsl_fft_backward
  90                 /;
  91 @EXPORT_OK =   (
  92                 @EXPORT_real,
  93                 @EXPORT_complex,
  94                 @EXPORT_halfcomplex,
  95                 @EXPORT_vars,
  96                 );
  97 %EXPORT_TAGS = (
  98                 all         => \@EXPORT_OK,
  99                 real        => \@EXPORT_real,
 100                 complex     => \@EXPORT_complex,
 101                 halfcomplex => \@EXPORT_halfcomplex,
 102                 vars        => \@EXPORT_vars,
 103                );
 104 __END__
 105
 106 =head1 NAME
 107
 108 Math::GSL::FFT - Functions for performing Fast Fourier Transforms (FFT)
 109
 110 =head1 SYNOPSIS
 111
 112 use Math::GSL::FFT qw /:all/;
 113
 114 =head1 DESCRIPTION
 115
 116 Here is a list of all the functions in this module :
 117
 118 =over
 119
 120 =item * C<gsl_fft_complex_radix2_forward($data, $stride, $n) > - This function computes the forward FFTs of length $n with stride $stride, on the array reference $data using an in-place radix-2 decimation-in-time algorithm. The length of the transform $n is restricted to powers of two. For the transform version of the function the sign argument can be either forward (-1) or backward (+1). The functions return a value of $GSL_SUCCESS if no errors were detected, or $GSL_EDOM if the length of the data $n is not a power of two. The complex functions of the FFT module are not yet fully implemented.
 121
 122 =item * C<gsl_fft_complex_radix2_backward >
 123
 124 =item * C<gsl_fft_complex_radix2_inverse >
 125
 126 =item * C<gsl_fft_complex_radix2_transform >
 127
 128 =item * C<gsl_fft_complex_radix2_dif_forward >
 129
 130 =item * C<gsl_fft_complex_radix2_dif_backward >
 131
 132 =item * C<gsl_fft_complex_radix2_dif_inverse >
 133
 134 =item * C<gsl_fft_complex_radix2_dif_transform >
 135
 136 =item * C<gsl_fft_complex_wavetable_alloc($n)> - This function prepares a trigonometric lookup table for a complex FFT of length $n. The function returns a pointer to the newly allocated gsl_fft_complex_wavetable if no errors were detected, and a null pointer in the case of error. The length $n is factorized into a product of subtransforms, and the factors and their trigonometric coefficients are stored in the wavetable. The trigonometric coefficients are computed using direct calls to sin and cos, for accuracy. Recursion relations could be used to compute the lookup table faster, but if an application performs many FFTs of the same length then this computation is a one-off overhead which does not affect the final throughput. The wavetable structure can be used repeatedly for any transform of the same length. The table is not modified by calls to any of the other FFT functions. The same wavetable can be used for both forward and backward (or inverse) transforms of a given length.
 137
 138 =item * C<gsl_fft_complex_wavetable_free($wavetable)> - This function frees the memory associated with the wavetable $wavetable. The wavetable can be freed if no further FFTs of the same length will be needed.
 139
 140 =item * C<gsl_fft_complex_workspace_alloc($n)> - This function allocates a workspace for a complex transform of length $n.
 141
 142 =item * C<gsl_fft_complex_workspace_free($workspace) > - This function frees the memory associated with the workspace $workspace. The workspace can be freed if no further FFTs of the same length will be needed.
 143
 144 =item * C<gsl_fft_complex_memcpy >
 145
 146 =item * C<gsl_fft_complex_forward >
 147
 148 =item * C<gsl_fft_complex_backward >
 149
 150 =item * C<gsl_fft_complex_inverse >
 151
 152 =item * C<gsl_fft_complex_transform >
 153
 154 =item * C<gsl_fft_halfcomplex_radix2_backward($data, $stride, $n)> - This function computes the backwards in-place radix-2 FFT of length $n and stride $stride on the half-complex sequence data stored according the output scheme used by gsl_fft_real_radix2. The result is a real array stored in natural order.
 155
 156 =item * C<gsl_fft_halfcomplex_radix2_inverse($data, $stride, $n)> - This function computes the inverse in-place radix-2 FFT of length $n and stride $stride on the half-complex sequence data stored according the output scheme used by gsl_fft_real_radix2. The result is a real array stored in natural order.
 157
 158 =item * C<gsl_fft_halfcomplex_radix2_transform>
 159
 160 =item * C<gsl_fft_halfcomplex_wavetable_alloc($n)> - This function prepares trigonometric lookup tables for an FFT of size $n real elements. The functions return a pointer to the newly allocated struct if no errors were detected, and a null pointer in the case of error. The length $n is factorized into a product of subtransforms, and the factors and their trigonometric coefficients are stored in the wavetable. The trigonometric coefficients are computed using direct calls to sin and cos, for accuracy. Recursion relations could be used to compute the lookup table faster, but if an application performs many FFTs of the same length then computing the wavetable is a one-off overhead which does not affect the final throughput.
 161     The wavetable structure can be used repeatedly for any transform of the same length. The table is not modified by calls to any of the other FFT functions. The appropriate type of wavetable must be used for forward real or inverse half-complex transforms.
 162
 163 =item * C<gsl_fft_halfcomplex_wavetable_free($wavetable)> - This function frees the memory associated with the wavetable $wavetable. The wavetable can be freed if no further FFTs of the same length will be needed.
 164
 165 =item * C<gsl_fft_halfcomplex_backward >
 166
 167 =item * C<gsl_fft_halfcomplex_inverse >
 168
 169 =item * C<gsl_fft_halfcomplex_transform >
 170
 171 =item * C<gsl_fft_halfcomplex_unpack >
 172
 173 =item * C<gsl_fft_halfcomplex_radix2_unpack >
 174
 175 =item * C<gsl_fft_real_radix2_transform($data, $stride, $n) > - This function computes an in-place radix-2 FFT of length $n and stride $stride on the real array reference $data. The output is a half-complex sequence, which is stored in-place. The arrangement of the half-complex terms uses the following scheme: for k < N/2 the real part of the k-th term is stored in location k, and the corresponding imaginary part is stored in location N-k. Terms with k > N/2 can be reconstructed using the symmetry z_k = z^*_{N-k}. The terms for k=0 and k=N/2 are both purely real, and count as a special case. Their real parts are stored in locations 0 and N/2 respectively, while their imaginary parts which are zero are not stored. The following table shows the correspondence between the output data and the equivalent results obtained by considering the input data as a complex sequence with zero imaginary part,
 176
 177           complex[0].real    =    data[0]
 178           complex[0].imag    =    0
 179           complex[1].real    =    data[1]
 180           complex[1].imag    =    data[N-1]
 181           ...............         ................
 182           complex[k].real    =    data[k]
 183           complex[k].imag    =    data[N-k]
 184           ...............         ................
 185           complex[N/2].real  =    data[N/2]
 186           complex[N/2].imag  =    0
 187           ...............         ................
 188           complex[k'].real   =    data[k]        k' = N - k
 189           complex[k'].imag   =   -data[N-k]
 190           ...............         ................
 191           complex[N-1].real  =    data[1]
 192           complex[N-1].imag  =   -data[N-1]
 193
 194 Note that the output data can be converted into the full complex sequence using the function gsl_fft_halfcomplex_unpack.
 195
 196 =item * C<gsl_fft_real_wavetable_alloc($n)> - This function prepares trigonometric lookup tables for an FFT of size $n real elements. The functions return a pointer to the newly allocated struct if no errors were detected, and a null pointer in the case of error. The length $n is factorized into a product of subtransforms, and the factors and their trigonometric coefficients are stored in the wavetable. The trigonometric coefficients are computed using direct calls to sin and cos, for accuracy. Recursion relations could be used to compute the lookup table faster, but if an application performs many FFTs of the same length then computing the wavetable is a one-off overhead which does not affect the final throughput.
 197     The wavetable structure can be used repeatedly for any transform of the same length. The table is not modified by calls to any of the other FFT functions. The appropriate type of wavetable must be used for forward real or inverse half-complex transforms.
 198
 199 =item * C<gsl_fft_real_wavetable_free($wavetable)> - This function frees the memory associated with the wavetable $wavetable. The wavetable can be freed if no further FFTs of the same length will be needed.
 200
 201 =item * C<gsl_fft_real_workspace_alloc($n)> - This function allocates a workspace for a real transform of length $n. The same workspace can be used for both forward real and inverse halfcomplex transforms.
 202
 203 =item * C<gsl_fft_real_workspace_free($workspace)> - This function frees the memory associated with the workspace $workspace. The workspace can be freed if no further FFTs of the same length will be needed.
 204
 205 =item * C<gsl_fft_real_transform >
 206
 207 =item * C<gsl_fft_real_unpack >
 208
 209 =back
 210
 211 This module also includes the following constants :
 212
 213 =over
 214
 215 =item * C<$gsl_fft_forward>
 216
 217 =item * C<$gsl_fft_backward>
 218
 219 =back
 220
 221 For more informations on the functions, we refer you to the GSL offcial
 222 documentation: L<http://www.gnu.org/software/gsl/manual/html_node/>
 223
 224  Tip : search on google: site:http://www.gnu.org/software/gsl/manual/html_node/ name_of_the_function_you_want
 225
 226 =head1 EXAMPLE
 227
 228 This is a very simple example of how to format the array reference needed by various functions in this module:
 229
 230     use Math::GSL::FFT;
 231     use Math::GSL::Errno qw/:all/;
 232     my $data = [ (1) x 10, (0) x 236, (1) x 10 ];
 233     my ($status, $fft) = gsl_fft_complex_radix2_forward ($data, 1, 128);
 234
 235 =head1 AUTHORS
 236
 237 Jonathan Leto <jonathan@leto.net> and Thierry Moisan <thierry.moisan@gmail.com>
 238
 239 =head1 COPYRIGHT AND LICENSE
 240
 241 Copyright (C) 2008 Jonathan Leto and Thierry Moisan
 242
 243 This program is free software; you can redistribute it and/or modify it
 244 under the same terms as Perl itself.
 245
 246 =cut
 247
 248 %}