1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
|
%perlcode %{
@EXPORT_complex = qw/
gsl_fft_complex_radix2_forward
gsl_fft_complex_radix2_backward
gsl_fft_complex_radix2_inverse
gsl_fft_complex_radix2_transform
gsl_fft_complex_radix2_dif_forward
gsl_fft_complex_radix2_dif_backward
gsl_fft_complex_radix2_dif_inverse
gsl_fft_complex_radix2_dif_transform
gsl_fft_complex_wavetable_alloc
gsl_fft_complex_wavetable_free
gsl_fft_complex_workspace_alloc
gsl_fft_complex_workspace_free
gsl_fft_complex_memcpy
gsl_fft_complex_forward
gsl_fft_complex_backward
gsl_fft_complex_inverse
gsl_fft_complex_transform
/;
@EXPORT_halfcomplex = qw/
gsl_fft_halfcomplex_radix2_backward
gsl_fft_halfcomplex_radix2_inverse
gsl_fft_halfcomplex_radix2_transform
gsl_fft_halfcomplex_wavetable_alloc
gsl_fft_halfcomplex_wavetable_free
gsl_fft_halfcomplex_backward
gsl_fft_halfcomplex_inverse
gsl_fft_halfcomplex_transform
gsl_fft_halfcomplex_unpack
gsl_fft_halfcomplex_radix2_unpack
/;
@EXPORT_real = qw/
gsl_fft_real_radix2_transform
gsl_fft_real_wavetable_alloc
gsl_fft_real_wavetable_free
gsl_fft_real_workspace_alloc
gsl_fft_real_workspace_free
gsl_fft_real_transform
gsl_fft_real_unpack
/;
@EXPORT_vars = qw/
$gsl_fft_forward
$gsl_fft_backward
/;
@EXPORT_OK = (
@EXPORT_real,
@EXPORT_complex,
@EXPORT_halfcomplex,
@EXPORT_vars,
);
%EXPORT_TAGS = (
all => \@EXPORT_OK,
real => \@EXPORT_real,
complex => \@EXPORT_complex,
halfcomplex => \@EXPORT_halfcomplex,
vars => \@EXPORT_vars,
);
__END__
=encoding utf8
=head1 NAME
Math::GSL::FFT - Fast Fourier Transforms (FFT)
=head1 SYNOPSIS
use Math::GSL::FFT qw /:all/;
my $input1 = [ 0 .. 7 ];
my $N1 = @$input1;
my ($status1, $output1) = gsl_fft_real_radix2_transform ($input, 1, $N1);
my ($status2, $output2) = gsl_fft_halfcomplex_radix2_inverse($output2, 1, $N1);
# $input1 == $output2
my $input2 = [ 0 .. 6 ];
my $N2 = @$input;
my $workspace1 = gsl_fft_real_workspace_alloc($N2);
my $wavetable1 = gsl_fft_real_wavetable_alloc($N2);
my ($status3,$output3) = gsl_fft_real_transform ($input, 1, $N2, $wavetable1, $workspace1);
my $wavetable4 = gsl_fft_halfcomplex_wavetable_alloc($N2);
my $workspace4 = gsl_fft_real_workspace_alloc($N2);
my ($status4,$output4) = gsl_fft_halfcomplex_inverse($output, 1, $N2, $wavetable4, $workspace4);
# $input2 == $output4
=head1 DESCRIPTION
=over
=item * C<gsl_fft_complex_radix2_forward($data, $stride, $n) >
This function computes the forward FFTs of length $n with stride $stride, on
the array reference $data using an in-place radix-2 decimation-in-time
algorithm. The length of the transform $n is restricted to powers of two. For
the transform version of the function the sign argument can be either forward
(-1) or backward (+1). The functions return a value of $GSL_SUCCESS if no
errors were detected, or $GSL_EDOM if the length of the data $n is not a power
of two.
=item * C<gsl_fft_complex_radix2_backward >
=item * C<gsl_fft_complex_radix2_inverse >
=item * C<gsl_fft_complex_radix2_transform >
=item * C<gsl_fft_complex_radix2_dif_forward >
=item * C<gsl_fft_complex_radix2_dif_backward >
=item * C<gsl_fft_complex_radix2_dif_inverse >
=item * C<gsl_fft_complex_radix2_dif_transform >
=item * C<gsl_fft_complex_wavetable_alloc($n)>
This function prepares a trigonometric lookup table for a complex FFT of length
$n. The function returns a pointer to the newly allocated
gsl_fft_complex_wavetable if no errors were detected, and a null pointer in the
case of error. The length $n is factorized into a product of subtransforms, and
the factors and their trigonometric coefficients are stored in the wavetable.
The trigonometric coefficients are computed using direct calls to sin and cos,
for accuracy. Recursion relations could be used to compute the lookup table
faster, but if an application performs many FFTs of the same length then this
computation is a one-off overhead which does not affect the final throughput.
The wavetable structure can be used repeatedly for any transform of the same
length. The table is not modified by calls to any of the other FFT functions.
The same wavetable can be used for both forward and backward (or inverse)
transforms of a given length.
=item * C<gsl_fft_complex_wavetable_free($wavetable)>
This function frees the memory associated with the wavetable $wavetable. The
wavetable can be freed if no further FFTs of the same length will be needed.
=item * C<gsl_fft_complex_workspace_alloc($n)>
This function allocates a workspace for a complex transform of length $n.
=item * C<gsl_fft_complex_workspace_free($workspace) >
This function frees the memory associated with the workspace $workspace. The
workspace can be freed if no further FFTs of the same length will be needed.
=item * C<gsl_fft_complex_memcpy >
=item * C<gsl_fft_complex_forward >
=item * C<gsl_fft_complex_backward >
=item * C<gsl_fft_complex_inverse >
=item * C<gsl_fft_complex_transform >
=item * C<gsl_fft_halfcomplex_radix2_backward($data, $stride, $n)>
This function computes the backwards in-place radix-2 FFT of length $n and
stride $stride on the half-complex sequence data stored according the output
scheme used by gsl_fft_real_radix2. The result is a real array stored in
natural order.
=item * C<gsl_fft_halfcomplex_radix2_inverse($data, $stride, $n)>
This function computes the inverse in-place radix-2 FFT of length $n and stride
$stride on the half-complex sequence data stored according the output scheme
used by gsl_fft_real_radix2. The result is a real array stored in natural
order.
=item * C<gsl_fft_halfcomplex_radix2_transform>
=item * C<gsl_fft_halfcomplex_wavetable_alloc($n)>
This function prepares trigonometric lookup tables for an FFT of size $n real
elements. The functions return a pointer to the newly allocated struct if no
errors were detected, and a null pointer in the case of error. The length $n is
factorized into a product of subtransforms, and the factors and their
trigonometric coefficients are stored in the wavetable. The trigonometric
coefficients are computed using direct calls to sin and cos, for accuracy.
Recursion relations could be used to compute the lookup table faster, but if an
application performs many FFTs of the same length then computing the wavetable
is a one-off overhead which does not affect the final throughput. The
wavetable structure can be used repeatedly for any transform of the same
length. The table is not modified by calls to any of the other FFT functions.
The appropriate type of wavetable must be used for forward real or inverse
half-complex transforms.
=item * C<gsl_fft_halfcomplex_wavetable_free($wavetable)>
This function frees the memory associated with the wavetable $wavetable. The
wavetable can be freed if no further FFTs of the same length will be needed.
=item * C<gsl_fft_halfcomplex_backward >
=item * C<gsl_fft_halfcomplex_inverse >
=item * C<gsl_fft_halfcomplex_transform >
=item * C<gsl_fft_halfcomplex_unpack >
=item * C<gsl_fft_halfcomplex_radix2_unpack >
=item * C<gsl_fft_real_radix2_transform($data, $stride, $n) >
This function computes an in-place radix-2 FFT of length $n and stride $stride
on the real array reference $data. The output is a half-complex sequence, which
is stored in-place. The arrangement of the half-complex terms uses the
following scheme: for k < N/2 the real part of the k-th term is stored in
location k, and the corresponding imaginary part is stored in location N-k.
Terms with k > N/2 can be reconstructed using the symmetry z_k = z^*_{N-k}. The
terms for k=0 and k=N/2 are both purely real, and count as a special case.
Their real parts are stored in locations 0 and N/2 respectively, while their
imaginary parts which are zero are not stored. The following table shows the
correspondence between the output data and the equivalent results obtained by
considering the input data as a complex sequence with zero imaginary part,
complex[0].real = data[0]
complex[0].imag = 0
complex[1].real = data[1]
complex[1].imag = data[N-1]
............... ................
complex[k].real = data[k]
complex[k].imag = data[N-k]
............... ................
complex[N/2].real = data[N/2]
complex[N/2].imag = 0
............... ................
complex[k'].real = data[k] k' = N - k
complex[k'].imag = -data[N-k]
............... ................
complex[N-1].real = data[1]
complex[N-1].imag = -data[N-1]
=for notyou #'
Note that the output data can be converted into the full complex sequence using
the function gsl_fft_halfcomplex_unpack.
=item * C<gsl_fft_real_wavetable_alloc($n)>
This function prepares trigonometric lookup tables for an FFT of size $n real
elements. The functions return a pointer to the newly allocated struct if no
errors were detected, and a null pointer in the case of error. The length $n is
factorized into a product of subtransforms, and the factors and their
trigonometric coefficients are stored in the wavetable. The trigonometric
coefficients are computed using direct calls to sin and cos, for accuracy.
Recursion relations could be used to compute the lookup table faster, but if an
application performs many FFTs of the same length then computing the wavetable
is a one-off overhead which does not affect the final throughput. The
wavetable structure can be used repeatedly for any transform of the same
length. The table is not modified by calls to any of the other FFT functions.
The appropriate type of wavetable must be used for forward real or inverse
half-complex transforms.
=item * C<gsl_fft_real_wavetable_free($wavetable)>
This function frees the memory associated with the wavetable $wavetable. The
wavetable can be freed if no further FFTs of the same length will be needed.
=item * C<gsl_fft_real_workspace_alloc($n)>
This function allocates a workspace for a real transform of length $n. The same
workspace can be used for both forward real and inverse halfcomplex transforms.
=item * C<gsl_fft_real_workspace_free($workspace)>
This function frees the memory associated with the workspace $workspace. The
workspace can be freed if no further FFTs of the same length will be needed.
=item * C<gsl_fft_real_transform >
=item * C<gsl_fft_real_unpack >
=back
This module also includes the following constants :
=over
=item * C<$gsl_fft_forward>
=item * C<$gsl_fft_backward>
=back
For more information on the functions, we refer you to the GSL official
documentation: L<http://www.gnu.org/software/gsl/manual/html_node/>
=head1 AUTHORS
Jonathan "Duke" Leto <jonathan@leto.net> and Thierry Moisan <thierry.moisan@gmail.com>
=head1 COPYRIGHT AND LICENSE
Copyright (C) 2008-2024 Jonathan "Duke" Leto and Thierry Moisan
This program is free software; you can redistribute it and/or modify it
under the same terms as Perl itself.
=cut
%}
|
