FFTW3.pd

=head1 NAME

PDL::FFTW3 - PDL interface to the Fastest Fourier Transform in the West v3

=cut

# -*- cperl -*-

##### General layout of the module #####
#
# Each type of transform that is supported by this module has a plain,
# unthreaded perl entry point the user calls. This entry point makes sure the
# FFTW plan exists (or makes it). Then it calls the THREADED PP function to
# actually compute the transform

use strict;
use warnings;

# I generate code for up to 10-dimensional FFTs
my $maxrank = 10;

our $VERSION = '0.19';
pp_setversion($VERSION);

pp_addpm( {At => 'Top'}, slurp('README.pod') . <<'EOF' );

use strict;
use warnings;

EOF

pp_addhdr( '
#include <fftw3.h>
' );


# I want to be able to say $X = fft1($x); rank is required. 'fft()' is ambiguous
# about whether threading is desired or if a large fft is desired. Old PDL::FFTW
# did one thing, matlab does another, so I do not include this function at all




# I define up to rank-10 FFTs. This is annoyingly arbitrary, but hopefully
# should be sufficient
for my $rank (1..$maxrank)
{
  generateDefinitions($rank);
}
pp_export_nothing();

pp_addxs('', slurp('compute_plan_template.xs'));
pp_addpm( {At => 'Middle'}, slurp('FFTW3_header_include.pm') );

for my $rank (1..$maxrank)
{
  my $shapestr = sprintf(q{$a->shape->slice('1:%d')->prodover},$rank);
  my $rshapestr = sprintf(q{$a->shape->slice('0:%d')->prodover},$rank-1);

  pp_addpm({At => 'Bot'}, pp_line_numbers(__LINE__, <<EOF));
sub fft$rank { __fft_internal( "fft$rank",\@_ ); }
*PDL::fft$rank = \\&fft$rank;

sub ifft$rank {
  my \$a = __fft_internal( "ifft$rank", \@_ );
  \$a /= \$_[0]->type->real ? $shapestr : $rshapestr;
  \$a;
}
*PDL::ifft$rank = \\&ifft$rank;

sub rfft$rank { __fft_internal( "rfft$rank", \@_ ); }
*PDL::rfft$rank = \\&rfft$rank;

sub rNfft$rank { __fft_internal( "rNfft$rank", \@_ ); }
*PDL::rNfft$rank = \\&rNfft$rank;

sub irfft$rank { my \$a = __fft_internal( "irfft$rank", \@_ ); \$a /= $rshapestr; \$a; }
*PDL::irfft$rank = \\&irfft$rank;
EOF

  pp_add_exported( map "${_}fft$rank", '', 'i', 'r', 'rN', 'ir' );
}


##########
# Generate the fftn case.  This should probably be done more prettily; for now it's just 
# a springboard that jumps into __fft_internal.
pp_addpm({At=> 'Bot'}, pp_line_numbers(__LINE__, sprintf <<'EOF', $maxrank));
sub _rank_springboard {
  my ($name, $source, $rank, @rest) = @_;
  my $inverse = ($name =~ m/^i/);
  my $real    = ($name =~ m/r/) || !$source->type->real;

  unless(defined $rank) {
    die "${name}n: second argument must be the rank of the transform you want";
  }
  $rank = 0+$rank;  # force numeric context
  unless($rank>=1 ) {
    die "${name}n: second argument (rank) must be between 1 and %d";
  }

  my $active_lo = ($real ? 0 : 1);
  my $active_hi = ($real ? $rank-1 : $rank);

  unless($source->ndims > $active_hi) {
    die "${name}n: rank is $rank but input has only ".($active_hi-$active_lo)." active dims!";
  }

  my $out = __fft_internal( $name.$rank, $source, @rest );

  if($inverse) {
    $out /= $out->shape->slice("$active_lo:$active_hi")->prodover;
  }
  return $out;
}

sub fftn    { _rank_springboard( "fft",      @_ ) }
sub ifftn   { _rank_springboard( "ifft",     @_ ) }
sub rfftn   { _rank_springboard( "rfft",  @_ ) }
sub irfftn  { _rank_springboard( "irfft", @_ ) }

*PDL::fftn   = \&fftn;
*PDL::ifftn  = \&ifftn;
*PDL::rfftn  = \&rfftn;
*PDL::irfftn = \&irfftn;

EOF
pp_add_exported( map "${_}fftn", '','i','r','ir' );

pp_done();


sub generateDefinitions
{
  my $rank = shift;

  my %pp_def = (
    HandleBad => 0,
    GenericTypes => [qw(F D)],
    OtherPars    => 'IV plan', # comes from pre-fft code not user
    # this is a private function so I don't want to create
    # user-visible documentation or exports
    Doc          => undef,
    PMFunc       => ''
  );

  ################################################################################
  ####### first I generate the definitions for the simple complex-complex FFT case
  # make dimension string 'n0=2,n1,n2,n3,n4...'. The leading 2 is for the
  # (real,imag) complex pair
  my @dims = map "n$_", 1..$rank;
  unshift @dims, 'n0=2';
  my $dims_string = join(',', @dims);
  $pp_def{Pars} = "in($dims_string); [o]out($dims_string);";
  $pp_def{Code} = slurp('template_complex.c');
  pp_def( "__fft$rank", %pp_def );

  ##################################################################################
  ####### now I generate the definitions for the real-complex and complex-real cases
  my @dims_real    = @dims;
  my @dims_complex = @dims;
  shift @dims_real; # get rid of the (real,imag) dimension for the real numbers
  $dims_complex[1] = 'nhalf'; # first complex dim is real->dim(0)/2+1
  my $dims_real_string    = join(',', @dims_real);
  my $dims_complex_string = join(',', @dims_complex);
  my $code_real          =  slurp('template_real.c');
  $code_real             =~ s/RANK/$rank/ge;
  my $code_real_forward  =  $code_real;
  my $code_real_backward =  $code_real;
  $code_real_forward  =~ s/INVERSE/0/g;
  $code_real_backward =~ s/INVERSE/1/g;

  # forward
  # I have the real dimensions, but not nhalf
  $pp_def{RedoDimsCode} = <<'EOF';
if( $PDL(complexv)->ndims <= 1 || $PDL(complexv)->dims[1] <= 0 )
  $SIZE(nhalf) = (int)( $PDL(real)->dims[0]/2 ) + 1;
EOF
  $pp_def{Pars} = "real($dims_real_string); [o]complexv($dims_complex_string);";
  $pp_def{Code} = $code_real_forward;
  pp_def( "__rfft$rank", %pp_def );

  # backward
  # I have the complex dimensions. Have nhalf, but not n1
  #
  # if we're not given an output, there's an ambiguity. I want
  # int($out->dim(0)/2) + 1 != $in->dim(1),
  # however this could mean that
  #  $out->dim(0) = 2*$in->dim(1) - 2
  # or
  #  $out->dim(0) = 2*$in->dim(1) - 1
  #
  # WITHOUT ANY OTHER INFORMATION, I ASSUME EVEN INPUT SIZES, SO I ASSUME
  #  $out->dim(0) = 2*$in->dim(1) - 2
  $pp_def{RedoDimsCode} = <<'EOF';
if( $PDL(real)->dims[0] <= 0 )
  $SIZE(n1) = 2*$PDL(complexv)->dims[1] - 2;
EOF
  $pp_def{Pars} = "complexv($dims_complex_string); [o]real($dims_real_string);";
  $pp_def{Code} = $code_real_backward;
  pp_def( "__irfft$rank", %pp_def );

  ################################################################################
  ####### now native-complex version
  shift @dims; # drop the (real,imag) dim
  $dims_string = join(',', @dims);
  delete $pp_def{RedoDimsCode};
  $pp_def{Pars} = "in($dims_string); [o]out($dims_string);";
  $pp_def{Code} = slurp('template_complex.c');
  $pp_def{Code} =~ s/TFD/TGC/g;
  $pp_def{Code} =~ s/\*2//g; # the sizeof-doubling
  $pp_def{GenericTypes} = [qw(G C)];
  pp_def( "__Nfft$rank", %pp_def );

  ##################################################################################
  ####### real-native complex and native complex-real
  @dims_real    = @dims;
  @dims_complex = @dims;
  $dims_complex[0] = 'nhalf'; # first complex dim is real->dim(0)/2+1
  $dims_real_string    = join(',', @dims_real);
  $dims_complex_string = join(',', @dims_complex);
  $code_real          =  slurp('template_real.c');
  $code_real             =~ s/RANK/$rank-1/ge;
  $code_real_forward  =  $code_real;
  $code_real_backward =  $code_real;
  $code_real_forward  =~ s/INVERSE/0/g;
  $code_real_backward =~ s/INVERSE/1/g;

  # forward
  $pp_def{RedoDimsCode} = <<'EOF';
if( $PDL(complexv)->ndims <= 1 || $PDL(complexv)->dims[1] <= 0 )
  $SIZE(nhalf) = (int)( $PDL(real)->dims[0]/2 ) + 1;
EOF
  $pp_def{Pars} = "real($dims_real_string); complex [o]complexv($dims_complex_string);";
  $pp_def{Code} = $code_real_forward;
  $pp_def{GenericTypes} = [qw(F D)];
  pp_def( "__rNfft$rank", %pp_def );

  # backward
  $pp_def{RedoDimsCode} = <<'EOF';
if( $PDL(real)->dims[0] <= 0 )
  $SIZE(n1) = 2*$PDL(complexv)->dims[0] - 2;
EOF
  $pp_def{Pars} = "complexv($dims_complex_string); real [o]real($dims_real_string);";
  $code_real_backward =~ s/TFD/TGC/g;
  $code_real_backward =~ s/\*2//g; # the sizeof-doubling
  $pp_def{Code} = $code_real_backward;
  $pp_def{GenericTypes} = [qw(G C)];
  pp_def( "__irNfft$rank", %pp_def );
}

sub slurp
{
  my $filename = shift;
  open FD, '<', $filename or die "Couldn't open '$filename' for rading";

  local $/ = undef;
  my $contents = <FD>;
  close FD;
  return qq{\n#line 0 "$filename"\n\n} . $contents;
}