Audacity 3.2.0
RealFFTf.cpp
Go to the documentation of this file.
1/*
2* Program: REALFFTF.C
3* Author: Philip Van Baren
4* Date: 2 September 1993
5*
6* Description: These routines perform an FFT on real data to get a conjugate-symmetric
7* output, and an inverse FFT on conjugate-symmetric input to get a real
8* output sequence.
9*
10* This code is for floating point data.
11*
12* Modified 8/19/1998 by Philip Van Baren
13* - made the InitializeFFT and EndFFT routines take a structure
14* holding the length and pointers to the BitReversed and SinTable
15* tables.
16* Modified 5/23/2009 by Philip Van Baren
17* - Added GetFFT and ReleaseFFT routines to retain common SinTable
18* and BitReversed tables so they don't need to be reallocated
19* and recomputed on every call.
20* - Added Reorder* functions to undo the bit-reversal
21*
22* Copyright (C) 2009 Philip VanBaren
23*
24* This program is free software; you can redistribute it and/or modify
25* it under the terms of the GNU General Public License as published by
26* the Free Software Foundation; either version 2 of the License, or
27* (at your option) any later version.
28*
29* This program is distributed in the hope that it will be useful,
30* but WITHOUT ANY WARRANTY; without even the implied warranty of
31* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32* GNU General Public License for more details.
33*
34* You should have received a copy of the GNU General Public License
35* along with this program; if not, write to the Free Software
36* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
37*/
38
39#include "RealFFTf.h"
40
41#include <vector>
42#include <stdlib.h>
43#include <stdio.h>
44#include <math.h>
45
46#include <wx/thread.h>
47
48#ifndef M_PI
49#define M_PI 3.14159265358979323846 /* pi */
50#endif
51
52/*
53* Initialize the Sine table and Twiddle pointers (bit-reversed pointers)
54* for the FFT routine.
55*/
56HFFT InitializeFFT(size_t fftlen)
57{
58 int temp;
60
61 /*
62 * FFT size is only half the number of data points
63 * The full FFT output can be reconstructed from this FFT's output.
64 * (This optimization can be made since the data is real.)
65 */
66 h->Points = fftlen / 2;
67
68 h->SinTable.reinit(2*h->Points);
69
70 h->BitReversed.reinit(h->Points);
71
72 for(size_t i = 0; i < h->Points; i++)
73 {
74 temp = 0;
75 for(size_t mask = h->Points / 2; mask > 0; mask >>= 1)
76 temp = (temp >> 1) + (i & mask ? h->Points : 0);
77
78 h->BitReversed[i] = temp;
79 }
80
81 for(size_t i = 0; i < h->Points; i++)
82 {
83 h->SinTable[h->BitReversed[i] ]=(fft_type)-sin(2*M_PI*i/(2*h->Points));
84 h->SinTable[h->BitReversed[i]+1]=(fft_type)-cos(2*M_PI*i/(2*h->Points));
85 }
86
87#ifdef EXPERIMENTAL_EQ_SSE_THREADED
88 // NEW SSE FFT routines work on live data
89 for(size_t i = 0; i < 32; i++)
90 if((1 << i) & fftlen)
91 h->pow2Bits = i;
92#endif
93
94 return h;
95}
96
97enum : size_t { MAX_HFFT = 10 };
98
99// Maintain a pool:
100static std::vector< std::unique_ptr<FFTParam> > hFFTArray(MAX_HFFT);
101wxCriticalSection getFFTMutex;
102
103/* Get a handle to the FFT tables of the desired length */
104/* This version keeps common tables rather than allocating a NEW table every time */
105HFFT GetFFT(size_t fftlen)
106{
107 // To do: smarter policy about when to retain in the pool and when to
108 // allocate a unique instance.
109
110 wxCriticalSectionLocker locker{ getFFTMutex };
111
112 size_t h = 0;
113 auto n = fftlen/2;
114 auto size = hFFTArray.size();
115 for(;
116 (h < size) && hFFTArray[h] && (n != hFFTArray[h]->Points);
117 h++)
118 ;
119 if(h < size) {
120 if(hFFTArray[h] == NULL) {
121 hFFTArray[h].reset( InitializeFFT(fftlen).release() );
122 }
123 return HFFT{ hFFTArray[h].get() };
124 } else {
125 // All buffers used, so fall back to allocating a NEW set of tables
126 return InitializeFFT(fftlen);
127 }
128}
129
130/* Release a previously requested handle to the FFT tables */
132{
133 wxCriticalSectionLocker locker{ getFFTMutex };
134
135 auto it = hFFTArray.begin(), end = hFFTArray.end();
136 while (it != end && it->get() != hFFT)
137 ++it;
138 if ( it != end )
139 ;
140 else
141 delete hFFT;
142}
143
144/*
145* Forward FFT routine. Must call GetFFT(fftlen) first!
146*
147* Note: Output is BIT-REVERSED! so you must use the BitReversed to
148* get legible output, (i.e. Real_i = buffer[ h->BitReversed[i] ]
149* Imag_i = buffer[ h->BitReversed[i]+1 ] )
150* Input is in normal order.
151*
152* Output buffer[0] is the DC bin, and output buffer[1] is the Fs/2 bin
153* - this can be done because both values will always be real only
154* - this allows us to not have to allocate an extra complex value for the Fs/2 bin
155*
156* Note: The scaling on this is done according to the standard FFT definition,
157* so a unit amplitude DC signal will output an amplitude of (N)
158* (Older revisions would progressively scale the input, so the output
159* values would be similar in amplitude to the input values, which is
160* good when using fixed point arithmetic)
161*/
162void RealFFTf(fft_type *buffer, const FFTParam *h)
163{
164 fft_type *A,*B;
165 const fft_type *sptr;
166 const fft_type *endptr1,*endptr2;
167 const int *br1,*br2;
168 fft_type HRplus,HRminus,HIplus,HIminus;
169 fft_type v1,v2,sin,cos;
170
171 auto ButterfliesPerGroup = h->Points/2;
172
173 /*
174 * Butterfly:
175 * Ain-----Aout
176 * \ /
177 * / \
178 * Bin-----Bout
179 */
180
181 endptr1 = buffer + h->Points * 2;
182
183 while(ButterfliesPerGroup > 0)
184 {
185 A = buffer;
186 B = buffer + ButterfliesPerGroup * 2;
187 sptr = h->SinTable.get();
188
189 while(A < endptr1)
190 {
191 sin = *sptr;
192 cos = *(sptr+1);
193 endptr2 = B;
194 while(A < endptr2)
195 {
196 v1 = *B * cos + *(B + 1) * sin;
197 v2 = *B * sin - *(B + 1) * cos;
198 *B = (*A + v1);
199 *(A++) = *(B++) - 2 * v1;
200 *B = (*A - v2);
201 *(A++) = *(B++) + 2 * v2;
202 }
203 A = B;
204 B += ButterfliesPerGroup * 2;
205 sptr += 2;
206 }
207 ButterfliesPerGroup >>= 1;
208 }
209 /* Massage output to get the output for a real input sequence. */
210 br1 = h->BitReversed.get() + 1;
211 br2 = h->BitReversed.get() + h->Points - 1;
212
213 while(br1<br2)
214 {
215 sin=h->SinTable[*br1];
216 cos=h->SinTable[*br1+1];
217 A=buffer+*br1;
218 B=buffer+*br2;
219 HRplus = (HRminus = *A - *B ) + (*B * 2);
220 HIplus = (HIminus = *(A+1) - *(B+1)) + (*(B+1) * 2);
221 v1 = (sin*HRminus - cos*HIplus);
222 v2 = (cos*HRminus + sin*HIplus);
223 *A = (HRplus + v1) * (fft_type)0.5;
224 *B = *A - v1;
225 *(A+1) = (HIminus + v2) * (fft_type)0.5;
226 *(B+1) = *(A+1) - HIminus;
227
228 br1++;
229 br2--;
230 }
231 /* Handle the center bin (just need a conjugate) */
232 A=buffer+*br1+1;
233 *A=-*A;
234 /* Handle DC bin separately - and ignore the Fs/2 bin
235 buffer[0]+=buffer[1];
236 buffer[1]=(fft_type)0;*/
237 /* Handle DC and Fs/2 bins separately */
238 /* Put the Fs/2 value into the imaginary part of the DC bin */
239 v1=buffer[0]-buffer[1];
240 buffer[0]+=buffer[1];
241 buffer[1]=v1;
242}
243
244
245/* Description: This routine performs an inverse FFT to real data.
246* This code is for floating point data.
247*
248* Note: Output is BIT-REVERSED! so you must use the BitReversed to
249* get legible output, (i.e. wave[2*i] = buffer[ BitReversed[i] ]
250* wave[2*i+1] = buffer[ BitReversed[i]+1 ] )
251* Input is in normal order, interleaved (real,imaginary) complex data
252* You must call GetFFT(fftlen) first to initialize some buffers!
253*
254* Input buffer[0] is the DC bin, and input buffer[1] is the Fs/2 bin
255* - this can be done because both values will always be real only
256* - this allows us to not have to allocate an extra complex value for the Fs/2 bin
257*
258* Note: The scaling on this is done according to the standard FFT definition,
259* so a unit amplitude DC signal will output an amplitude of (N)
260* (Older revisions would progressively scale the input, so the output
261* values would be similar in amplitude to the input values, which is
262* good when using fixed point arithmetic)
263*/
264void InverseRealFFTf(fft_type *buffer, const FFTParam *h)
265{
266 fft_type *A,*B;
267 const fft_type *sptr;
268 const fft_type *endptr1,*endptr2;
269 const int *br1;
270 fft_type HRplus,HRminus,HIplus,HIminus;
271 fft_type v1,v2,sin,cos;
272
273 auto ButterfliesPerGroup = h->Points / 2;
274
275 /* Massage input to get the input for a real output sequence. */
276 A = buffer + 2;
277 B = buffer + h->Points * 2 - 2;
278 br1 = h->BitReversed.get() + 1;
279 while(A<B)
280 {
281 sin=h->SinTable[*br1];
282 cos=h->SinTable[*br1+1];
283 HRplus = (HRminus = *A - *B ) + (*B * 2);
284 HIplus = (HIminus = *(A+1) - *(B+1)) + (*(B+1) * 2);
285 v1 = (sin*HRminus + cos*HIplus);
286 v2 = (cos*HRminus - sin*HIplus);
287 *A = (HRplus + v1) * (fft_type)0.5;
288 *B = *A - v1;
289 *(A+1) = (HIminus - v2) * (fft_type)0.5;
290 *(B+1) = *(A+1) - HIminus;
291
292 A+=2;
293 B-=2;
294 br1++;
295 }
296 /* Handle center bin (just need conjugate) */
297 *(A+1)=-*(A+1);
298 /* Handle DC bin separately - this ignores any Fs/2 component
299 buffer[1]=buffer[0]=buffer[0]/2;*/
300 /* Handle DC and Fs/2 bins specially */
301 /* The DC bin is passed in as the real part of the DC complex value */
302 /* The Fs/2 bin is passed in as the imaginary part of the DC complex value */
303 /* (v1+v2) = buffer[0] == the DC component */
304 /* (v1-v2) = buffer[1] == the Fs/2 component */
305 v1=0.5f*(buffer[0]+buffer[1]);
306 v2=0.5f*(buffer[0]-buffer[1]);
307 buffer[0]=v1;
308 buffer[1]=v2;
309
310 /*
311 * Butterfly:
312 * Ain-----Aout
313 * \ /
314 * / \
315 * Bin-----Bout
316 */
317
318 endptr1 = buffer + h->Points * 2;
319
320 while(ButterfliesPerGroup > 0)
321 {
322 A = buffer;
323 B = buffer + ButterfliesPerGroup * 2;
324 sptr = h->SinTable.get();
325
326 while(A < endptr1)
327 {
328 sin = *(sptr++);
329 cos = *(sptr++);
330 endptr2 = B;
331 while(A < endptr2)
332 {
333 v1 = *B * cos - *(B + 1) * sin;
334 v2 = *B * sin + *(B + 1) * cos;
335 *B = (*A + v1) * (fft_type)0.5;
336 *(A++) = *(B++) - v1;
337 *B = (*A + v2) * (fft_type)0.5;
338 *(A++) = *(B++) - v2;
339 }
340 A = B;
341 B += ButterfliesPerGroup * 2;
342 }
343 ButterfliesPerGroup >>= 1;
344 }
345}
346
347void ReorderToFreq(const FFTParam *hFFT, const fft_type *buffer,
348 fft_type *RealOut, fft_type *ImagOut)
349{
350 // Copy the data into the real and imaginary outputs
351 for(size_t i = 1; i < hFFT->Points; i++) {
352 RealOut[i] = buffer[hFFT->BitReversed[i] ];
353 ImagOut[i] = buffer[hFFT->BitReversed[i]+1];
354 }
355 RealOut[0] = buffer[0]; // DC component
356 ImagOut[0] = 0;
357 RealOut[hFFT->Points] = buffer[1]; // Fs/2 component
358 ImagOut[hFFT->Points] = 0;
359}
360
361void ReorderToTime(const FFTParam *hFFT, const fft_type *buffer, fft_type *TimeOut)
362{
363 // Copy the data into the real outputs
364 for(size_t i = 0; i < hFFT->Points; i++) {
365 TimeOut[i*2 ]=buffer[hFFT->BitReversed[i] ];
366 TimeOut[i*2+1]=buffer[hFFT->BitReversed[i]+1];
367 }
368}
#define safenew
Definition: MemoryX.h:10
float fft_type
Definition: RealFFTf48x.h:6
void RealFFTf(fft_type *buffer, const FFTParam *h)
Definition: RealFFTf.cpp:162
static std::vector< std::unique_ptr< FFTParam > > hFFTArray(MAX_HFFT)
void InverseRealFFTf(fft_type *buffer, const FFTParam *h)
Definition: RealFFTf.cpp:264
HFFT InitializeFFT(size_t fftlen)
Definition: RealFFTf.cpp:56
void ReorderToTime(const FFTParam *hFFT, const fft_type *buffer, fft_type *TimeOut)
Definition: RealFFTf.cpp:361
@ MAX_HFFT
Definition: RealFFTf.cpp:97
HFFT GetFFT(size_t fftlen)
Definition: RealFFTf.cpp:105
wxCriticalSection getFFTMutex
Definition: RealFFTf.cpp:101
void ReorderToFreq(const FFTParam *hFFT, const fft_type *buffer, fft_type *RealOut, fft_type *ImagOut)
Definition: RealFFTf.cpp:347
#define M_PI
Definition: RealFFTf.cpp:49
std::unique_ptr< FFTParam, FFTDeleter > HFFT
Definition: RealFFTf.h:22
#define A(N)
Definition: ToChars.cpp:62
auto end(const Ptr< Type, BaseDeleter > &p)
Enables range-for.
Definition: PackedArray.h:159
void operator()(FFTParam *p) const
Definition: RealFFTf.cpp:131
size_t Points
Definition: RealFFTf.h:10
ArrayOf< int > BitReversed
Definition: RealFFTf.h:8
ArrayOf< fft_type > SinTable
Definition: RealFFTf.h:9