~ubuntu-branches/ubuntu/maverick/pdl/maverick

« back to all changes in this revision

Viewing changes to Lib/FFT/fft.pd

• browse files at revision 16
• compare with another revision
• download diff
• download tarball
• view history from revision 16

Show diffs side-by-side

added

removed

Lib/FFT/fft.pd

27

27

        $kernel = kernctr($image,$smallk);

28

28

        fftconvolve($image,$kernel);

29

29

 

 

30

=head1 DATA TYPES

 

31

 

 

32

The underlying C library upon which this module is based performs FFTs

 

33

on both single precision and double precision floating point piddles.

 

34

Performing FFTs on integer data types is not reliable.  Consider the

 

35

following FFT on piddles of type 'double':

 

36

 

 

37

        $r = pdl(0,1,0,1);

 

38

        $i = zeroes($r);

 

39

        fft($r,$i);

 

40

        print $r,$i;

 

41

        [2 0 -2 0] [0 0 0 0]

 

42

 

 

43

But if $r and $i are unsigned short integers (ushorts):

 

44

 

 

45

        $r = pdl(ushort,0,1,0,1);

 

46

        $i = zeroes($r);

 

47

        fft($r,$i);

 

48

        print $r,$i;

 

49

        [2 0 65534 0] [0 0 0 0]

 

50

 

 

51

This used to occur because L<PDL::PP|PDL::PP> converts the ushort

 

52

piddles to floats or doubles, performs the FFT on them, and then

 

53

converts them back to ushort, causing the overflow where the amplitude

 

54

of the frequency should be -2.

 

55

 

 

56

Therefore, if you pass in a piddle of integer datatype (byte, short,

 

57

ushort, long) to any of the routines in PDL::FFT, your data will be

 

58

promoted to a double-precision piddle.  If you pass in a float, the

 

59

single-precision FFT will be performed.

 

60

 

 

61

=head1 FREQUENCIES

 

62

 

 

63

For even-sized input arrays, the frequencies are packed like normal

 

64

for FFTs (where N is the size of the array and D is the physical step

 

65

size between elements):

 

66

 

 

67

0, 1/ND, 2/ND, ..., (N/2-1)/ND, 1/2D, -(N/2-1)/ND, ..., -1/ND.

 

68

 

 

69

which can easily be obtained (taking the Nyquist frequency to be

 

70

positive) using

 

71

 

 

72

$kx = $real->xlinvals(-($N/2-1)/$N/$D,1/2/$D)->rotate(-($N/2 -1));

 

73

 

 

74

For odd-sized input arrays the Nyquist frequency is not directly

 

75

acessible, and the frequencies are

 

76

 

 

77

0, 1/ND, 2/ND, ..., (N/2-0.5)/ND, -(N/2-0.5)/ND, ..., -1/ND.

 

78

 

 

79

which can easily be obtained using

 

80

 

 

81

$kx = $real->xlinvals(-($N/2-0.5)/$N/$D,($N/2-0.5)/$N/$D)->rotate(-($N-1)/2);

 

82

 

 

83

 

30

84

=head1 ALTERNATIVE FFT PACKAGES

31

85

 

32

86

Various other modules - such as 

56

110

   OUTPUT:

57

111

     RETVAL

58

112

');

59

 

pp_def('fft',

60

 

        Pars => '[o,nc]real(n); [o,nc]imag(n);',

61

 

        GenericTypes => [F,D],

62

 

        Code => '$TFD(fftnf,fftn)

63

 

          ($SIZE(n), NULL , $P(real),$P(imag), 1, 1.);',

64

 

        Doc=>'Complex FFT of the "real" and "imag" arrays [inplace]'

65

 

);

66

 

 

67

 

pp_def('ifft',

68

 

        Pars => '[o,nc]real(n); [o,nc]imag(n);',

69

 

        GenericTypes => [F,D],

70

 

        Code => '$TFD(fftnf,fftn)

71

 

          ($SIZE(n), NULL , $P(real),$P(imag), -1, -1.);',

72

 

        Doc=>'Complex Inverse FFT of the "real" and "imag" arrays [inplace]'

73

 

);

74

 

 

75

 

pp_add_exported('',"fftnd ifftnd fftconvolve realfft realifft kernctr");

 

113

pp_def('_fft',

 

114

        Pars => '[o,nc]real(n); [o,nc]imag(n);',

 

115

        GenericTypes => [F,D],

 

116

        Code => '$TFD(fftnf,fftn)

 

117

        ($SIZE(n), NULL , $P(real),$P(imag), 1, 1.);',

 

118

        Doc=>undef

 

119

);

 

120

 

 

121

pp_def('_ifft',

 

122

        Pars => '[o,nc]real(n); [o,nc]imag(n);',

 

123

        GenericTypes => [F,D],

 

124

        Code => '$TFD(fftnf,fftn)

 

125

        ($SIZE(n), NULL , $P(real),$P(imag), -1, -1.);',

 

126

        Doc=>undef

 

127

);

 

128

 

 

129

pp_add_exported('',"fft ifft fftnd ifftnd fftconvolve realfft realifft kernctr");

76

130

 

77

131

pp_addpm(<<'EOD');

78

132

 

88

142

  fft_free();

89

143

}

90

144

 

91

 

sub todouble {

 

145

sub todecimal {

92

146

    my ($arg) = @_;

93

 

    $arg = $arg->double if $arg->get_datatype != $PDL_D;

 

147

    $arg = $arg->double if (($arg->get_datatype != $PDL_F) && 

 

148

                           ($arg->get_datatype != $PDL_D));

94

149

    $_[0] = $arg;

95

150

1;}

96

151

 

 

152

=head2 fft()

 

153

 

 

154

=for ref

 

155

 

 

156

Complex FFT of the "real" and "imag" arrays [inplace].

 

157

 

 

158

=for usage

 

159

 

 

160

fft($real,$imag);

 

161

 

 

162

=cut

 

163

 

 

164

*fft = \&PDL::fft;

 

165

 

 

166

sub PDL::fft {

 

167

todecimal($_[0]);

 

168

todecimal($_[1]);

 

169

_fft($_[0],$_[1]);

 

170

}

 

171

 

 

172

 

 

173

=head2 ifft()

 

174

 

 

175

=for ref

 

176

 

 

177

Complex inverse FFT of the "real" and "imag" arrays [inplace].

 

178

 

 

179

=for usage

 

180

 

 

181

ifft($real,$imag);

 

182

 

 

183

=cut

 

184

 

 

185

*ifft = \&PDL::ifft;

 

186

 

 

187

sub PDL::ifft {

 

188

todecimal($_[0]);

 

189

todecimal($_[1]);

 

190

_ifft($_[0],$_[1]);

 

191

}

 

192

 

97

193

=head2 realfft()

98

194

 

99

195

=for ref

115

211

sub PDL::realfft {

116

212

    barf("Usage: realfft(real(*)") if $#_ != 0;

117

213

    my ($a) = @_;

118

 

    todouble($a);

 

214

    todecimal($a);

119

215

# FIX: could eliminate $b

120

216

    my ($b) = 0*$a;

121

217

    fft($a,$b);

142

238

    use PDL::Ufunc 'max';

143

239

    barf("Usage: realifft(xfm(*)") if $#_ != 0;

144

240

    my ($a) = @_;

145

 

    todouble($a);

 

241

    todecimal($a);

146

242

    my ($n) = int((($a->dims)[0]-1)/2); my($t);

147

243

# FIX: could eliminate $b

148

244

    my ($b) = 0*$a;

176

272

    barf "Dimensions of real and imag must be the same for fft"

177

273

        if ($n != $i->getndims);

178

274

    $n--;

179

 

    todouble($r);

180

 

    todouble($i);

 

275

    todecimal($r);

 

276

    todecimal($i);

181

277

    # need the copy in case $r and $i point to same memory

182

278

    $i = $i->copy;

183

279

    foreach (0..$n) {

209

305

    my ($n) = $r->getndims;

210

306

    barf "Dimensions of real and imag must be the same for ifft"

211

307

        if ($n != $i->getndims);

212

 

    todouble($r);

213

 

    todouble($i);

 

308

    todecimal($r);

 

309

    todecimal($i);

214

310

    # need the copy in case $r and $i point to same memory

215

311

    $i = $i->copy;

216

312

    $n--;

327

423

    barf "Must have image & kernel for fftconvolve" if $#_ != 1;

328

424

    my ($hr, $hi) = @_;

329

425

    my ($n) = $hr->getndims;

330

 

    todouble($hr);

331

 

    todouble($hi);

 

426

    todecimal($hr);

 

427

    todecimal($hi);

332

428

    # need the copy in case $r and $i point to same memory

333

429

    $hi = $hi->copy;

334

430

    $hr = $hr->copy;

355

451

    ($hr,$hi);

356

452

}

357

453

 

358

 

 

359

 

=cut

360

 

 

361

454

EOD

362

455

 

363

456

# convmath does local part of the maths necessary to handle a,b which

• Older »

Loggerhead is a web-based interface for Breezy
Version: 2.0.1