File: Mixed_002dradix-FFT-routines-for-real-data.html

package info (click to toggle)
gsl-ref-html 2.3-1
  • links: PTS
  • area: non-free
  • in suites: bullseye, buster, sid
  • size: 6,876 kB
  • ctags: 4,574
  • sloc: makefile: 35
file content (358 lines) | stat: -rw-r--r-- 16,745 bytes parent folder | download
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
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/html4/loose.dtd">
<html>
<!-- Copyright (C) 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014, 2015, 2016 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.3 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 the freedom to copy and modify this
GNU Manual." -->
<!-- Created by GNU Texinfo 5.1, http://www.gnu.org/software/texinfo/ -->
<head>
<title>GNU Scientific Library &ndash; Reference Manual: Mixed-radix FFT routines for real data</title>

<meta name="description" content="GNU Scientific Library &ndash; Reference Manual: Mixed-radix FFT routines for real data">
<meta name="keywords" content="GNU Scientific Library &ndash; Reference Manual: Mixed-radix FFT routines for real data">
<meta name="resource-type" content="document">
<meta name="distribution" content="global">
<meta name="Generator" content="makeinfo">
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<link href="index.html#Top" rel="start" title="Top">
<link href="Function-Index.html#Function-Index" rel="index" title="Function Index">
<link href="Fast-Fourier-Transforms.html#Fast-Fourier-Transforms" rel="up" title="Fast Fourier Transforms">
<link href="FFT-References-and-Further-Reading.html#FFT-References-and-Further-Reading" rel="next" title="FFT References and Further Reading">
<link href="Radix_002d2-FFT-routines-for-real-data.html#Radix_002d2-FFT-routines-for-real-data" rel="previous" title="Radix-2 FFT routines for real data">
<style type="text/css">
<!--
a.summary-letter {text-decoration: none}
blockquote.smallquotation {font-size: smaller}
div.display {margin-left: 3.2em}
div.example {margin-left: 3.2em}
div.indentedblock {margin-left: 3.2em}
div.lisp {margin-left: 3.2em}
div.smalldisplay {margin-left: 3.2em}
div.smallexample {margin-left: 3.2em}
div.smallindentedblock {margin-left: 3.2em; font-size: smaller}
div.smalllisp {margin-left: 3.2em}
kbd {font-style:oblique}
pre.display {font-family: inherit}
pre.format {font-family: inherit}
pre.menu-comment {font-family: serif}
pre.menu-preformatted {font-family: serif}
pre.smalldisplay {font-family: inherit; font-size: smaller}
pre.smallexample {font-size: smaller}
pre.smallformat {font-family: inherit; font-size: smaller}
pre.smalllisp {font-size: smaller}
span.nocodebreak {white-space:nowrap}
span.nolinebreak {white-space:nowrap}
span.roman {font-family:serif; font-weight:normal}
span.sansserif {font-family:sans-serif; font-weight:normal}
ul.no-bullet {list-style: none}
-->
</style>


</head>

<body lang="en" bgcolor="#FFFFFF" text="#000000" link="#0000FF" vlink="#800080" alink="#FF0000">
<a name="Mixed_002dradix-FFT-routines-for-real-data"></a>
<div class="header">
<p>
Next: <a href="FFT-References-and-Further-Reading.html#FFT-References-and-Further-Reading" accesskey="n" rel="next">FFT References and Further Reading</a>, Previous: <a href="Radix_002d2-FFT-routines-for-real-data.html#Radix_002d2-FFT-routines-for-real-data" accesskey="p" rel="previous">Radix-2 FFT routines for real data</a>, Up: <a href="Fast-Fourier-Transforms.html#Fast-Fourier-Transforms" accesskey="u" rel="up">Fast Fourier Transforms</a> &nbsp; [<a href="Function-Index.html#Function-Index" title="Index" rel="index">Index</a>]</p>
</div>
<hr>
<a name="Mixed_002dradix-FFT-routines-for-real-data-1"></a>
<h3 class="section">16.7 Mixed-radix FFT routines for real data</h3>
<a name="index-FFT-of-real-data_002c-mixed_002dradix-algorithm"></a>
<a name="index-Mixed_002dradix-FFT_002c-real-data"></a>

<p>This section describes mixed-radix FFT algorithms for real data.  The
mixed-radix functions work for FFTs of any length.  They are a
reimplementation of the real-FFT routines in the Fortran <small>FFTPACK</small> library
by Paul Swarztrauber.  The theory behind the algorithm is explained in
the article <cite>Fast Mixed-Radix Real Fourier Transforms</cite> by Clive
Temperton.  The routines here use the same indexing scheme and basic
algorithms as <small>FFTPACK</small>.
</p>
<p>The functions use the <small>FFTPACK</small> storage convention for half-complex
sequences.  In this convention the half-complex transform of a real
sequence is stored with frequencies in increasing order, starting at
zero, with the real and imaginary parts of each frequency in neighboring
locations.  When a value is known to be real the imaginary part is not
stored.  The imaginary part of the zero-frequency component is never
stored.  It is known to be zero (since the zero frequency component is
simply the sum of the input data (all real)).  For a sequence of even
length the imaginary part of the frequency <em>n/2</em> is not stored
either, since the symmetry 
<em>z_k = z_{n-k}^*</em> implies that this is
purely real too.
</p>
<p>The storage scheme is best shown by some examples.  The table below
shows the output for an odd-length sequence, <em>n=5</em>.  The two columns
give the correspondence between the 5 values in the half-complex
sequence returned by <code>gsl_fft_real_transform</code>, <var>halfcomplex</var>[] and the
values <var>complex</var>[] that would be returned if the same real input
sequence were passed to <code>gsl_fft_complex_backward</code> as a complex
sequence (with imaginary parts set to <code>0</code>),
</p>
<div class="example">
<pre class="example">complex[0].real  =  halfcomplex[0] 
complex[0].imag  =  0
complex[1].real  =  halfcomplex[1] 
complex[1].imag  =  halfcomplex[2]
complex[2].real  =  halfcomplex[3]
complex[2].imag  =  halfcomplex[4]
complex[3].real  =  halfcomplex[3]
complex[3].imag  = -halfcomplex[4]
complex[4].real  =  halfcomplex[1]
complex[4].imag  = -halfcomplex[2]
</pre></div>

<p>The upper elements of the <var>complex</var> array, <code>complex[3]</code> and
<code>complex[4]</code> are filled in using the symmetry condition.  The
imaginary part of the zero-frequency term <code>complex[0].imag</code> is
known to be zero by the symmetry.
</p>
<p>The next table shows the output for an even-length sequence, <em>n=6</em>.
In the even case there are two values which are purely real,
</p>
<div class="example">
<pre class="example">complex[0].real  =  halfcomplex[0]
complex[0].imag  =  0
complex[1].real  =  halfcomplex[1] 
complex[1].imag  =  halfcomplex[2] 
complex[2].real  =  halfcomplex[3] 
complex[2].imag  =  halfcomplex[4] 
complex[3].real  =  halfcomplex[5] 
complex[3].imag  =  0 
complex[4].real  =  halfcomplex[3] 
complex[4].imag  = -halfcomplex[4]
complex[5].real  =  halfcomplex[1] 
complex[5].imag  = -halfcomplex[2] 
</pre></div>

<p>The upper elements of the <var>complex</var> array, <code>complex[4]</code> and
<code>complex[5]</code> are filled in using the symmetry condition.  Both
<code>complex[0].imag</code> and <code>complex[3].imag</code> are known to be zero.
</p>
<p>All these functions are declared in the header files
<samp>gsl_fft_real.h</samp> and <samp>gsl_fft_halfcomplex.h</samp>.
</p>
<dl>
<dt><a name="index-gsl_005ffft_005freal_005fwavetable_005falloc"></a>Function: <em>gsl_fft_real_wavetable *</em> <strong>gsl_fft_real_wavetable_alloc</strong> <em>(size_t <var>n</var>)</em></dt>
<dt><a name="index-gsl_005ffft_005fhalfcomplex_005fwavetable_005falloc"></a>Function: <em>gsl_fft_halfcomplex_wavetable *</em> <strong>gsl_fft_halfcomplex_wavetable_alloc</strong> <em>(size_t <var>n</var>)</em></dt>
<dd><a name="index-gsl_005ffft_005freal_005fwavetable"></a>
<a name="index-gsl_005ffft_005fhalfcomplex_005fwavetable"></a>
<p>These functions prepare trigonometric lookup tables for an FFT of size
<em>n</em> 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 <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
computing the wavetable is a one-off overhead which does not affect the
final throughput.
</p>
<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 appropriate type of wavetable must be used for forward
real or inverse half-complex transforms.
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005ffft_005freal_005fwavetable_005ffree"></a>Function: <em>void</em> <strong>gsl_fft_real_wavetable_free</strong> <em>(gsl_fft_real_wavetable * <var>wavetable</var>)</em></dt>
<dt><a name="index-gsl_005ffft_005fhalfcomplex_005fwavetable_005ffree"></a>Function: <em>void</em> <strong>gsl_fft_halfcomplex_wavetable_free</strong> <em>(gsl_fft_halfcomplex_wavetable * <var>wavetable</var>)</em></dt>
<dd><p>These functions free 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></dd></dl>

<p>The mixed radix algorithms require additional working space to hold
the intermediate steps of the transform,
</p>
<dl>
<dt><a name="index-gsl_005ffft_005freal_005fworkspace_005falloc"></a>Function: <em>gsl_fft_real_workspace *</em> <strong>gsl_fft_real_workspace_alloc</strong> <em>(size_t <var>n</var>)</em></dt>
<dd><a name="index-gsl_005ffft_005freal_005fworkspace"></a>
<p>This function allocates a workspace for a real transform of length
<var>n</var>.  The same workspace can be used for both forward real and inverse
halfcomplex transforms.
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005ffft_005freal_005fworkspace_005ffree"></a>Function: <em>void</em> <strong>gsl_fft_real_workspace_free</strong> <em>(gsl_fft_real_workspace * <var>workspace</var>)</em></dt>
<dd><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></dd></dl>

<p>The following functions compute the transforms of real and half-complex
data,
</p>
<dl>
<dt><a name="index-gsl_005ffft_005freal_005ftransform"></a>Function: <em>int</em> <strong>gsl_fft_real_transform</strong> <em>(double <var>data</var>[], size_t <var>stride</var>, size_t <var>n</var>, const gsl_fft_real_wavetable * <var>wavetable</var>, gsl_fft_real_workspace * <var>work</var>)</em></dt>
<dt><a name="index-gsl_005ffft_005fhalfcomplex_005ftransform"></a>Function: <em>int</em> <strong>gsl_fft_halfcomplex_transform</strong> <em>(double <var>data</var>[], size_t <var>stride</var>, size_t <var>n</var>, const gsl_fft_halfcomplex_wavetable * <var>wavetable</var>, gsl_fft_real_workspace * <var>work</var>)</em></dt>
<dd><p>These functions compute the FFT of <var>data</var>, a real or half-complex
array of length <var>n</var>, using a mixed radix decimation-in-frequency
algorithm.  For <code>gsl_fft_real_transform</code> <var>data</var> is an array of
time-ordered real data.  For <code>gsl_fft_halfcomplex_transform</code>
<var>data</var> contains Fourier coefficients in the half-complex ordering
described above.  There is no restriction on the length <var>n</var>.
Efficient modules are provided for subtransforms of length 2, 3, 4 and
5.  Any remaining factors are computed with a slow, <em>O(n^2)</em>,
general-n module.  The caller must supply a <var>wavetable</var> containing
trigonometric lookup tables and a workspace <var>work</var>. 
</p></dd></dl>

<dl>
<dt><a name="index-gsl_005ffft_005freal_005funpack"></a>Function: <em>int</em> <strong>gsl_fft_real_unpack</strong> <em>(const double <var>real_coefficient</var>[], gsl_complex_packed_array <var>complex_coefficient</var>, size_t <var>stride</var>, size_t <var>n</var>)</em></dt>
<dd>
<p>This function converts a single real array, <var>real_coefficient</var> into
an equivalent complex array, <var>complex_coefficient</var>, (with imaginary
part set to zero), suitable for <code>gsl_fft_complex</code> routines.  The
algorithm for the conversion is simply,
</p>
<div class="example">
<pre class="example">for (i = 0; i &lt; n; i++)
  {
    complex_coefficient[i*stride].real 
      = real_coefficient[i*stride];
    complex_coefficient[i*stride].imag 
      = 0.0;
  }
</pre></div>
</dd></dl>

<dl>
<dt><a name="index-gsl_005ffft_005fhalfcomplex_005funpack"></a>Function: <em>int</em> <strong>gsl_fft_halfcomplex_unpack</strong> <em>(const double <var>halfcomplex_coefficient</var>[], gsl_complex_packed_array <var>complex_coefficient</var>, size_t <var>stride</var>, size_t <var>n</var>)</em></dt>
<dd>
<p>This function converts <var>halfcomplex_coefficient</var>, an array of
half-complex coefficients as returned by <code>gsl_fft_real_transform</code>, into an
ordinary complex array, <var>complex_coefficient</var>.  It fills in the
complex array using the symmetry 
<em>z_k = z_{n-k}^*</em>
to reconstruct the redundant elements.  The algorithm for the conversion
is,
</p>
<div class="example">
<pre class="example">complex_coefficient[0].real 
  = halfcomplex_coefficient[0];
complex_coefficient[0].imag 
  = 0.0;

for (i = 1; i &lt; n - i; i++)
  {
    double hc_real 
      = halfcomplex_coefficient[(2 * i - 1)*stride];
    double hc_imag 
      = halfcomplex_coefficient[(2 * i)*stride];
    complex_coefficient[i*stride].real = hc_real;
    complex_coefficient[i*stride].imag = hc_imag;
    complex_coefficient[(n - i)*stride].real = hc_real;
    complex_coefficient[(n - i)*stride].imag = -hc_imag;
  }

if (i == n - i)
  {
    complex_coefficient[i*stride].real 
      = halfcomplex_coefficient[(n - 1)*stride];
    complex_coefficient[i*stride].imag 
      = 0.0;
  }
</pre></div>
</dd></dl>


<p>Here is an example program using <code>gsl_fft_real_transform</code> and
<code>gsl_fft_halfcomplex_inverse</code>.  It generates a real signal in the
shape of a square pulse.  The pulse is Fourier transformed to frequency
space, and all but the lowest ten frequency components are removed from
the array of Fourier coefficients returned by
<code>gsl_fft_real_transform</code>.
</p>
<p>The remaining Fourier coefficients are transformed back to the
time-domain, to give a filtered version of the square pulse.  Since
Fourier coefficients are stored using the half-complex symmetry both
positive and negative frequencies are removed and the final filtered
signal is also real.
</p>
<div class="example">
<pre class="verbatim">#include &lt;stdio.h&gt;
#include &lt;math.h&gt;
#include &lt;gsl/gsl_errno.h&gt;
#include &lt;gsl/gsl_fft_real.h&gt;
#include &lt;gsl/gsl_fft_halfcomplex.h&gt;

int
main (void)
{
  int i, n = 100;
  double data[n];

  gsl_fft_real_wavetable * real;
  gsl_fft_halfcomplex_wavetable * hc;
  gsl_fft_real_workspace * work;

  for (i = 0; i &lt; n; i++)
    {
      data[i] = 0.0;
    }

  for (i = n / 3; i &lt; 2 * n / 3; i++)
    {
      data[i] = 1.0;
    }

  for (i = 0; i &lt; n; i++)
    {
      printf (&quot;%d: %e\n&quot;, i, data[i]);
    }
  printf (&quot;\n&quot;);

  work = gsl_fft_real_workspace_alloc (n);
  real = gsl_fft_real_wavetable_alloc (n);

  gsl_fft_real_transform (data, 1, n, 
                          real, work);

  gsl_fft_real_wavetable_free (real);

  for (i = 11; i &lt; n; i++)
    {
      data[i] = 0;
    }

  hc = gsl_fft_halfcomplex_wavetable_alloc (n);

  gsl_fft_halfcomplex_inverse (data, 1, n, 
                               hc, work);
  gsl_fft_halfcomplex_wavetable_free (hc);

  for (i = 0; i &lt; n; i++)
    {
      printf (&quot;%d: %e\n&quot;, i, data[i]);
    }

  gsl_fft_real_workspace_free (work);
  return 0;
}
</pre></div>


<hr>
<div class="header">
<p>
Next: <a href="FFT-References-and-Further-Reading.html#FFT-References-and-Further-Reading" accesskey="n" rel="next">FFT References and Further Reading</a>, Previous: <a href="Radix_002d2-FFT-routines-for-real-data.html#Radix_002d2-FFT-routines-for-real-data" accesskey="p" rel="previous">Radix-2 FFT routines for real data</a>, Up: <a href="Fast-Fourier-Transforms.html#Fast-Fourier-Transforms" accesskey="u" rel="up">Fast Fourier Transforms</a> &nbsp; [<a href="Function-Index.html#Function-Index" title="Index" rel="index">Index</a>]</p>
</div>



</body>
</html>