# # GENERATED WITH PDL::PP! Don't modify! # package PDL::FFT; our @EXPORT_OK = qw(fft ifft fftnd ifftnd fftconvolve realfft realifft kernctr ); our %EXPORT_TAGS = (Func=>\@EXPORT_OK); use PDL::Core; use PDL::Exporter; use DynaLoader; our @ISA = ( 'PDL::Exporter','DynaLoader' ); push @PDL::Core::PP, __PACKAGE__; bootstrap PDL::FFT ; #line 7 "fft.pd" =head1 NAME PDL::FFT - FFTs for PDL =head1 DESCRIPTION !!!!!!!!!!!!!!!!!!!!!!!!!!WARNING!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! As of PDL-2.006_04, the direction of the FFT/IFFT has been reversed to match the usage in the FFTW library and the convention in use generally. !!!!!!!!!!!!!!!!!!!!!!!!!!WARNING!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! FFTs for PDL. These work for arrays of any dimension, although ones with small prime factors are likely to be the quickest. The forward FFT is unnormalized while the inverse FFT is normalized so that the IFFT of the FFT returns the original values. For historical reasons, these routines work in-place and do not recognize the in-place flag. That should be fixed. =head1 SYNOPSIS use PDL::FFT qw/:Func/; fft($real, $imag); ifft($real, $imag); realfft($real); realifft($real); fftnd($real,$imag); ifftnd($real,$imag); $kernel = kernctr($image,$smallk); fftconvolve($image,$kernel); =head1 DATA TYPES The underlying C library upon which this module is based performs FFTs on both single precision and double precision floating point ndarrays. The PP functions are defined to only take those data types. Therefore, if you pass in an ndarray of integer datatype (byte, short, ushort, long) to any of the routines in PDL::FFT, your data will be promoted to a double-precision ndarray. If you pass in a float, the single-precision FFT will be performed. =head1 FREQUENCIES For even-sized input arrays, the frequencies are packed like normal for FFTs (where N is the size of the array and D is the physical step size between elements): 0, 1/ND, 2/ND, ..., (N/2-1)/ND, 1/2D, -(N/2-1)/ND, ..., -1/ND. which can easily be obtained (taking the Nyquist frequency to be positive) using C<< $kx = $real->xlinvals(-($N/2-1)/$N/$D,1/2/$D)->rotate(-($N/2 -1)); >> For odd-sized input arrays the Nyquist frequency is not directly acessible, and the frequencies are 0, 1/ND, 2/ND, ..., (N/2-0.5)/ND, -(N/2-0.5)/ND, ..., -1/ND. which can easily be obtained using C<< $kx = $real->xlinvals(-($N/2-0.5)/$N/$D,($N/2-0.5)/$N/$D)->rotate(-($N-1)/2); >> =head1 ALTERNATIVE FFT PACKAGES Various other modules - such as L<PDL::FFTW3> and L<PDL::Slatec> - contain FFT routines. However, unlike PDL::FFT, these modules are optional, and so may not be installed. =cut #line 102 "FFT.pm" =head1 FUNCTIONS =cut #line 1059 "/home/osboxes/.perlbrew/libs/perl-5.32.0@normal/lib/perl5/x86_64-linux/PDL/PP.pm" =head2 fft =for sig Signature: ([io]real(n); [io]imag(n)) =for ref Complex 1-D FFT of the "real" and "imag" arrays [inplace]. A single cfloat/cdouble input ndarray can also be used. =for usage fft($real,$imag); fft($complex); =for bad fft does not process bad values. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. =cut #line 143 "FFT.pm" #line 1060 "/home/osboxes/.perlbrew/libs/perl-5.32.0@normal/lib/perl5/x86_64-linux/PDL/PP.pm" sub PDL::fft { # Convert the first argument to decimal and check for trouble. my ($re, $im) = @_; if (!$re->type->real) { $im=$re->im; $re=$re->re; } eval { todecimal($re); }; if ($@) { $@ =~ s/ at .*//s; barf("Error in FFT with first argument: $@"); } # Convert the second argument to decimal and check for trouble. eval { todecimal($im); }; if ($@) { $@ =~ s/ at .*//s; my $message = "Error in FFT with second argument: $@"; $message .= '. Did you forget to supply the second (imaginary) ndarray?' if ($message =~ /undefined value/); barf($message); } PDL::_fft_int($re,$im); if (!$_[0]->type->real) { $_[0]= czip($re, $im); } else { $_[0]=$re,$_[1]=$im; } } #line 176 "FFT.pm" #line 1061 "/home/osboxes/.perlbrew/libs/perl-5.32.0@normal/lib/perl5/x86_64-linux/PDL/PP.pm" *fft = \&PDL::fft; #line 182 "FFT.pm" #line 1059 "/home/osboxes/.perlbrew/libs/perl-5.32.0@normal/lib/perl5/x86_64-linux/PDL/PP.pm" =head2 ifft =for sig Signature: ([io]real(n); [io]imag(n)) =for ref Complex inverse 1-D FFT of the "real" and "imag" arrays [inplace]. A single cfloat/cdouble input ndarray can also be used. =for usage ifft($real,$imag); ifft($complex); =for bad ifft does not process bad values. It will set the bad-value flag of all output ndarrays if the flag is set for any of the input ndarrays. =cut #line 213 "FFT.pm" #line 1060 "/home/osboxes/.perlbrew/libs/perl-5.32.0@normal/lib/perl5/x86_64-linux/PDL/PP.pm" sub PDL::ifft { # Convert the first argument to decimal and check for trouble. my ($re, $im) = @_; if (!$re->type->real) { $im=$re->im; $re=$re->re; } eval { todecimal($re); }; if ($@) { $@ =~ s/ at .*//s; barf("Error in FFT with first argument: $@"); } # Convert the second argument to decimal and check for trouble. eval { todecimal($im); }; if ($@) { $@ =~ s/ at .*//s; my $message = "Error in FFT with second argument: $@"; $message .= '. Did you forget to supply the second (imaginary) ndarray?' if ($message =~ /undefined value/); barf($message); } PDL::_ifft_int($re,$im); if (!$_[0]->type->real) { $_[0]= czip($re, $im); } else { $_[0]=$re,$_[1]=$im; } } #line 246 "FFT.pm" #line 1061 "/home/osboxes/.perlbrew/libs/perl-5.32.0@normal/lib/perl5/x86_64-linux/PDL/PP.pm" *ifft = \&PDL::ifft; #line 252 "FFT.pm" #line 186 "fft.pd" use Carp; use PDL::Core qw/:Func/; use PDL::Basic qw/:Func/; use PDL::Types; use PDL::ImageND qw/kernctr/; # moved to ImageND since FFTW uses it too use PDL::Ops qw/czip/; sub todecimal { my ($arg) = @_; $arg = $arg->double if $arg->type->integer; $_[0] = $arg; 1;} =head2 realfft() =for ref One-dimensional FFT of real function [inplace]. The real part of the transform ends up in the first half of the array and the imaginary part of the transform ends up in the second half of the array. =for usage realfft($real); =cut *realfft = \&PDL::realfft; sub PDL::realfft { barf("Usage: realfft(real(*)") if $#_ != 0; my ($x) = @_; todecimal($x); # FIX: could eliminate $y my ($y) = 0*$x; fft($x,$y); my ($n) = int((($x->dims)[0]-1)/2); my($t); ($t=$x->slice("-$n:-1")) .= $y->slice("1:$n"); undef; } =head2 realifft() =for ref Inverse of one-dimensional realfft routine [inplace]. =for usage realifft($real); =cut *realifft = \&PDL::realifft; sub PDL::realifft { use PDL::Ufunc 'max'; barf("Usage: realifft(xfm(*)") if $#_ != 0; my ($x) = @_; todecimal($x); my ($n) = int((($x->dims)[0]-1)/2); my($t); # FIX: could eliminate $y my ($y) = 0*$x; ($t=$y->slice("1:$n")) .= $x->slice("-$n:-1"); ($t=$x->slice("-$n:-1")) .= $x->slice("$n:1"); ($t=$y->slice("-$n:-1")) .= -$y->slice("$n:1"); ifft($x,$y); # Sanity check -- shouldn't happen carp "Bad inverse transform in realifft" if max(abs($y)) > 1e-6*max(abs($x)); undef; } =head2 fftnd() =for ref N-dimensional FFT over all pdl dims of input (inplace) =for example fftnd($real,$imag); =cut *fftnd = \&PDL::fftnd; sub PDL::fftnd { my ($r,$i) = @_; barf "Must have real and imaginary parts or complex for fftnd" if $r->type->real and @_ != 2; if (!$r->type->real) { $i=$r->im; $r=$r->re; } my ($n) = $r->getndims; barf "Dimensions of real and imag must be the same for fft" if ($n != $i->getndims); $n--; todecimal($r); todecimal($i); # need the copy in case $r and $i point to same memory $i = $i->copy; foreach (0..$n) { fft($r,$i); $r = $r->mv(0,$n); $i = $i->mv(0,$n); } if (!$_[0]->type->real) { $_[0]= czip($r, $i); } else { $_[0] = $r; $_[1] = $i; } undef; } =head2 ifftnd() =for ref N-dimensional inverse FFT over all pdl dims of input (inplace) =for example ifftnd($real,$imag); =cut *ifftnd = \&PDL::ifftnd; sub PDL::ifftnd { my ($r,$i) = @_; barf "Must have real and imaginary parts or complex for ifftnd" if $r->type->real and @_ != 2; if (!$r->type->real) { $i=$r->im; $r=$r->re; } my ($n) = $r->getndims; barf "Dimensions of real and imag must be the same for ifft" if ($n != $i->getndims); todecimal($r); todecimal($i); # need the copy in case $r and $i point to same memory $i = $i->copy; $n--; foreach (0..$n) { ifft($r,$i); $r = $r->mv(0,$n); $i = $i->mv(0,$n); } if (!$_[0]->type->real) { $_[0]= czip($r, $i); } else { $_[0] = $r; $_[1] = $i; } undef; } =head2 fftconvolve() =for ref N-dimensional convolution with periodic boundaries (FFT method) =for usage $kernel = kernctr($image,$smallk); fftconvolve($image,$kernel); fftconvolve works inplace, and returns an error array in kernel as an accuracy check -- all the values in it should be negligible. See also L<PDL::ImageND::convolveND|PDL::ImageND/convolveND>, which performs speed-optimized convolution with a variety of boundary conditions. The sizes of the image and the kernel must be the same. L<kernctr|PDL::ImageND/kernctr> centres a small kernel to emulate the behaviour of the direct convolution routines. The speed cross-over between using straight convolution (L<PDL::Image2D::conv2d()|PDL::Image2D/conv2d>) and these fft routines is for kernel sizes roughly 7x7. =cut *fftconvolve = \&PDL::fftconvolve; sub PDL::fftconvolve { barf "Must have image & kernel for fftconvolve" if $#_ != 1; my ($im, $k) = map $_->r2C, @_; fftnd($im); fftnd($k); my $c = $im * $k; ifftnd($c); $_[0] = $c->re->sever; $_[1] = $c->im->sever; @_; } #line 457 "FFT.pm" #line 389 "fft.pd" =head1 BUGS Where the source is marked `FIX', could re-implement using phase-shift factors on the transforms and some real-space bookkeeping, to save some temporary space and redundant transforms. =head1 AUTHOR This file copyright (C) 1997, 1998 R.J.R. Williams (rjrw@ast.leeds.ac.uk), Karl Glazebrook (kgb@aaoepp.aao.gov.au), Tuomas J. Lukka, (lukka@husc.harvard.edu). 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. =cut #line 481 "FFT.pm" # Exit with OK status 1;