Audacity  3.0.3
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 */
56 HFFT InitializeFFT(size_t fftlen)
57 {
58  int temp;
59  HFFT h{ safenew FFTParam };
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 
97 enum : size_t { MAX_HFFT = 10 };
98 
99 // Maintain a pool:
100 static std::vector< std::unique_ptr<FFTParam> > hFFTArray(MAX_HFFT);
101 wxCriticalSection 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 */
105 HFFT 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 */
162 void 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 */
264 void 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 
347 void 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 
361 void 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 }
size
size_t size
Definition: ffmpeg-2.3.6-single-header.h:412
MAX_HFFT
@ MAX_HFFT
Definition: RealFFTf.cpp:97
FFTParam::SinTable
ArrayOf< fft_type > SinTable
Definition: RealFFTf.h:9
RealFFTf
void RealFFTf(fft_type *buffer, const FFTParam *h)
Definition: RealFFTf.cpp:162
RealFFTf.h
ReorderToTime
void ReorderToTime(const FFTParam *hFFT, const fft_type *buffer, fft_type *TimeOut)
Definition: RealFFTf.cpp:361
M_PI
#define M_PI
Definition: RealFFTf.cpp:49
InitializeFFT
HFFT InitializeFFT(size_t fftlen)
Definition: RealFFTf.cpp:56
FFTDeleter::operator()
void operator()(FFTParam *p) const
Definition: RealFFTf.cpp:131
InverseRealFFTf
void InverseRealFFTf(fft_type *buffer, const FFTParam *h)
Definition: RealFFTf.cpp:264
FFTParam::Points
size_t Points
Definition: RealFFTf.h:10
FFTParam
Definition: RealFFTf.h:7
hFFTArray
static std::vector< std::unique_ptr< FFTParam > > hFFTArray(MAX_HFFT)
HFFT
std::unique_ptr< FFTParam, FFTDeleter > HFFT
Definition: RealFFTf.h:22
FFTParam::BitReversed
ArrayOf< int > BitReversed
Definition: RealFFTf.h:8
GetFFT
HFFT GetFFT(size_t fftlen)
Definition: RealFFTf.cpp:105
getFFTMutex
wxCriticalSection getFFTMutex
Definition: RealFFTf.cpp:101
safenew
#define safenew
Definition: MemoryX.h:10
ReorderToFreq
void ReorderToFreq(const FFTParam *hFFT, const fft_type *buffer, fft_type *RealOut, fft_type *ImagOut)
Definition: RealFFTf.cpp:347
fft_type
float fft_type
Definition: RealFFTf48x.h:6