"Fossies" - the Fresh Open Source Software Archive

Member "PDL-2.080/GENERATED/PDL/FFT.pm" (28 May 2022, 10387 Bytes) of package /linux/misc/PDL-2.080.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) Perl source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "FFT.pm" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 2.079_vs_2.080.

    1 #
    2 # GENERATED WITH PDL::PP! Don't modify!
    3 #
    4 package PDL::FFT;
    5 
    6 our @EXPORT_OK = qw(fft ifft fftnd ifftnd fftconvolve realfft realifft kernctr );
    7 our %EXPORT_TAGS = (Func=>\@EXPORT_OK);
    8 
    9 use PDL::Core;
   10 use PDL::Exporter;
   11 use DynaLoader;
   12 
   13 
   14    
   15    our @ISA = ( 'PDL::Exporter','DynaLoader' );
   16    push @PDL::Core::PP, __PACKAGE__;
   17    bootstrap PDL::FFT ;
   18 
   19 
   20 
   21 
   22 
   23 
   24 #line 6 "fft.pd"
   25 
   26 =head1 NAME
   27 
   28 PDL::FFT - FFTs for PDL
   29 
   30 =head1 DESCRIPTION
   31 
   32 !!!!!!!!!!!!!!!!!!!!!!!!!!WARNING!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   33 As of PDL-2.006_04, the direction of the FFT/IFFT has been
   34 reversed to match the usage in the FFTW library and the convention
   35 in use generally.
   36 !!!!!!!!!!!!!!!!!!!!!!!!!!WARNING!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   37 
   38 FFTs for PDL.  These work for arrays of any dimension, although ones
   39 with small prime factors are likely to be the quickest.  The forward
   40 FFT is unnormalized while the inverse FFT is normalized so that the
   41 IFFT of the FFT returns the original values.
   42 
   43 For historical reasons, these routines work in-place and do not recognize
   44 the in-place flag.  That should be fixed.
   45 
   46 =head1 SYNOPSIS
   47 
   48         use PDL::FFT qw/:Func/;
   49 
   50     fft($real, $imag);
   51     ifft($real, $imag);
   52     realfft($real);
   53     realifft($real);
   54 
   55     fftnd($real,$imag);
   56     ifftnd($real,$imag);
   57 
   58     $kernel = kernctr($image,$smallk);
   59     fftconvolve($image,$kernel);
   60 
   61 =head1 DATA TYPES
   62 
   63 The underlying C library upon which this module is based performs FFTs
   64 on both single precision and double precision floating point ndarrays.
   65 The PP functions are defined to only take those data types.
   66 Therefore, if you pass in an ndarray of integer datatype (byte, short,
   67 ushort, long) to any of the routines in PDL::FFT, your data will be
   68 promoted to a double-precision ndarray.  If you pass in a float, the
   69 single-precision FFT will be performed.
   70 
   71 =head1 FREQUENCIES
   72 
   73 For even-sized input arrays, the frequencies are packed like normal
   74 for FFTs (where N is the size of the array and D is the physical step
   75 size between elements):
   76 
   77  0, 1/ND, 2/ND, ..., (N/2-1)/ND, 1/2D, -(N/2-1)/ND, ..., -1/ND.
   78 
   79 which can easily be obtained (taking the Nyquist frequency to be
   80 positive) using
   81 
   82 C<< $kx = $real->xlinvals(-($N/2-1)/$N/$D,1/2/$D)->rotate(-($N/2 -1)); >>
   83 
   84 For odd-sized input arrays the Nyquist frequency is not directly
   85 acessible, and the frequencies are
   86 
   87  0, 1/ND, 2/ND, ..., (N/2-0.5)/ND, -(N/2-0.5)/ND, ..., -1/ND.
   88 
   89 which can easily be obtained using
   90 
   91 C<< $kx = $real->xlinvals(-($N/2-0.5)/$N/$D,($N/2-0.5)/$N/$D)->rotate(-($N-1)/2); >>
   92 
   93 
   94 =head1 ALTERNATIVE FFT PACKAGES
   95 
   96 Various other modules - such as L<PDL::FFTW3> and L<PDL::Slatec> -
   97 contain FFT routines.
   98 However, unlike PDL::FFT, these modules are optional,
   99 and so may not be installed.
  100 
  101 =cut
  102 #line 103 "FFT.pm"
  103 
  104 
  105 
  106 
  107 
  108 
  109 =head1 FUNCTIONS
  110 
  111 =cut
  112 
  113 
  114 
  115 
  116 #line 948 "../../blib/lib/PDL/PP.pm"
  117 
  118 
  119 
  120 =head2 fft
  121 
  122 =for sig
  123 
  124   Signature: ([io]real(n); [io]imag(n))
  125 
  126 =for ref
  127 
  128 Complex 1-D FFT of the "real" and "imag" arrays [inplace]. A single
  129 cfloat/cdouble input ndarray can also be used.
  130 
  131 =for usage
  132 
  133   fft($real,$imag);
  134   fft($complex);
  135 
  136 
  137 =for bad
  138 
  139 fft does not process bad values.
  140 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  141 
  142 
  143 =cut
  144 #line 145 "FFT.pm"
  145 
  146 
  147 
  148 #line 949 "../../blib/lib/PDL/PP.pm"
  149 
  150 sub PDL::fft {
  151     # Convert the first argument to decimal and check for trouble.
  152     my ($re, $im) = @_;
  153     if (!$re->type->real) {
  154         $im=$re->im;
  155         $re=$re->re;
  156     }
  157     eval {  todecimal($re); };
  158     if ($@) {
  159         $@ =~ s/ at .*//s;
  160         barf("Error in FFT with first argument: $@");
  161     }
  162     # Convert the second argument to decimal and check for trouble.
  163     eval {  todecimal($im); };
  164     if ($@) {
  165         $@ =~ s/ at .*//s;
  166         my $message = "Error in FFT with second argument: $@";
  167         $message .= '. Did you forget to supply the second (imaginary) ndarray?'
  168             if ($message =~ /undefined value/);
  169         barf($message);
  170     }
  171     PDL::_fft_int($re,$im);
  172     if (!$_[0]->type->real) {
  173         $_[0]= czip($re, $im);
  174     } else {
  175         $_[0]=$re,$_[1]=$im;
  176     }
  177 }
  178 #line 179 "FFT.pm"
  179 
  180 
  181 
  182 #line 950 "../../blib/lib/PDL/PP.pm"
  183 
  184 *fft = \&PDL::fft;
  185 #line 186 "FFT.pm"
  186 
  187 
  188 
  189 #line 948 "../../blib/lib/PDL/PP.pm"
  190 
  191 
  192 
  193 =head2 ifft
  194 
  195 =for sig
  196 
  197   Signature: ([io]real(n); [io]imag(n))
  198 
  199 =for ref
  200 
  201 Complex inverse 1-D FFT of the "real" and "imag" arrays [inplace]. A single
  202 cfloat/cdouble input ndarray can also be used.
  203 
  204 =for usage
  205 
  206   ifft($real,$imag);
  207   ifft($complex);
  208 
  209 
  210 =for bad
  211 
  212 ifft does not process bad values.
  213 It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays.
  214 
  215 
  216 =cut
  217 #line 218 "FFT.pm"
  218 
  219 
  220 
  221 #line 949 "../../blib/lib/PDL/PP.pm"
  222 
  223 sub PDL::ifft {
  224     # Convert the first argument to decimal and check for trouble.
  225     my ($re, $im) = @_;
  226     if (!$re->type->real) {
  227         $im=$re->im;
  228         $re=$re->re;
  229     }
  230     eval {  todecimal($re); };
  231     if ($@) {
  232         $@ =~ s/ at .*//s;
  233         barf("Error in FFT with first argument: $@");
  234     }
  235     # Convert the second argument to decimal and check for trouble.
  236     eval {  todecimal($im); };
  237     if ($@) {
  238         $@ =~ s/ at .*//s;
  239         my $message = "Error in FFT with second argument: $@";
  240         $message .= '. Did you forget to supply the second (imaginary) ndarray?'
  241             if ($message =~ /undefined value/);
  242         barf($message);
  243     }
  244     PDL::_ifft_int($re,$im);
  245     if (!$_[0]->type->real) {
  246         $_[0]= czip($re, $im);
  247     } else {
  248         $_[0]=$re,$_[1]=$im;
  249     }
  250 }
  251 #line 252 "FFT.pm"
  252 
  253 
  254 
  255 #line 950 "../../blib/lib/PDL/PP.pm"
  256 
  257 *ifft = \&PDL::ifft;
  258 #line 259 "FFT.pm"
  259 
  260 
  261 
  262 #line 185 "fft.pd"
  263 
  264 use Carp;
  265 use PDL::Core qw/:Func/;
  266 use PDL::Basic qw/:Func/;
  267 use PDL::Types;
  268 use PDL::ImageND qw/kernctr/; # moved to ImageND since FFTW uses it too
  269 use PDL::Ops qw/czip/;
  270 
  271 sub todecimal {
  272     my ($arg) = @_;
  273     $arg = $arg->double if $arg->type->integer;
  274     $_[0] = $arg;
  275 1;}
  276 
  277 =head2 realfft()
  278 
  279 =for ref
  280 
  281 One-dimensional FFT of real function [inplace].
  282 
  283 The real part of the transform ends up in the first half of the array
  284 and the imaginary part of the transform ends up in the second half of
  285 the array.
  286 
  287 =for usage
  288 
  289     realfft($real);
  290 
  291 =cut
  292 
  293 *realfft = \&PDL::realfft;
  294 
  295 sub PDL::realfft {
  296     barf("Usage: realfft(real(*)") if $#_ != 0;
  297     my ($x) = @_;
  298     todecimal($x);
  299 # FIX: could eliminate $y
  300     my ($y) = 0*$x;
  301     fft($x,$y);
  302     my ($n) = int((($x->dims)[0]-1)/2); my($t);
  303     ($t=$x->slice("-$n:-1")) .= $y->slice("1:$n");
  304     undef;
  305 }
  306 
  307 =head2 realifft()
  308 
  309 =for ref
  310 
  311 Inverse of one-dimensional realfft routine [inplace].
  312 
  313 =for usage
  314 
  315     realifft($real);
  316 
  317 =cut
  318 
  319 *realifft = \&PDL::realifft;
  320 
  321 sub PDL::realifft {
  322     use PDL::Ufunc 'max';
  323     barf("Usage: realifft(xfm(*)") if $#_ != 0;
  324     my ($x) = @_;
  325     todecimal($x);
  326     my ($n) = int((($x->dims)[0]-1)/2); my($t);
  327 # FIX: could eliminate $y
  328     my ($y) = 0*$x;
  329     ($t=$y->slice("1:$n")) .= $x->slice("-$n:-1");
  330     ($t=$x->slice("-$n:-1")) .= $x->slice("$n:1");
  331     ($t=$y->slice("-$n:-1")) .= -$y->slice("$n:1");
  332     ifft($x,$y);
  333 # Sanity check -- shouldn't happen
  334     carp "Bad inverse transform in realifft" if max(abs($y)) > 1e-6*max(abs($x));
  335     undef;
  336 }
  337 
  338 =head2 fftnd()
  339 
  340 =for ref
  341 
  342 N-dimensional FFT over all pdl dims of input (inplace) 
  343 
  344 =for example
  345 
  346     fftnd($real,$imag);
  347 
  348 =cut
  349 
  350 *fftnd = \&PDL::fftnd;
  351 
  352 sub PDL::fftnd {
  353     my ($r,$i) = @_;
  354     barf "Must have real and imaginary parts or complex for fftnd"
  355       if $r->type->real and @_ != 2;
  356     if (!$r->type->real) {
  357     $i=$r->im;
  358     $r=$r->re;
  359     }
  360     my ($n) = $r->getndims;
  361     barf "Dimensions of real and imag must be the same for fft"
  362         if ($n != $i->getndims);
  363     $n--;
  364     todecimal($r);
  365     todecimal($i);
  366     # need the copy in case $r and $i point to same memory
  367     $i = $i->copy;
  368     foreach (0..$n) {
  369       fft($r,$i);
  370       $r = $r->mv(0,$n);
  371       $i = $i->mv(0,$n);
  372     }
  373     if (!$_[0]->type->real) {
  374     $_[0]= czip($r, $i);
  375     } else {
  376     $_[0] = $r; $_[1] = $i;
  377     }
  378     undef;
  379 }
  380 
  381 =head2 ifftnd()
  382 
  383 =for ref
  384 
  385 N-dimensional inverse FFT over all pdl dims of input (inplace) 
  386 
  387 =for example
  388 
  389     ifftnd($real,$imag);
  390 
  391 =cut
  392 
  393 *ifftnd = \&PDL::ifftnd;
  394 
  395 sub PDL::ifftnd {
  396     my ($r,$i) = @_;
  397     barf "Must have real and imaginary parts or complex for ifftnd"
  398       if $r->type->real and @_ != 2;
  399     if (!$r->type->real) {
  400     $i=$r->im;
  401     $r=$r->re;
  402     }
  403     my ($n) = $r->getndims;
  404     barf "Dimensions of real and imag must be the same for ifft"
  405         if ($n != $i->getndims);
  406     todecimal($r);
  407     todecimal($i);
  408     # need the copy in case $r and $i point to same memory
  409     $i = $i->copy;
  410     $n--;
  411     foreach (0..$n) {
  412       ifft($r,$i);
  413       $r = $r->mv(0,$n);
  414       $i = $i->mv(0,$n);
  415     }
  416     if (!$_[0]->type->real) {
  417     $_[0]= czip($r, $i);
  418     } else {
  419     $_[0] = $r; $_[1] = $i;
  420     }
  421     undef;
  422 }
  423 
  424 =head2 fftconvolve()
  425 
  426 =for ref
  427 
  428 N-dimensional convolution with periodic boundaries (FFT method)
  429 
  430 =for usage
  431 
  432     $kernel = kernctr($image,$smallk);
  433     fftconvolve($image,$kernel);
  434 
  435 fftconvolve works inplace, and returns an error array in kernel as an
  436 accuracy check -- all the values in it should be negligible.
  437 
  438 See also L<PDL::ImageND::convolveND|PDL::ImageND/convolveND>, which 
  439 performs speed-optimized convolution with a variety of boundary conditions.
  440 
  441 The sizes of the image and the kernel must be the same.
  442 L<kernctr|PDL::ImageND/kernctr> centres a small kernel to emulate the
  443 behaviour of the direct convolution routines.
  444 
  445 The speed cross-over between using straight convolution 
  446 (L<PDL::Image2D::conv2d()|PDL::Image2D/conv2d>) and
  447 these fft routines is for kernel sizes roughly 7x7.
  448 
  449 =cut
  450 
  451 *fftconvolve = \&PDL::fftconvolve;
  452 
  453 sub PDL::fftconvolve {
  454     barf "Must have image & kernel for fftconvolve" if $#_ != 1;
  455     my ($im, $k) = map $_->r2C, @_;
  456     fftnd($im);
  457     fftnd($k);
  458     my $c = $im * $k;
  459     ifftnd($c);
  460     $_[0] = $c->re->sever;
  461     $_[1] = $c->im->sever;
  462     @_;
  463 }
  464 #line 465 "FFT.pm"
  465 
  466 
  467 
  468 
  469 
  470 #line 388 "fft.pd"
  471 
  472 =head1 BUGS
  473 
  474 Where the source is marked `FIX', could re-implement using phase-shift
  475 factors on the transforms and some real-space bookkeeping, to save
  476 some temporary space and redundant transforms.
  477 
  478 =head1 AUTHOR
  479 
  480 This file copyright (C) 1997, 1998 R.J.R. Williams
  481 (rjrw@ast.leeds.ac.uk), Karl Glazebrook (kgb@aaoepp.aao.gov.au),
  482 Tuomas J. Lukka, (lukka@husc.harvard.edu).  All rights reserved. There
  483 is no warranty. You are allowed to redistribute this software /
  484 documentation under certain conditions. For details, see the file
  485 COPYING in the PDL distribution. If this file is separated from the
  486 PDL distribution, the copyright notice should be included in the file.
  487 
  488 =cut
  489 #line 490 "FFT.pm"
  490 
  491 
  492 
  493 
  494 # Exit with OK status
  495 
  496 1;