Christian Soeller > PDL-2.4.1 > PDL::FFTW

Download:
PDL-2.4.1.tar.gz

Annotate this POD

CPAN RT

Open  0
View/Report Bugs
Source   Latest Release: PDL-2.006

NAME ^

PDL::FFTW - PDL interface to the Fastest Fourier Transform in the West v2.x

DESCRIPTION ^

This is a means to interface PDL with the FFTW library. It's similar to the standard FFT routine but it's usually faster and has support for real transforms. It works well for the types of PDLs for which was the library was compiled (otherwise it must do conversions).

SYNOPSIS ^

    use PDL::FFTW

    load_wisdom("file_name");

    out_cplx_pdl =  fftw(in_cplx_pdl);
    out_cplx_pdl = ifftw(in_cplx_pdl);

    out_cplx_pdl =  rfftw(in_real_pdl);
    out_real_pdl = irfftw(in_cplx_pdl);

    cplx_pdl = nfftw(cplx_pdl);
    cplx_pdl = infftw(cplx_pdl);

    cplx_pdl = Cmul(a_cplx_pdl, b_cplx_pdl);
    cplx_pdl = Cconj(a_cplx_pdl);
    real_pdl = Cmod(a_cplx_pdl); 
    real_pdl = Cmod2(a_cplx_pdl); 

FFTW documentation ^

Please refer to the FFTW documentation for a better understanding of these functions.

Note that complex numbers are represented as piddles with leading dimension size 2 (real/imaginary pairs).

load_wisdom

Loads the wisdom from a file for better FFTW performance.

The wisdom is automatically saved when the program ends. It will be automagically called when the variable $PDL::FFT::wisdom is set to a file name. For example, the following is a useful idiom to have in your .perldlrc file:

  $PDL::FFT::wisdom = "$ENV{HOME}/.fftwisdom"; # save fftw wisdom in this file

Explicit usage:

  load_wisdom($fname);

fftw

calculate the complex FFT of a real piddle (complex input, complex output)

   $pdl_cplx = fftw $pdl_cplx;

ifftw

Complex inverse FFT (complex input, complex output).

   $pdl_cplx = ifftw $pdl_cplx;

nfftw

Complex inplace FFT (complex input, complex output).

   $pdl_cplx = nfftw $pdl_cplx;

infftw

Complex inplace inverse FFT (complex input, complex output).

   $pdl_cplx = infftw $pdl_cplx;

rfftw

Real FFT. For an input piddle of dimensions [n1,n2,...] the output is [2,(n1/2)+1,n2,...] (real input, complex output).

  $pdl_cplx = fftw $pdl_real;

irfftw

Real inverse FFT. Have a look at rfftw to understand the format. USE ONLY an even n1! (complex input, real output)

  $pdl_real = ifftw $pdl_cplx;

nrfftw

Real inplace FFT. If you want a transformation on a piddle with dimensions [n1,n2,....] you MUST pass in a piddle with dimensions [2*(n1/2+1),n2,...] (real input, complex output).

Use with care due to dimension restrictions mentioned below. For details check the html docs that come with the fftw library.

  $pdl_cplx = nrfftw $pdl_real;

inrfftw

Real inplace inverse FFT. Have a look at nrfftw to understand the format. USE ONLY an even first dimension size! (complex input, real output)

  $pdl_real = infftw $pdl_cplx;

rfftwconv

ND convolution using real ffts from the FFTW library

  $conv = rfftwconv $im, kernctr $im, $k; # kernctr is from PDL::FFT

fftwconv

ND convolution using complex ffts from the FFTW library

Assumes real input!

  $conv = fftwconv $im, kernctr $im, $k; # kernctr is from PDL::FFT

Complex multiplication

   $out_pdl_cplx = Cmul($a_pdl_cplx,$b_pdl_cplx);

EOD );

pp_def('Cscale', Pars => 'a(n); b(); [o]c(n);', Code => ' if ($SIZE(n)!=2) barf("This function works only on complex\n"); threadloop %{ $c(n => 0)= $a(n => 0)*$b(); $c(n => 1)= $a(n => 1)*$b(); %} ', Doc => <<'EOD', =for ref

Complex by real multiplation.

  Cscale($a_pdl_cplx,$b_pdl_real);
EOD
      );

pp_def('Cdiv', Pars => 'a(n); b(n); [o]c(n);', Code => ' $GENERIC() divi;

if ($SIZE(n)!=2) barf("This function works only on complex\n"); threadloop %{ divi = $b(n => 0)*$b(n => 0) + $b(n => 1)*$b(n => 1);

$c(n => 0)= ( $a(n => 0)*$b(n => 0) + $a(n => 1)*$b(n => 1) ) / divi; $c(n => 1)= ( $a(n => 1)*$b(n => 0) - $a(n => 0)*$b(n => 1) ) / divi; %} ', Doc => <<'EOD', =for ref

Complex division.

  $out_pdl_cplx = Cdiv($a_pdl_cplx,$b_pdl_cplx);

EOD );

pp_def('Cbmul', Pars => 'a(n); b(n);', Code => ' $GENERIC() tmp1,tmp;

if ($SIZE(n)!=2) barf("This function works only on complex\n"); threadloop %{ tmp = $a(n => 0); tmp1 = $a(n => 1); $a(n => 0)= tmp * $b(n => 0) - tmp1 * $b(n => 1); $a(n => 1)= tmp * $b(n => 1) + tmp1 * $b(n => 0); %} ', Doc => <<'EOD', =for ref

Complex inplace multiplication.

   Cbmul($a_pdl_cplx,$b_pdl_cplx);
EOD
      );

pp_def('Cbscale', Pars => 'a(n); b();', Code => ' $GENERIC() tmp1,tmp,divi;

if ($SIZE(n)!=2) barf("This function works only on complex\n"); threadloop %{ $a(n => 0) *= $b(); $a(n => 1) *= $b(); %} ', Doc => <<'EOD', =for ref

Complex inplace multiplaction by real.

   Cbscale($a_pdl_cplx,$b_pdl_real);
EOD
      );

pp_def('Cconj', Pars => 'a(n); [o]c(n);', Code => ' if ($SIZE(n)!=2) barf("This function works only on complex\n"); threadloop %{ $c(n => 0)= $a(n => 0); $c(n => 1)= -$a(n => 1); %} ', Doc => <<'EOD', =for ref

Complex conjugate.

   $out_pdl_cplx = Cconj($a_pdl_cplx);
EOD
      );

pp_def('Cbconj', Pars => 'a(n);', Code => ' if ($SIZE(n)!=2) barf("This function works only on complex\n"); threadloop %{ $a(n => 1)= -$a(n => 1); %} ', Doc => <<'EOD', =for ref

Complex inplace conjugate.

   Cbconj($a_pdl_cplx);
EOD
      );

pp_def('Cexp', Pars => 'a(n); [o]c(n);', Code => ' if ($SIZE(n)!=2) barf("This function works only on complex\n"); threadloop %{ $c(n => 0)= exp($a(n => 0)) * cos($a(n => 1)); $c(n => 1)= exp($a(n => 0)) * sin($a(n => 1)); %} ', Doc => <<'EOD', =for ref

Complex exponentation.

   $out_pdl_cplx = Cexp($a_pdl_cplx);
EOD
      );

pp_def('Cbexp', Pars => 'a(n);', Code => ' if ($SIZE(n)!=2) barf("This function works only on complex\n"); threadloop %{ double re = $a(n => 0); double im = $a(n => 1);

$a(n => 0)= exp(re) * cos(im); $a(n => 1)= exp(re) * sin(im); %} ', Doc => <<'EOD', =for ref

Complex inplace exponentation.

   $out_pdl_cplx = Cexp($a_pdl_cplx);
EOD
      );

pp_def('Cmod', Pars => 'a(n); [o]c();', Code => ' if ($SIZE(n)!=2) barf("This function works only on complex\n"); threadloop %{ $c()= sqrt ( $a(n => 0)*$a(n => 0) + $a(n => 1)*$a(n => 1) ); %} ', Doc => <<'EOD', =for ref

modulus of a complex piddle.

  $out_pdl_real = Cmod($a_pdl_cplx);
EOD
      );

pp_def('Carg', Pars => 'a(n); [o]c();', Code => ' if ($SIZE(n)!=2) barf("This function works only on complex\n"); threadloop %{ $c()= atan2($a(n=>1),$a(n=>0)); %} ', Doc => <<'EOD', =for ref

argument of a complex number.

   $out_pdl_real = Carg($a_pdl_cplx);
EOD
      );

pp_def('Cmod2', Pars => 'a(n); [o]c();', Code => ' if ($SIZE(n)!=2) barf("This function works only on complex\n"); threadloop %{ $c()= ( $a(n => 0)*$a(n => 0) + $a(n => 1)*$a(n => 1) ); %} ', Doc => <<'EOD', =for ref

Returns squared modulus of a complex number.

   $out_pdl_real = Cmod2($a_pdl_cplx);
EOD
      );

### real fftw

pp_addxs('',' MODULE = PDL::FFTW PACKAGE = PDL::FFTW

int PDL_rfftwnd_create_plan(dims, dir, flag) pdl* dims int dir int flag CODE: { fftw_direction fdir=0; int fflag=FFTW_USE_WISDOM;

  if (dims->ndims != 1) {barf("Only 1d input dimesions make sense");} 
  if (dims->data == NULL) {barf("input piddles must be physical");} 
  if (dims->datatype != PDL_L) {barf("Only integers please");} 

  if (dir) {
    fdir=FFTW_COMPLEX_TO_REAL;
  }
  else {
    fdir=FFTW_REAL_TO_COMPLEX;
  }
  if (flag & 1 ) { 
    fflag |= FFTW_ESTIMATE;
  }
  else {
    fflag |= FFTW_MEASURE;
  }
  if (flag & 2 ) { 
    fflag |= FFTW_IN_PLACE;
  }
  else {
    fflag |= FFTW_OUT_OF_PLACE;
  }
 
  RETVAL = 
  (int) rfftwnd_create_plan( dims->dims[0], 
                            ( int *) dims->data,
                            fdir,
                            fflag);

 }
OUTPUT:
 RETVAL
'
);

pp_addxs('',' MODULE = PDL::FFTW PACKAGE = PDL::FFTW

void PDL_rfftwnd_one_real_to_complex(plan, in, out) int plan pdl* in pdl* out CODE: if (in->data==NULL || out->data==NULL) {barf("Need a physical pdl!");} if (in->datatype != PDL_MYTYPE || out->datatype != PDL_MYTYPE) {barf("Bad Type");} rfftwnd_one_real_to_complex( (rfftwnd_plan) plan, (fftw_real *) in->data, (fftw_complex *) out->data); ');

pp_addxs('',' MODULE = PDL::FFTW PACKAGE = PDL::FFTW

void PDL_rfftwnd_one_complex_to_real(plan, in, out) int plan pdl* in pdl* out CODE: if (in->data==NULL || out->data==NULL) {barf("Need a physical pdl!");} if (in->datatype != PDL_MYTYPE || out->datatype != PDL_MYTYPE) {barf("Bad type");} rfftwnd_one_complex_to_real( (rfftwnd_plan) plan, (fftw_complex *) in->data, (fftw_real *) out->data);

');

# NOTE: BUG! the inplace code below will not work with slices! # no backpropagation of results! pp_addxs('',' MODULE = PDL::FFTW PACKAGE = PDL::FFTW

void PDL_inplace_rfftwnd_one_real_to_complex(plan, in) int plan pdl* in CODE: if (in->data==NULL) {barf("Need a physical pdl!");} if (in->datatype != PDL_MYTYPE) {barf("Bad Type");} PDL->children_changesoon(in, PDL_PARENTDATACHANGED); rfftwnd_one_real_to_complex( (rfftwnd_plan) plan, (fftw_real *) in->data, NULL); /* this call is crucial to propagate changes back if slices are given as arguments * Note: must not used vaffinechanged (!) since any slice has physical data */ PDL->changed( in , PDL_PARENTDATACHANGED , 0 ); ');

pp_addxs('',' MODULE = PDL::FFTW PACKAGE = PDL::FFTW

void PDL_inplace_rfftwnd_one_complex_to_real(plan, in) int plan pdl* in CODE: if (in->data==NULL) {barf("Need a physical pdl!");} if (in->datatype != PDL_MYTYPE) {barf("Bad type");} PDL->children_changesoon(in, PDL_PARENTDATACHANGED); rfftwnd_one_complex_to_real( (rfftwnd_plan) plan, (fftw_complex *) in->data, NULL); /* this call is crucial to propagate changes back if slices are given as arguments * Note: must not used vaffinechanged (!) since any slice has physical data */ PDL->changed( in , PDL_PARENTDATACHANGED , 0 );

');

### complex fftw

pp_addxs('',' MODULE = PDL::FFTW PACKAGE = PDL::FFTW

int PDL_fftwnd_create_plan(dims, dir, flag) pdl* dims int dir int flag CODE: { fftw_direction fdir=0; int fflag=FFTW_USE_WISDOM;

  if (dims->ndims != 1) {barf("Only 1d input dimesions make sense");} 
  if (dims->data == NULL) {barf("input piddles must be physical");} 
  if (dims->datatype != PDL_L) {barf("Only integers please");} 

  if (dir) {
    fdir=FFTW_BACKWARD;
  }
  else {
    fdir=FFTW_FORWARD;
  }
  if (flag & 1 ) { 
    fflag |= FFTW_ESTIMATE;
  }
  else {
    fflag |= FFTW_MEASURE;
  }
  if (flag & 2 ) { 
    fflag |= FFTW_IN_PLACE;
  }
  else {
    fflag |= FFTW_OUT_OF_PLACE;
  }
 
  RETVAL = 
  (int) fftwnd_create_plan( dims->dims[0], 
                            ( int *) dims->data,
                            fdir,
                            fflag);

 }
OUTPUT:
 RETVAL
'
);

pp_addxs('',' MODULE = PDL::FFTW PACKAGE = PDL::FFTW

void PDL_fftwnd_one(plan, in, out) int plan pdl* in pdl* out CODE: if (in->data==NULL || out->data==NULL) {barf("Need a physical pdl!");} if (in->datatype != PDL_MYTYPE || out->datatype != PDL_MYTYPE) {barf("Bad type!");} fftwnd_one( (fftwnd_plan) plan, (fftw_complex *) in->data, (fftw_complex *) out->data); ');

pp_addxs('',' MODULE = PDL::FFTW PACKAGE = PDL::FFTW

void PDL_inplace_fftwnd_one(plan, in) int plan pdl* in CODE: if (in->data==NULL) {barf("Need a physical pdl!");} if (in->datatype != PDL_MYTYPE) {barf("Only float please");} PDL->children_changesoon(in, PDL_PARENTDATACHANGED); fftwnd_one( (fftwnd_plan) plan, (fftw_complex *) in->data, NULL); /* this call is crucial to propagate changes back if slices are given as arguments * Note: must not used vaffinechanged (!) since any slice has physical data */ PDL->changed( in , PDL_PARENTDATACHANGED , 0 ); ');

### wisdom stuff

pp_addxs('',' MODULE = PDL::FFTW PACKAGE = PDL::FFTW

int PDL_fftw_import_wisdom_from_string (wisdom) char* wisdom CODE: RETVAL = ( fftw_import_wisdom_from_string(wisdom) == FFTW_SUCCESS ); OUTPUT: RETVAL ');

pp_addxs('',' MODULE = PDL::FFTW PACKAGE = PDL::FFTW

char* PDL_fftw_export_wisdom_to_string () CODE: RETVAL = fftw_export_wisdom_to_string(); OUTPUT: RETVAL ');

pp_add_exported('','load_wisdom save_wisdom rfftw irfftw fftw ifftw nfftw infftw nrfftw inrfftw fftwconv rfftwconv kernctr');

# I don't see a point in exporting these (CS) # 'PDL_rfftwnd_create_plan PDL_rfftwnd_one_real_to_complex PDL_rfftwnd_one_complex_to_real PDL_fftw_export_wisdom_to_string PDL_fftw_import_wisdom_from_string PDL_inplace_fftwnd_one PDL_fftwnd_one PDL_fftwnd_create_plan PDL_inplace_rfftwnd_one_real_to_complex PDL_inplace_rfftwnd_one_complex_to_real';

pp_addpm({At => 'Bot'},<< 'EOD'); =head1 AUTHOR

Copyright (C) 1999 Christian Pellegrin, 2000 Christian Soeller. All rights reserved. There is no warranty. You are allowed to redistribute this software / documentation under certain conditions. For details, see the file COPYING in the PDL distribution. If this file is separated from the PDL distribution, the copyright notice should be included in the file.

syntax highlighting: