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
|
<html lang="en">
<head>
<title>Mixed-radix FFT routines for complex data - GNU Scientific Library -- Reference Manual</title>
<meta http-equiv="Content-Type" content="text/html">
<meta name="description" content="GNU Scientific Library -- Reference Manual">
<meta name="generator" content="makeinfo 4.8">
<link title="Top" rel="start" href="index.html#Top">
<link rel="up" href="Fast-Fourier-Transforms.html" title="Fast Fourier Transforms">
<link rel="prev" href="Radix_002d2-FFT-routines-for-complex-data.html" title="Radix-2 FFT routines for complex data">
<link rel="next" href="Overview-of-real-data-FFTs.html" title="Overview of real data FFTs">
<link href="http://www.gnu.org/software/texinfo/" rel="generator-home" title="Texinfo Homepage">
<!--
Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007 The GSL Team.
Permission is granted to copy, distribute and/or modify this document
under the terms of the GNU Free Documentation License, Version 1.2 or
any later version published by the Free Software Foundation; with the
Invariant Sections being ``GNU General Public License'' and ``Free Software
Needs Free Documentation'', the Front-Cover text being ``A GNU Manual'',
and with the Back-Cover Text being (a) (see below). A copy of the
license is included in the section entitled ``GNU Free Documentation
License''.
(a) The Back-Cover Text is: ``You have freedom to copy and modify this
GNU Manual, like GNU software.''-->
<meta http-equiv="Content-Style-Type" content="text/css">
<style type="text/css"><!--
pre.display { font-family:inherit }
pre.format { font-family:inherit }
pre.smalldisplay { font-family:inherit; font-size:smaller }
pre.smallformat { font-family:inherit; font-size:smaller }
pre.smallexample { font-size:smaller }
pre.smalllisp { font-size:smaller }
span.sc { font-variant:small-caps }
span.roman { font-family:serif; font-weight:normal; }
span.sansserif { font-family:sans-serif; font-weight:normal; }
--></style>
</head>
<body>
<div class="node">
<p>
<a name="Mixed-radix-FFT-routines-for-complex-data"></a>
<a name="Mixed_002dradix-FFT-routines-for-complex-data"></a>
Next: <a rel="next" accesskey="n" href="Overview-of-real-data-FFTs.html">Overview of real data FFTs</a>,
Previous: <a rel="previous" accesskey="p" href="Radix_002d2-FFT-routines-for-complex-data.html">Radix-2 FFT routines for complex data</a>,
Up: <a rel="up" accesskey="u" href="Fast-Fourier-Transforms.html">Fast Fourier Transforms</a>
<hr>
</div>
<h3 class="section">15.4 Mixed-radix FFT routines for complex data</h3>
<p><a name="index-FFT-of-complex-data_002c-mixed_002dradix-algorithm-1407"></a><a name="index-Mixed_002dradix-FFT_002c-complex-data-1408"></a>
This section describes mixed-radix FFT algorithms for complex data. The
mixed-radix functions work for FFTs of any length. They are a
reimplementation of Paul Swarztrauber's Fortran <span class="sc">fftpack</span> library.
The theory is explained in the review article <cite>Self-sorting
Mixed-radix FFTs</cite> by Clive Temperton. The routines here use the same
indexing scheme and basic algorithms as <span class="sc">fftpack</span>.
<p>The mixed-radix algorithm is based on sub-transform modules—highly
optimized small length FFTs which are combined to create larger FFTs.
There are efficient modules for factors of 2, 3, 4, 5, 6 and 7. The
modules for the composite factors of 4 and 6 are faster than combining
the modules for 2*2 and 2*3.
<p>For factors which are not implemented as modules there is a fall-back to
a general length-n module which uses Singleton's method for
efficiently computing a DFT. This module is O(n^2), and slower
than a dedicated module would be but works for any length n. Of
course, lengths which use the general length-n module will still
be factorized as much as possible. For example, a length of 143 will be
factorized into 11*13. Large prime factors are the worst case
scenario, e.g. as found in n=2*3*99991, and should be avoided
because their O(n^2) scaling will dominate the run-time (consult
the document <cite>GSL FFT Algorithms</cite> included in the GSL distribution
if you encounter this problem).
<p>The mixed-radix initialization function <code>gsl_fft_complex_wavetable_alloc</code>
returns the list of factors chosen by the library for a given length
N. It can be used to check how well the length has been
factorized, and estimate the run-time. To a first approximation the
run-time scales as N \sum f_i, where the f_i are the
factors of N. For programs under user control you may wish to
issue a warning that the transform will be slow when the length is
poorly factorized. If you frequently encounter data lengths which
cannot be factorized using the existing small-prime modules consult
<cite>GSL FFT Algorithms</cite> for details on adding support for other
factors.
<!-- First, the space for the trigonometric lookup tables and scratch area is -->
<!-- allocated by a call to one of the @code{alloc} functions. We -->
<!-- call the combination of factorization, scratch space and trigonometric -->
<!-- lookup arrays a @dfn{wavetable}. It contains the sine and cosine -->
<!-- waveforms for the all the frequencies that will be used in the FFT. -->
<!-- The wavetable is initialized by a call to the corresponding @code{init} -->
<!-- function. It factorizes the data length, using the implemented -->
<!-- subtransforms as preferred factors wherever possible. The trigonometric -->
<!-- lookup table for the chosen factorization is also computed. -->
<!-- An FFT is computed by a call to one of the @code{forward}, -->
<!-- @code{backward} or @code{inverse} functions, with the data, length and -->
<!-- wavetable as arguments. -->
<p>All the functions described in this section are declared in the header
file <samp><span class="file">gsl_fft_complex.h</span></samp>.
<div class="defun">
— Function: gsl_fft_complex_wavetable * <b>gsl_fft_complex_wavetable_alloc</b> (<var>size_t n</var>)<var><a name="index-gsl_005ffft_005fcomplex_005fwavetable_005falloc-1409"></a></var><br>
<blockquote><p>This function prepares a trigonometric lookup table for a complex FFT of
length <var>n</var>. The function returns a pointer to the newly allocated
<code>gsl_fft_complex_wavetable</code> if no errors were detected, and a null
pointer in the case of error. The length <var>n</var> 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 <code>sin</code> and <code>cos</code>, 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.
<p>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.
</p></blockquote></div>
<div class="defun">
— Function: void <b>gsl_fft_complex_wavetable_free</b> (<var>gsl_fft_complex_wavetable * wavetable</var>)<var><a name="index-gsl_005ffft_005fcomplex_005fwavetable_005ffree-1410"></a></var><br>
<blockquote><p>This function frees the memory associated with the wavetable
<var>wavetable</var>. The wavetable can be freed if no further FFTs of the
same length will be needed.
</p></blockquote></div>
<p class="noindent">These functions operate on a <code>gsl_fft_complex_wavetable</code> structure
which contains internal parameters for the FFT. It is not necessary to
set any of the components directly but it can sometimes be useful to
examine them. For example, the chosen factorization of the FFT length
is given and can be used to provide an estimate of the run-time or
numerical error. The wavetable structure is declared in the header file
<samp><span class="file">gsl_fft_complex.h</span></samp>.
<div class="defun">
— Data Type: <b>gsl_fft_complex_wavetable</b><var><a name="index-gsl_005ffft_005fcomplex_005fwavetable-1411"></a></var><br>
<blockquote><p>This is a structure that holds the factorization and trigonometric
lookup tables for the mixed radix fft algorithm. It has the following
components:
<dl>
<dt><code>size_t n</code><dd>This is the number of complex data points
<br><dt><code>size_t nf</code><dd>This is the number of factors that the length <code>n</code> was decomposed into.
<br><dt><code>size_t factor[64]</code><dd>This is the array of factors. Only the first <code>nf</code> elements are
used.
<!-- (FIXME: This is a fixed length array and therefore probably in -->
<!-- violation of the GNU Coding Standards). -->
<br><dt><code>gsl_complex * trig</code><dd>This is a pointer to a preallocated trigonometric lookup table of
<code>n</code> complex elements.
<br><dt><code>gsl_complex * twiddle[64]</code><dd>This is an array of pointers into <code>trig</code>, giving the twiddle
factors for each pass.
</dl>
</p></blockquote></div>
<p class="noindent">The mixed radix algorithms require additional working space to hold
the intermediate steps of the transform.
<div class="defun">
— Function: gsl_fft_complex_workspace * <b>gsl_fft_complex_workspace_alloc</b> (<var>size_t n</var>)<var><a name="index-gsl_005ffft_005fcomplex_005fworkspace_005falloc-1412"></a></var><br>
<blockquote><p>This function allocates a workspace for a complex transform of length
<var>n</var>.
</p></blockquote></div>
<div class="defun">
— Function: void <b>gsl_fft_complex_workspace_free</b> (<var>gsl_fft_complex_workspace * workspace</var>)<var><a name="index-gsl_005ffft_005fcomplex_005fworkspace_005ffree-1413"></a></var><br>
<blockquote><p>This function frees the memory associated with the workspace
<var>workspace</var>. The workspace can be freed if no further FFTs of the
same length will be needed.
</p></blockquote></div>
<!-- @deftp {Data Type} gsl_fft_complex_workspace -->
<!-- This is a structure that holds the workspace for the mixed radix fft -->
<!-- algorithm. It has the following components: -->
<!-- @table @code -->
<!-- @item gsl_complex * scratch -->
<!-- This is a pointer to a workspace of @code{n} complex elements, -->
<!-- capable of holding intermediate copies of the original data set. -->
<!-- @end table -->
<!-- @end deftp -->
<p class="noindent">The following functions compute the transform,
<div class="defun">
— Function: int <b>gsl_fft_complex_forward</b> (<var>gsl_complex_packed_array data, size_t stride, size_t n, const gsl_fft_complex_wavetable * wavetable, gsl_fft_complex_workspace * work</var>)<var><a name="index-gsl_005ffft_005fcomplex_005fforward-1414"></a></var><br>
— Function: int <b>gsl_fft_complex_transform</b> (<var>gsl_complex_packed_array data, size_t stride, size_t n, const gsl_fft_complex_wavetable * wavetable, gsl_fft_complex_workspace * work, gsl_fft_direction sign</var>)<var><a name="index-gsl_005ffft_005fcomplex_005ftransform-1415"></a></var><br>
— Function: int <b>gsl_fft_complex_backward</b> (<var>gsl_complex_packed_array data, size_t stride, size_t n, const gsl_fft_complex_wavetable * wavetable, gsl_fft_complex_workspace * work</var>)<var><a name="index-gsl_005ffft_005fcomplex_005fbackward-1416"></a></var><br>
— Function: int <b>gsl_fft_complex_inverse</b> (<var>gsl_complex_packed_array data, size_t stride, size_t n, const gsl_fft_complex_wavetable * wavetable, gsl_fft_complex_workspace * work</var>)<var><a name="index-gsl_005ffft_005fcomplex_005finverse-1417"></a></var><br>
<blockquote>
<p>These functions compute forward, backward and inverse FFTs of length
<var>n</var> with stride <var>stride</var>, on the packed complex array
<var>data</var>, using a mixed radix decimation-in-frequency algorithm.
There is no restriction on the length <var>n</var>. Efficient modules are
provided for subtransforms of length 2, 3, 4, 5, 6 and 7. Any remaining
factors are computed with a slow, O(n^2), general-n
module. The caller must supply a <var>wavetable</var> containing the
trigonometric lookup tables and a workspace <var>work</var>. For the
<code>transform</code> version of the function the <var>sign</var> argument can be
either <code>forward</code> (-1) or <code>backward</code> (+1).
<p>The functions return a value of <code>0</code> if no errors were detected. The
following <code>gsl_errno</code> conditions are defined for these functions:
<dl>
<dt><code>GSL_EDOM</code><dd>The length of the data <var>n</var> is not a positive integer (i.e. <var>n</var>
is zero).
<br><dt><code>GSL_EINVAL</code><dd>The length of the data <var>n</var> and the length used to compute the given
<var>wavetable</var> do not match.
</dl>
</p></blockquote></div>
<!-- @node Example of using mixed-radix FFT routines for complex data -->
<!-- @subsection Example of using mixed-radix FFT routines for complex data -->
<p>Here is an example program which computes the FFT of a short pulse in a
sample of length 630 (=2*3*3*5*7) using the mixed-radix
algorithm.
<pre class="example"><pre class="verbatim"> #include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fft_complex.h>
#define REAL(z,i) ((z)[2*(i)])
#define IMAG(z,i) ((z)[2*(i)+1])
int
main (void)
{
int i;
const int n = 630;
double data[2*n];
gsl_fft_complex_wavetable * wavetable;
gsl_fft_complex_workspace * workspace;
for (i = 0; i < n; i++)
{
REAL(data,i) = 0.0;
IMAG(data,i) = 0.0;
}
data[0] = 1.0;
for (i = 1; i <= 10; i++)
{
REAL(data,i) = REAL(data,n-i) = 1.0;
}
for (i = 0; i < n; i++)
{
printf ("%d: %e %e\n", i, REAL(data,i),
IMAG(data,i));
}
printf ("\n");
wavetable = gsl_fft_complex_wavetable_alloc (n);
workspace = gsl_fft_complex_workspace_alloc (n);
for (i = 0; i < wavetable->nf; i++)
{
printf ("# factor %d: %d\n", i,
wavetable->factor[i]);
}
gsl_fft_complex_forward (data, 1, n,
wavetable, workspace);
for (i = 0; i < n; i++)
{
printf ("%d: %e %e\n", i, REAL(data,i),
IMAG(data,i));
}
gsl_fft_complex_wavetable_free (wavetable);
gsl_fft_complex_workspace_free (workspace);
return 0;
}
</pre></pre>
<p class="noindent">Note that we have assumed that the program is using the default
<code>gsl</code> error handler (which calls <code>abort</code> for any errors). If
you are not using a safe error handler you would need to check the
return status of all the <code>gsl</code> routines.
<hr>The GNU Scientific Library - a free numerical library licensed under the GNU GPL<br>Back to the <a href="/software/gsl/">GNU Scientific Library Homepage</a></body></html>
|
