swig/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 - Fast Fourier Transforms (FFT)
 109
 110 =head1 SYNOPSIS
 111
 112     use Math::GSL::FFT qw /:all/;
 113     # alternating elements are real/imaginary part, hence 256 element array
 114     my $data = [ (1) x 10, (0) x 236, (1) x 10 ];
 115
 116     # use every element of the array
 117     my $stride = 1;
 118
 119     # But it contains 128 complex numbers
 120     my ($status, $fft) = gsl_fft_complex_radix2_forward ($data, $stride, 128);
 121
 122 =head1 DESCRIPTION
 123
 124 This module and this documentation is still in a very early state. Danger Will Robinson!
 125 An OO interface will evolve soon.
 126
 127 =over
 128
 129 =item * C<gsl_fft_complex_radix2_forward($data, $stride, $n) >
 130
 131 This function computes the forward FFTs of length $n with stride $stride, on
 132 the array reference $data using an in-place radix-2 decimation-in-time
 133 algorithm. The length of the transform $n is restricted to powers of two. For
 134 the transform version of the function the sign argument can be either forward
 135 (-1) or backward (+1). The functions return a value of $GSL_SUCCESS if no
 136 errors were detected, or $GSL_EDOM if the length of the data $n is not a power
 137 of two. The complex functions of the FFT module are not yet fully implemented.
 138
 139 =item * C<gsl_fft_complex_radix2_backward >
 140
 141 =item * C<gsl_fft_complex_radix2_inverse >
 142
 143 =item * C<gsl_fft_complex_radix2_transform >
 144
 145 =item * C<gsl_fft_complex_radix2_dif_forward >
 146
 147 =item * C<gsl_fft_complex_radix2_dif_backward >
 148
 149 =item * C<gsl_fft_complex_radix2_dif_inverse >
 150
 151 =item * C<gsl_fft_complex_radix2_dif_transform >
 152
 153 =item * C<gsl_fft_complex_wavetable_alloc($n)>
 154
 155 This function prepares a trigonometric lookup table for a complex FFT of length
 156 $n. The function returns a pointer to the newly allocated
 157 gsl_fft_complex_wavetable if no errors were detected, and a null pointer in the
 158 case of error. The length $n is factorized into a product of subtransforms, and
 159 the factors and their trigonometric coefficients are stored in the wavetable.
 160 The trigonometric coefficients are computed using direct calls to sin and cos,
 161 for accuracy. Recursion relations could be used to compute the lookup table
 162 faster, but if an application performs many FFTs of the same length then this
 163 computation is a one-off overhead which does not affect the final throughput.
 164 The wavetable structure can be used repeatedly for any transform of the same
 165 length. The table is not modified by calls to any of the other FFT functions.
 166 The same wavetable can be used for both forward and backward (or inverse)
 167 transforms of a given length.
 168
 169 =item * C<gsl_fft_complex_wavetable_free($wavetable)>
 170
 171 This function frees the memory associated with the wavetable $wavetable. The
 172 wavetable can be freed if no further FFTs of the same length will be needed.
 173
 174 =item * C<gsl_fft_complex_workspace_alloc($n)>
 175
 176 This function allocates a workspace for a complex transform of length $n.
 177
 178 =item * C<gsl_fft_complex_workspace_free($workspace) >
 179
 180 This function frees the memory associated with the workspace $workspace. The
 181 workspace can be freed if no further FFTs of the same length will be needed.
 182
 183 =item * C<gsl_fft_complex_memcpy >
 184
 185 =item * C<gsl_fft_complex_forward >
 186
 187 =item * C<gsl_fft_complex_backward >
 188
 189 =item * C<gsl_fft_complex_inverse >
 190
 191 =item * C<gsl_fft_complex_transform >
 192
 193 =item * C<gsl_fft_halfcomplex_radix2_backward($data, $stride, $n)>
 194
 195 This function computes the backwards in-place radix-2 FFT of length $n and
 196 stride $stride on the half-complex sequence data stored according the output
 197 scheme used by gsl_fft_real_radix2. The result is a real array stored in
 198 natural order.
 199
 200 =item * C<gsl_fft_halfcomplex_radix2_inverse($data, $stride, $n)>
 201
 202 This function computes the inverse in-place radix-2 FFT of length $n and stride
 203 $stride on the half-complex sequence data stored according the output scheme
 204 used by gsl_fft_real_radix2. The result is a real array stored in natural
 205 order.
 206
 207 =item * C<gsl_fft_halfcomplex_radix2_transform>
 208
 209 =item * C<gsl_fft_halfcomplex_wavetable_alloc($n)>
 210
 211 This function prepares trigonometric lookup tables for an FFT of size $n real
 212 elements. The functions return a pointer to the newly allocated struct if no
 213 errors were detected, and a null pointer in the case of error. The length $n is
 214 factorized into a product of subtransforms, and the factors and their
 215 trigonometric coefficients are stored in the wavetable. The trigonometric
 216 coefficients are computed using direct calls to sin and cos, for accuracy.
 217 Recursion relations could be used to compute the lookup table faster, but if an
 218 application performs many FFTs of the same length then computing the wavetable
 219 is a one-off overhead which does not affect the final throughput.  The
 220 wavetable structure can be used repeatedly for any transform of the same
 221 length. The table is not modified by calls to any of the other FFT functions.
 222 The appropriate type of wavetable must be used for forward real or inverse
 223 half-complex transforms.
 224
 225 =item * C<gsl_fft_halfcomplex_wavetable_free($wavetable)>
 226
 227 This function frees the memory associated with the wavetable $wavetable. The
 228 wavetable can be freed if no further FFTs of the same length will be needed.
 229
 230 =item * C<gsl_fft_halfcomplex_backward >
 231
 232 =item * C<gsl_fft_halfcomplex_inverse >
 233
 234 =item * C<gsl_fft_halfcomplex_transform >
 235
 236 =item * C<gsl_fft_halfcomplex_unpack >
 237
 238 =item * C<gsl_fft_halfcomplex_radix2_unpack >
 239
 240 =item * C<gsl_fft_real_radix2_transform($data, $stride, $n) >
 241
 242 This function computes an in-place radix-2 FFT of length $n and stride $stride
 243 on the real array reference $data. The output is a half-complex sequence, which
 244 is stored in-place. The arrangement of the half-complex terms uses the
 245 following scheme: for k < N/2 the real part of the k-th term is stored in
 246 location k, and the corresponding imaginary part is stored in location N-k.
 247 Terms with k > N/2 can be reconstructed using the symmetry z_k = z^*_{N-k}. The
 248 terms for k=0 and k=N/2 are both purely real, and count as a special case.
 249 Their real parts are stored in locations 0 and N/2 respectively, while their
 250 imaginary parts which are zero are not stored. The following table shows the
 251 correspondence between the output data and the equivalent results obtained by
 252 considering the input data as a complex sequence with zero imaginary part,
 253
 254           complex[0].real    =    data[0]
 255           complex[0].imag    =    0
 256           complex[1].real    =    data[1]
 257           complex[1].imag    =    data[N-1]
 258           ...............         ................
 259           complex[k].real    =    data[k]
 260           complex[k].imag    =    data[N-k]
 261           ...............         ................
 262           complex[N/2].real  =    data[N/2]
 263           complex[N/2].imag  =    0
 264           ...............         ................
 265           complex[k'].real   =    data[k]        k' = N - k
 266           complex[k'].imag   =   -data[N-k]
 267           ...............         ................
 268           complex[N-1].real  =    data[1]
 269           complex[N-1].imag  =   -data[N-1]
 270
 271 =for notyou #'
 272
 273 Note that the output data can be converted into the full complex sequence using
 274 the function gsl_fft_halfcomplex_unpack.
 275
 276 =item * C<gsl_fft_real_wavetable_alloc($n)>
 277
 278 This function prepares trigonometric lookup tables for an FFT of size $n real
 279 elements. The functions return a pointer to the newly allocated struct if no
 280 errors were detected, and a null pointer in the case of error. The length $n is
 281 factorized into a product of subtransforms, and the factors and their
 282 trigonometric coefficients are stored in the wavetable. The trigonometric
 283 coefficients are computed using direct calls to sin and cos, for accuracy.
 284 Recursion relations could be used to compute the lookup table faster, but if an
 285 application performs many FFTs of the same length then computing the wavetable
 286 is a one-off overhead which does not affect the final throughput.  The
 287 wavetable structure can be used repeatedly for any transform of the same
 288 length. The table is not modified by calls to any of the other FFT functions.
 289 The appropriate type of wavetable must be used for forward real or inverse
 290 half-complex transforms.
 291
 292 =item * C<gsl_fft_real_wavetable_free($wavetable)>
 293
 294 This function frees the memory associated with the wavetable $wavetable. The
 295 wavetable can be freed if no further FFTs of the same length will be needed.
 296
 297 =item * C<gsl_fft_real_workspace_alloc($n)>
 298
 299 This function allocates a workspace for a real transform of length $n. The same
 300 workspace can be used for both forward real and inverse halfcomplex transforms.
 301
 302 =item * C<gsl_fft_real_workspace_free($workspace)>
 303
 304 This function frees the memory associated with the workspace $workspace. The
 305 workspace can be freed if no further FFTs of the same length will be needed.
 306
 307 =item * C<gsl_fft_real_transform >
 308
 309 =item * C<gsl_fft_real_unpack >
 310
 311 =back
 312
 313 This module also includes the following constants :
 314
 315 =over
 316
 317 =item * C<$gsl_fft_forward>
 318
 319 =item * C<$gsl_fft_backward>
 320
 321 =back
 322
 323 For more informations on the functions, we refer you to the GSL offcial
 324 documentation: L<http://www.gnu.org/software/gsl/manual/html_node/>
 325
 326  Tip : search on google: site:http://www.gnu.org/software/gsl/manual/html_node/ name_of_the_function_you_want
 327
 328
 329 =head1 AUTHORS
 330
 331 Jonathan Leto <jonathan@leto.net> and Thierry Moisan <thierry.moisan@gmail.com>
 332
 333 =head1 COPYRIGHT AND LICENSE
 334
 335 Copyright (C) 2008 Jonathan Leto and Thierry Moisan
 336
 337 This program is free software; you can redistribute it and/or modify it
 338 under the same terms as Perl itself.
 339
 340 =cut
 341
 342 %}