1 /**********************************************************************
2
3 FFT.cpp
4
5 Dominic Mazzoni
6
7 September 2000
8
9 *******************************************************************//*!
10
11 \file FFT.cpp
12 \brief Fast Fourier Transform routines.
13
14 This file contains a few FFT routines, including a real-FFT
15 routine that is almost twice as fast as a normal complex FFT,
16 and a power spectrum routine when you know you don't care
17 about phase information.
18
19 Some of this code was based on a free implementation of an FFT
20 by Don Cross, available on the web at:
21
22 http://www.intersrv.com/~dcross/fft.html
23
24 The basic algorithm for his code was based on Numerican Recipes
25 in Fortran. I optimized his code further by reducing array
26 accesses, caching the bit reversal table, and eliminating
27 float-to-double conversions, and I added the routines to
28 calculate a real FFT and a real power spectrum.
29
30 *//*******************************************************************/
31 /*
32 Salvo Ventura - November 2006
33 Added more window functions:
34 * 4: Blackman
35 * 5: Blackman-Harris
36 * 6: Welch
37 * 7: Gaussian(a=2.5)
38 * 8: Gaussian(a=3.5)
39 * 9: Gaussian(a=4.5)
40 */
41
42 #include "FFT.h"
43
44 #include "Internat.h"
45
46 #include "SampleFormat.h"
47
48 #include <wx/wxcrtvararg.h>
49 #include <wx/intl.h>
50 #include <stdlib.h>
51 #include <stdio.h>
52 #include <math.h>
53
54 #include "RealFFTf.h"
55
56 static ArraysOf<int> gFFTBitTable;
57 static const size_t MaxFastBits = 16;
58
59 /* Declare Static functions */
60 static void InitFFT();
61
IsPowerOfTwo(size_t x)62 static bool IsPowerOfTwo(size_t x)
63 {
64 if (x < 2)
65 return false;
66
67 if (x & (x - 1)) /* Thanks to 'byang' for this cute trick! */
68 return false;
69
70 return true;
71 }
72
NumberOfBitsNeeded(size_t PowerOfTwo)73 static size_t NumberOfBitsNeeded(size_t PowerOfTwo)
74 {
75 if (PowerOfTwo < 2) {
76 wxFprintf(stderr, "Error: FFT called with size %ld\n", PowerOfTwo);
77 exit(1);
78 }
79
80 size_t i = 0;
81 while (PowerOfTwo > 1)
82 PowerOfTwo >>= 1, ++i;
83
84 return i;
85 }
86
ReverseBits(size_t index,size_t NumBits)87 int ReverseBits(size_t index, size_t NumBits)
88 {
89 size_t i, rev;
90
91 for (i = rev = 0; i < NumBits; i++) {
92 rev = (rev << 1) | (index & 1);
93 index >>= 1;
94 }
95
96 return rev;
97 }
98
InitFFT()99 void InitFFT()
100 {
101 gFFTBitTable.reinit(MaxFastBits);
102
103 size_t len = 2;
104 for (size_t b = 1; b <= MaxFastBits; b++) {
105 auto &array = gFFTBitTable[b - 1];
106 array.reinit(len);
107 for (size_t i = 0; i < len; i++)
108 array[i] = ReverseBits(i, b);
109
110 len <<= 1;
111 }
112 }
113
DeinitFFT()114 void DeinitFFT()
115 {
116 gFFTBitTable.reset();
117 }
118
FastReverseBits(size_t i,size_t NumBits)119 static inline size_t FastReverseBits(size_t i, size_t NumBits)
120 {
121 if (NumBits <= MaxFastBits)
122 return gFFTBitTable[NumBits - 1][i];
123 else
124 return ReverseBits(i, NumBits);
125 }
126
127 /*
128 * Complex Fast Fourier Transform
129 */
130
FFT(size_t NumSamples,bool InverseTransform,const float * RealIn,const float * ImagIn,float * RealOut,float * ImagOut)131 void FFT(size_t NumSamples,
132 bool InverseTransform,
133 const float *RealIn, const float *ImagIn,
134 float *RealOut, float *ImagOut)
135 {
136 double angle_numerator = 2.0 * M_PI;
137 double tr, ti; /* temp real, temp imaginary */
138
139 if (!IsPowerOfTwo(NumSamples)) {
140 wxFprintf(stderr, "%ld is not a power of two\n", NumSamples);
141 exit(1);
142 }
143
144 if (!gFFTBitTable)
145 InitFFT();
146
147 if (!InverseTransform)
148 angle_numerator = -angle_numerator;
149
150 /* Number of bits needed to store indices */
151 auto NumBits = NumberOfBitsNeeded(NumSamples);
152
153 /*
154 ** Do simultaneous data copy and bit-reversal ordering into outputs...
155 */
156
157 for (size_t i = 0; i < NumSamples; i++) {
158 auto j = FastReverseBits(i, NumBits);
159 RealOut[j] = RealIn[i];
160 ImagOut[j] = (ImagIn == NULL) ? 0.0 : ImagIn[i];
161 }
162
163 /*
164 ** Do the FFT itself...
165 */
166
167 size_t BlockEnd = 1;
168 for (size_t BlockSize = 2; BlockSize <= NumSamples; BlockSize <<= 1) {
169
170 double delta_angle = angle_numerator / (double) BlockSize;
171
172 double sm2 = sin(-2 * delta_angle);
173 double sm1 = sin(-delta_angle);
174 double cm2 = cos(-2 * delta_angle);
175 double cm1 = cos(-delta_angle);
176 double w = 2 * cm1;
177 double ar0, ar1, ar2, ai0, ai1, ai2;
178
179 for (size_t i = 0; i < NumSamples; i += BlockSize) {
180 ar2 = cm2;
181 ar1 = cm1;
182
183 ai2 = sm2;
184 ai1 = sm1;
185
186 for (size_t j = i, n = 0; n < BlockEnd; j++, n++) {
187 ar0 = w * ar1 - ar2;
188 ar2 = ar1;
189 ar1 = ar0;
190
191 ai0 = w * ai1 - ai2;
192 ai2 = ai1;
193 ai1 = ai0;
194
195 size_t k = j + BlockEnd;
196 tr = ar0 * RealOut[k] - ai0 * ImagOut[k];
197 ti = ar0 * ImagOut[k] + ai0 * RealOut[k];
198
199 RealOut[k] = RealOut[j] - tr;
200 ImagOut[k] = ImagOut[j] - ti;
201
202 RealOut[j] += tr;
203 ImagOut[j] += ti;
204 }
205 }
206
207 BlockEnd = BlockSize;
208 }
209
210 /*
211 ** Need to normalize if inverse transform...
212 */
213
214 if (InverseTransform) {
215 float denom = (float) NumSamples;
216
217 for (size_t i = 0; i < NumSamples; i++) {
218 RealOut[i] /= denom;
219 ImagOut[i] /= denom;
220 }
221 }
222 }
223
224 /*
225 * Real Fast Fourier Transform
226 *
227 * This is merely a wrapper of RealFFTf() from RealFFTf.h.
228 */
229
RealFFT(size_t NumSamples,const float * RealIn,float * RealOut,float * ImagOut)230 void RealFFT(size_t NumSamples, const float *RealIn, float *RealOut, float *ImagOut)
231 {
232 auto hFFT = GetFFT(NumSamples);
233 Floats pFFT{ NumSamples };
234 // Copy the data into the processing buffer
235 for(size_t i = 0; i < NumSamples; i++)
236 pFFT[i] = RealIn[i];
237
238 // Perform the FFT
239 RealFFTf(pFFT.get(), hFFT.get());
240
241 // Copy the data into the real and imaginary outputs
242 for (size_t i = 1; i<(NumSamples / 2); i++) {
243 RealOut[i]=pFFT[hFFT->BitReversed[i] ];
244 ImagOut[i]=pFFT[hFFT->BitReversed[i]+1];
245 }
246 // Handle the (real-only) DC and Fs/2 bins
247 RealOut[0] = pFFT[0];
248 RealOut[NumSamples / 2] = pFFT[1];
249 ImagOut[0] = ImagOut[NumSamples / 2] = 0;
250 // Fill in the upper half using symmetry properties
251 for(size_t i = NumSamples / 2 + 1; i < NumSamples; i++) {
252 RealOut[i] = RealOut[NumSamples-i];
253 ImagOut[i] = -ImagOut[NumSamples-i];
254 }
255 }
256
257 /*
258 * InverseRealFFT
259 *
260 * This function computes the inverse of RealFFT, above.
261 * The RealIn and ImagIn is assumed to be conjugate-symmetric
262 * and as a result the output is purely real.
263 * Only the first half of RealIn and ImagIn are used due to this
264 * symmetry assumption.
265 *
266 * This is merely a wrapper of InverseRealFFTf() from RealFFTf.h.
267 */
InverseRealFFT(size_t NumSamples,const float * RealIn,const float * ImagIn,float * RealOut)268 void InverseRealFFT(size_t NumSamples, const float *RealIn, const float *ImagIn,
269 float *RealOut)
270 {
271 auto hFFT = GetFFT(NumSamples);
272 Floats pFFT{ NumSamples };
273 // Copy the data into the processing buffer
274 for (size_t i = 0; i < (NumSamples / 2); i++)
275 pFFT[2*i ] = RealIn[i];
276 if(ImagIn == NULL) {
277 for (size_t i = 0; i < (NumSamples / 2); i++)
278 pFFT[2*i+1] = 0;
279 } else {
280 for (size_t i = 0; i < (NumSamples / 2); i++)
281 pFFT[2*i+1] = ImagIn[i];
282 }
283 // Put the fs/2 component in the imaginary part of the DC bin
284 pFFT[1] = RealIn[NumSamples / 2];
285
286 // Perform the FFT
287 InverseRealFFTf(pFFT.get(), hFFT.get());
288
289 // Copy the data to the (purely real) output buffer
290 ReorderToTime(hFFT.get(), pFFT.get(), RealOut);
291 }
292
293 /*
294 * PowerSpectrum
295 *
296 * This function uses RealFFTf() from RealFFTf.h to perform the real
297 * FFT computation, and then squares the real and imaginary part of
298 * each coefficient, extracting the power and throwing away the phase.
299 *
300 * For speed, it does not call RealFFT, but duplicates some
301 * of its code.
302 */
303
PowerSpectrum(size_t NumSamples,const float * In,float * Out)304 void PowerSpectrum(size_t NumSamples, const float *In, float *Out)
305 {
306 auto hFFT = GetFFT(NumSamples);
307 Floats pFFT{ NumSamples };
308 // Copy the data into the processing buffer
309 for (size_t i = 0; i<NumSamples; i++)
310 pFFT[i] = In[i];
311
312 // Perform the FFT
313 RealFFTf(pFFT.get(), hFFT.get());
314
315 // Copy the data into the real and imaginary outputs
316 for (size_t i = 1; i<NumSamples / 2; i++) {
317 Out[i]= (pFFT[hFFT->BitReversed[i] ]*pFFT[hFFT->BitReversed[i] ])
318 + (pFFT[hFFT->BitReversed[i]+1]*pFFT[hFFT->BitReversed[i]+1]);
319 }
320 // Handle the (real-only) DC and Fs/2 bins
321 Out[0] = pFFT[0]*pFFT[0];
322 Out[NumSamples / 2] = pFFT[1]*pFFT[1];
323 }
324
325 /*
326 * Windowing Functions
327 */
328
NumWindowFuncs()329 int NumWindowFuncs()
330 {
331 return eWinFuncCount;
332 }
333
WindowFuncName(int whichFunction)334 const TranslatableString WindowFuncName(int whichFunction)
335 {
336 switch (whichFunction) {
337 default:
338 case eWinFuncRectangular:
339 return XO("Rectangular");
340 case eWinFuncBartlett:
341 /* i18n-hint a proper name */
342 return XO("Bartlett");
343 case eWinFuncHamming:
344 /* i18n-hint a proper name */
345 return XO("Hamming");
346 case eWinFuncHann:
347 /* i18n-hint a proper name */
348 return XO("Hann");
349 case eWinFuncBlackman:
350 /* i18n-hint a proper name */
351 return XO("Blackman");
352 case eWinFuncBlackmanHarris:
353 /* i18n-hint two proper names */
354 return XO("Blackman-Harris");
355 case eWinFuncWelch:
356 /* i18n-hint a proper name */
357 return XO("Welch");
358 case eWinFuncGaussian25:
359 /* i18n-hint a mathematical function named for C. F. Gauss */
360 return XO("Gaussian(a=2.5)");
361 case eWinFuncGaussian35:
362 /* i18n-hint a mathematical function named for C. F. Gauss */
363 return XO("Gaussian(a=3.5)");
364 case eWinFuncGaussian45:
365 /* i18n-hint a mathematical function named for C. F. Gauss */
366 return XO("Gaussian(a=4.5)");
367 }
368 }
369
NewWindowFunc(int whichFunction,size_t NumSamplesIn,bool extraSample,float * in)370 void NewWindowFunc(int whichFunction, size_t NumSamplesIn, bool extraSample, float *in)
371 {
372 int NumSamples = (int)NumSamplesIn;
373 if (extraSample) {
374 wxASSERT(NumSamples > 0);
375 --NumSamples;
376 }
377 wxASSERT(NumSamples > 0);
378
379 switch (whichFunction) {
380 default:
381 wxFprintf(stderr, "FFT::WindowFunc - Invalid window function: %d\n", whichFunction);
382 break;
383 case eWinFuncRectangular:
384 // Multiply all by 1.0f -- do nothing
385 break;
386
387 case eWinFuncBartlett:
388 {
389 // Bartlett (triangular) window
390 const int nPairs = (NumSamples - 1) / 2; // whether even or odd NumSamples, this is correct
391 const float denom = NumSamples / 2.0f;
392 in[0] = 0.0f;
393 for (int ii = 1;
394 ii <= nPairs; // Yes, <=
395 ++ii) {
396 const float value = ii / denom;
397 in[ii] *= value;
398 in[NumSamples - ii] *= value;
399 }
400 // When NumSamples is even, in[half] should be multiplied by 1.0, so unchanged
401 // When odd, the value of 1.0 is not reached
402 }
403 break;
404 case eWinFuncHamming:
405 {
406 // Hamming
407 const double multiplier = 2 * M_PI / NumSamples;
408 static const double coeff0 = 0.54, coeff1 = -0.46;
409 for (int ii = 0; ii < NumSamples; ++ii)
410 in[ii] *= coeff0 + coeff1 * cos(ii * multiplier);
411 }
412 break;
413 case eWinFuncHann:
414 {
415 // Hann
416 const double multiplier = 2 * M_PI / NumSamples;
417 static const double coeff0 = 0.5, coeff1 = -0.5;
418 for (int ii = 0; ii < NumSamples; ++ii)
419 in[ii] *= coeff0 + coeff1 * cos(ii * multiplier);
420 }
421 break;
422 case eWinFuncBlackman:
423 {
424 // Blackman
425 const double multiplier = 2 * M_PI / NumSamples;
426 const double multiplier2 = 2 * multiplier;
427 static const double coeff0 = 0.42, coeff1 = -0.5, coeff2 = 0.08;
428 for (int ii = 0; ii < NumSamples; ++ii)
429 in[ii] *= coeff0 + coeff1 * cos(ii * multiplier) + coeff2 * cos(ii * multiplier2);
430 }
431 break;
432 case eWinFuncBlackmanHarris:
433 {
434 // Blackman-Harris
435 const double multiplier = 2 * M_PI / NumSamples;
436 const double multiplier2 = 2 * multiplier;
437 const double multiplier3 = 3 * multiplier;
438 static const double coeff0 = 0.35875, coeff1 = -0.48829, coeff2 = 0.14128, coeff3 = -0.01168;
439 for (int ii = 0; ii < NumSamples; ++ii)
440 in[ii] *= coeff0 + coeff1 * cos(ii * multiplier) + coeff2 * cos(ii * multiplier2) + coeff3 * cos(ii * multiplier3);
441 }
442 break;
443 case eWinFuncWelch:
444 {
445 // Welch
446 const float N = NumSamples;
447 for (int ii = 0; ii < NumSamples; ++ii) {
448 const float iOverN = ii / N;
449 in[ii] *= 4 * iOverN * (1 - iOverN);
450 }
451 }
452 break;
453 case eWinFuncGaussian25:
454 {
455 // Gaussian (a=2.5)
456 // Precalculate some values, and simplify the fmla to try and reduce overhead
457 static const double A = -2 * 2.5*2.5;
458 const float N = NumSamples;
459 for (int ii = 0; ii < NumSamples; ++ii) {
460 const float iOverN = ii / N;
461 // full
462 // in[ii] *= exp(-0.5*(A*((ii-NumSamples/2)/NumSamples/2))*(A*((ii-NumSamples/2)/NumSamples/2)));
463 // reduced
464 in[ii] *= exp(A * (0.25 + (iOverN * iOverN) - iOverN));
465 }
466 }
467 break;
468 case eWinFuncGaussian35:
469 {
470 // Gaussian (a=3.5)
471 static const double A = -2 * 3.5*3.5;
472 const float N = NumSamples;
473 for (int ii = 0; ii < NumSamples; ++ii) {
474 const float iOverN = ii / N;
475 in[ii] *= exp(A * (0.25 + (iOverN * iOverN) - iOverN));
476 }
477 }
478 break;
479 case eWinFuncGaussian45:
480 {
481 // Gaussian (a=4.5)
482 static const double A = -2 * 4.5*4.5;
483 const float N = NumSamples;
484 for (int ii = 0; ii < NumSamples; ++ii) {
485 const float iOverN = ii / N;
486 in[ii] *= exp(A * (0.25 + (iOverN * iOverN) - iOverN));
487 }
488 }
489 break;
490 }
491
492 if (extraSample && whichFunction != eWinFuncRectangular) {
493 double value = 0.0;
494 switch (whichFunction) {
495 case eWinFuncHamming:
496 value = 0.08;
497 break;
498 case eWinFuncGaussian25:
499 value = exp(-2 * 2.5 * 2.5 * 0.25);
500 break;
501 case eWinFuncGaussian35:
502 value = exp(-2 * 3.5 * 3.5 * 0.25);
503 break;
504 case eWinFuncGaussian45:
505 value = exp(-2 * 4.5 * 4.5 * 0.25);
506 break;
507 default:
508 break;
509 }
510 in[NumSamples] *= value;
511 }
512 }
513
514 // See cautions in FFT.h !
WindowFunc(int whichFunction,size_t NumSamples,float * in)515 void WindowFunc(int whichFunction, size_t NumSamples, float *in)
516 {
517 bool extraSample = false;
518 switch (whichFunction)
519 {
520 case eWinFuncHamming:
521 case eWinFuncHann:
522 case eWinFuncBlackman:
523 case eWinFuncBlackmanHarris:
524 extraSample = true;
525 break;
526 default:
527 break;
528 case eWinFuncBartlett:
529 // PRL: Do nothing here either
530 // But I want to comment that the old function did this case
531 // wrong in the second half of the array, in case NumSamples was odd
532 // but I think that never happened, so I am not bothering to preserve that
533 break;
534 }
535 NewWindowFunc(whichFunction, NumSamples, extraSample, in);
536 }
537
DerivativeOfWindowFunc(int whichFunction,size_t NumSamples,bool extraSample,float * in)538 void DerivativeOfWindowFunc(int whichFunction, size_t NumSamples, bool extraSample, float *in)
539 {
540 if (eWinFuncRectangular == whichFunction)
541 {
542 // Rectangular
543 // There are deltas at the ends
544 wxASSERT(NumSamples > 0);
545 --NumSamples;
546 // in[0] *= 1.0f;
547 for (int ii = 1; ii < (int)NumSamples; ++ii)
548 in[ii] = 0.0f;
549 in[NumSamples] *= -1.0f;
550 return;
551 }
552
553 if (extraSample) {
554 wxASSERT(NumSamples > 0);
555 --NumSamples;
556 }
557
558 wxASSERT(NumSamples > 0);
559
560 double A;
561 switch (whichFunction) {
562 case eWinFuncBartlett:
563 {
564 // Bartlett (triangular) window
565 // There are discontinuities in the derivative at the ends, and maybe at the midpoint
566 const int nPairs = (NumSamples - 1) / 2; // whether even or odd NumSamples, this is correct
567 const float value = 2.0f / NumSamples;
568 in[0] *=
569 // Average the two limiting values of discontinuous derivative
570 value / 2.0f;
571 for (int ii = 1;
572 ii <= nPairs; // Yes, <=
573 ++ii) {
574 in[ii] *= value;
575 in[NumSamples - ii] *= -value;
576 }
577 if (NumSamples % 2 == 0)
578 // Average the two limiting values of discontinuous derivative
579 in[NumSamples / 2] = 0.0f;
580 if (extraSample)
581 in[NumSamples] *=
582 // Average the two limiting values of discontinuous derivative
583 -value / 2.0f;
584 else
585 // Halve the multiplier previously applied
586 // Average the two limiting values of discontinuous derivative
587 in[NumSamples - 1] *= 0.5f;
588 }
589 break;
590 case eWinFuncHamming:
591 {
592 // Hamming
593 // There are deltas at the ends
594 const double multiplier = 2 * M_PI / NumSamples;
595 static const double coeff0 = 0.54, coeff1 = -0.46 * multiplier;
596 // TODO This code should be more explicit about the precision it intends.
597 // For now we get C4305 warnings, truncation from 'const double' to 'float'
598 in[0] *= coeff0;
599 if (!extraSample)
600 --NumSamples;
601 for (int ii = 0; ii < (int)NumSamples; ++ii)
602 in[ii] *= - coeff1 * sin(ii * multiplier);
603 if (extraSample)
604 in[NumSamples] *= - coeff0;
605 else
606 // slightly different
607 in[NumSamples] *= - coeff0 - coeff1 * sin(NumSamples * multiplier);
608 }
609 break;
610 case eWinFuncHann:
611 {
612 // Hann
613 const double multiplier = 2 * M_PI / NumSamples;
614 const double coeff1 = -0.5 * multiplier;
615 for (int ii = 0; ii < (int)NumSamples; ++ii)
616 in[ii] *= - coeff1 * sin(ii * multiplier);
617 if (extraSample)
618 in[NumSamples] = 0.0f;
619 }
620 break;
621 case eWinFuncBlackman:
622 {
623 // Blackman
624 const double multiplier = 2 * M_PI / NumSamples;
625 const double multiplier2 = 2 * multiplier;
626 const double coeff1 = -0.5 * multiplier, coeff2 = 0.08 * multiplier2;
627 for (int ii = 0; ii < (int)NumSamples; ++ii)
628 in[ii] *= - coeff1 * sin(ii * multiplier) - coeff2 * sin(ii * multiplier2);
629 if (extraSample)
630 in[NumSamples] = 0.0f;
631 }
632 break;
633 case eWinFuncBlackmanHarris:
634 {
635 // Blackman-Harris
636 const double multiplier = 2 * M_PI / NumSamples;
637 const double multiplier2 = 2 * multiplier;
638 const double multiplier3 = 3 * multiplier;
639 const double coeff1 = -0.48829 * multiplier,
640 coeff2 = 0.14128 * multiplier2, coeff3 = -0.01168 * multiplier3;
641 for (int ii = 0; ii < (int)NumSamples; ++ii)
642 in[ii] *= - coeff1 * sin(ii * multiplier) - coeff2 * sin(ii * multiplier2) - coeff3 * sin(ii * multiplier3);
643 if (extraSample)
644 in[NumSamples] = 0.0f;
645 }
646 break;
647 case eWinFuncWelch:
648 {
649 // Welch
650 const float N = NumSamples;
651 const float NN = NumSamples * NumSamples;
652 for (int ii = 0; ii < (int)NumSamples; ++ii) {
653 in[ii] *= 4 * (N - ii - ii) / NN;
654 }
655 if (extraSample)
656 in[NumSamples] = 0.0f;
657 // Average the two limiting values of discontinuous derivative
658 in[0] /= 2.0f;
659 in[NumSamples - 1] /= 2.0f;
660 }
661 break;
662 case eWinFuncGaussian25:
663 // Gaussian (a=2.5)
664 A = -2 * 2.5*2.5;
665 goto Gaussian;
666 case eWinFuncGaussian35:
667 // Gaussian (a=3.5)
668 A = -2 * 3.5*3.5;
669 goto Gaussian;
670 case eWinFuncGaussian45:
671 // Gaussian (a=4.5)
672 A = -2 * 4.5*4.5;
673 goto Gaussian;
674 Gaussian:
675 {
676 // Gaussian (a=2.5)
677 // There are deltas at the ends
678 const float invN = 1.0f / NumSamples;
679 const float invNN = invN * invN;
680 // Simplify formula from the loop for ii == 0, add term for the delta
681 in[0] *= exp(A * 0.25) * (1 - invN);
682 if (!extraSample)
683 --NumSamples;
684 for (int ii = 1; ii < (int)NumSamples; ++ii) {
685 const float iOverN = ii * invN;
686 in[ii] *= exp(A * (0.25 + (iOverN * iOverN) - iOverN)) * (2 * ii * invNN - invN);
687 }
688 if (extraSample)
689 in[NumSamples] *= exp(A * 0.25) * (invN - 1);
690 else {
691 // Slightly different
692 const float iOverN = NumSamples * invN;
693 in[NumSamples] *= exp(A * (0.25 + (iOverN * iOverN) - iOverN)) * (2 * NumSamples * invNN - invN - 1);
694 }
695 }
696 break;
697 default:
698 wxFprintf(stderr, "FFT::DerivativeOfWindowFunc - Invalid window function: %d\n", whichFunction);
699 }
700 }
701