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