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