OpenShot Audio Library | OpenShotAudio  0.3.1
juce_FFT.cpp
1 /*
2  ==============================================================================
3 
4  This file is part of the JUCE library.
5  Copyright (c) 2017 - ROLI Ltd.
6 
7  JUCE is an open source library subject to commercial or open-source
8  licensing.
9 
10  By using JUCE, you agree to the terms of both the JUCE 5 End-User License
11  Agreement and JUCE 5 Privacy Policy (both updated and effective as of the
12  27th April 2017).
13 
14  End User License Agreement: www.juce.com/juce-5-licence
15  Privacy Policy: www.juce.com/juce-5-privacy-policy
16 
17  Or: You may also use this code under the terms of the GPL v3 (see
18  www.gnu.org/licenses).
19 
20  JUCE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY, AND ALL WARRANTIES, WHETHER
21  EXPRESSED OR IMPLIED, INCLUDING MERCHANTABILITY AND FITNESS FOR PURPOSE, ARE
22  DISCLAIMED.
23 
24  ==============================================================================
25 */
26 
27 namespace juce
28 {
29 namespace dsp
30 {
31 
32 struct FFT::Instance
33 {
34  virtual ~Instance() {}
35  virtual void perform (const Complex<float>* input, Complex<float>* output, bool inverse) const noexcept = 0;
36  virtual void performRealOnlyForwardTransform (float*, bool) const noexcept = 0;
37  virtual void performRealOnlyInverseTransform (float*) const noexcept = 0;
38 };
39 
40 struct FFT::Engine
41 {
42  Engine (int priorityToUse) : enginePriority (priorityToUse)
43  {
44  auto& list = getEngines();
45  list.add (this);
46  std::sort (list.begin(), list.end(), [] (Engine* a, Engine* b) { return b->enginePriority < a->enginePriority; });
47  }
48 
49  virtual ~Engine() {}
50 
51  virtual FFT::Instance* create (int order) const = 0;
52 
53  //==============================================================================
54  static FFT::Instance* createBestEngineForPlatform (int order)
55  {
56  for (auto* engine : getEngines())
57  if (auto* instance = engine->create (order))
58  return instance;
59 
60  jassertfalse; // This should never happen as the fallback engine should always work!
61  return nullptr;
62  }
63 
64 private:
65  static Array<Engine*>& getEngines()
66  {
67  static Array<Engine*> engines;
68  return engines;
69  }
70 
71  int enginePriority; // used so that faster engines have priority over slower ones
72 };
73 
74 template <typename InstanceToUse>
75 struct FFT::EngineImpl : public FFT::Engine
76 {
77  EngineImpl() : FFT::Engine (InstanceToUse::priority) {}
78  FFT::Instance* create (int order) const override { return InstanceToUse::create (order); }
79 };
80 
81 //==============================================================================
82 //==============================================================================
83 struct FFTFallback : public FFT::Instance
84 {
85  // this should have the least priority of all engines
86  static constexpr int priority = -1;
87 
88  static FFTFallback* create (int order)
89  {
90  return new FFTFallback (order);
91  }
92 
93  FFTFallback (int order)
94  {
95  configForward.reset (new FFTConfig (1 << order, false));
96  configInverse.reset (new FFTConfig (1 << order, true));
97 
98  size = 1 << order;
99  }
100 
101  void perform (const Complex<float>* input, Complex<float>* output, bool inverse) const noexcept override
102  {
103  if (size == 1)
104  {
105  *output = *input;
106  return;
107  }
108 
109  const SpinLock::ScopedLockType sl(processLock);
110 
111  jassert (configForward != nullptr);
112 
113  if (inverse)
114  {
115  configInverse->perform (input, output);
116 
117  const float scaleFactor = 1.0f / size;
118 
119  for (int i = 0; i < size; ++i)
120  output[i] *= scaleFactor;
121  }
122  else
123  {
124  configForward->perform (input, output);
125  }
126  }
127 
128  const size_t maxFFTScratchSpaceToAlloca = 256 * 1024;
129 
130  void performRealOnlyForwardTransform (float* d, bool) const noexcept override
131  {
132  if (size == 1)
133  return;
134 
135  const size_t scratchSize = 16 + (size_t) size * sizeof (Complex<float>);
136 
137  if (scratchSize < maxFFTScratchSpaceToAlloca)
138  {
139  performRealOnlyForwardTransform (static_cast<Complex<float>*> (alloca (scratchSize)), d);
140  }
141  else
142  {
143  HeapBlock<char> heapSpace (scratchSize);
144  performRealOnlyForwardTransform (reinterpret_cast<Complex<float>*> (heapSpace.getData()), d);
145  }
146  }
147 
148  void performRealOnlyInverseTransform (float* d) const noexcept override
149  {
150  if (size == 1)
151  return;
152 
153  const size_t scratchSize = 16 + (size_t) size * sizeof (Complex<float>);
154 
155  if (scratchSize < maxFFTScratchSpaceToAlloca)
156  {
157  performRealOnlyInverseTransform (static_cast<Complex<float>*> (alloca (scratchSize)), d);
158  }
159  else
160  {
161  HeapBlock<char> heapSpace (scratchSize);
162  performRealOnlyInverseTransform (reinterpret_cast<Complex<float>*> (heapSpace.getData()), d);
163  }
164  }
165 
166  void performRealOnlyForwardTransform (Complex<float>* scratch, float* d) const noexcept
167  {
168  for (int i = 0; i < size; ++i)
169  scratch[i] = { d[i], 0 };
170 
171  perform (scratch, reinterpret_cast<Complex<float>*> (d), false);
172  }
173 
174  void performRealOnlyInverseTransform (Complex<float>* scratch, float* d) const noexcept
175  {
176  auto* input = reinterpret_cast<Complex<float>*> (d);
177 
178  for (auto i = size >> 1; i < size; ++i)
179  input[i] = std::conj (input[size - i]);
180 
181  perform (input, scratch, true);
182 
183  for (int i = 0; i < size; ++i)
184  {
185  d[i] = scratch[i].real();
186  d[i + size] = scratch[i].imag();
187  }
188  }
189 
190  //==============================================================================
191  struct FFTConfig
192  {
193  FFTConfig (int sizeOfFFT, bool isInverse)
194  : fftSize (sizeOfFFT), inverse (isInverse), twiddleTable ((size_t) sizeOfFFT)
195  {
196  auto inverseFactor = (inverse ? 2.0 : -2.0) * MathConstants<double>::pi / (double) fftSize;
197 
198  if (fftSize <= 4)
199  {
200  for (int i = 0; i < fftSize; ++i)
201  {
202  auto phase = i * inverseFactor;
203 
204  twiddleTable[i] = { (float) std::cos (phase),
205  (float) std::sin (phase) };
206  }
207  }
208  else
209  {
210  for (int i = 0; i < fftSize / 4; ++i)
211  {
212  auto phase = i * inverseFactor;
213 
214  twiddleTable[i] = { (float) std::cos (phase),
215  (float) std::sin (phase) };
216  }
217 
218  for (int i = fftSize / 4; i < fftSize / 2; ++i)
219  {
220  auto other = twiddleTable[i - fftSize / 4];
221 
222  twiddleTable[i] = { inverse ? -other.imag() : other.imag(),
223  inverse ? other.real() : -other.real() };
224  }
225 
226  twiddleTable[fftSize / 2].real (-1.0f);
227  twiddleTable[fftSize / 2].imag (0.0f);
228 
229  for (int i = fftSize / 2; i < fftSize; ++i)
230  {
231  auto index = fftSize / 2 - (i - fftSize / 2);
232  twiddleTable[i] = conj(twiddleTable[index]);
233  }
234  }
235 
236  auto root = (int) std::sqrt ((double) fftSize);
237  int divisor = 4, n = fftSize;
238 
239  for (int i = 0; i < numElementsInArray (factors); ++i)
240  {
241  while ((n % divisor) != 0)
242  {
243  if (divisor == 2) divisor = 3;
244  else if (divisor == 4) divisor = 2;
245  else divisor += 2;
246 
247  if (divisor > root)
248  divisor = n;
249  }
250 
251  n /= divisor;
252 
253  jassert (divisor == 1 || divisor == 2 || divisor == 4);
254  factors[i].radix = divisor;
255  factors[i].length = n;
256  }
257  }
258 
259  void perform (const Complex<float>* input, Complex<float>* output) const noexcept
260  {
261  perform (input, output, 1, 1, factors);
262  }
263 
264  const int fftSize;
265  const bool inverse;
266 
267  struct Factor { int radix, length; };
268  Factor factors[32];
269  HeapBlock<Complex<float>> twiddleTable;
270 
271  void perform (const Complex<float>* input, Complex<float>* output, int stride, int strideIn, const Factor* facs) const noexcept
272  {
273  auto factor = *facs++;
274  auto* originalOutput = output;
275  auto* outputEnd = output + factor.radix * factor.length;
276 
277  if (stride == 1 && factor.radix <= 5)
278  {
279  for (int i = 0; i < factor.radix; ++i)
280  perform (input + stride * strideIn * i, output + i * factor.length, stride * factor.radix, strideIn, facs);
281 
282  butterfly (factor, output, stride);
283  return;
284  }
285 
286  if (factor.length == 1)
287  {
288  do
289  {
290  *output++ = *input;
291  input += stride * strideIn;
292  }
293  while (output < outputEnd);
294  }
295  else
296  {
297  do
298  {
299  perform (input, output, stride * factor.radix, strideIn, facs);
300  input += stride * strideIn;
301  output += factor.length;
302  }
303  while (output < outputEnd);
304  }
305 
306  butterfly (factor, originalOutput, stride);
307  }
308 
309  void butterfly (const Factor factor, Complex<float>* data, int stride) const noexcept
310  {
311  switch (factor.radix)
312  {
313  case 1: break;
314  case 2: butterfly2 (data, stride, factor.length); return;
315  case 4: butterfly4 (data, stride, factor.length); return;
316  default: jassertfalse; break;
317  }
318 
319  auto* scratch = static_cast<Complex<float>*> (alloca ((size_t) factor.radix * sizeof (Complex<float>)));
320 
321  for (int i = 0; i < factor.length; ++i)
322  {
323  for (int k = i, q1 = 0; q1 < factor.radix; ++q1)
324  {
325  scratch[q1] = data[k];
326  k += factor.length;
327  }
328 
329  for (int k = i, q1 = 0; q1 < factor.radix; ++q1)
330  {
331  int twiddleIndex = 0;
332  data[k] = scratch[0];
333 
334  for (int q = 1; q < factor.radix; ++q)
335  {
336  twiddleIndex += stride * k;
337 
338  if (twiddleIndex >= fftSize)
339  twiddleIndex -= fftSize;
340 
341  data[k] += scratch[q] * twiddleTable[twiddleIndex];
342  }
343 
344  k += factor.length;
345  }
346  }
347  }
348 
349  void butterfly2 (Complex<float>* data, const int stride, const int length) const noexcept
350  {
351  auto* dataEnd = data + length;
352  auto* tw = twiddleTable.getData();
353 
354  for (int i = length; --i >= 0;)
355  {
356  auto s = *dataEnd;
357  s *= (*tw);
358  tw += stride;
359  *dataEnd++ = *data - s;
360  *data++ += s;
361  }
362  }
363 
364  void butterfly4 (Complex<float>* data, const int stride, const int length) const noexcept
365  {
366  auto lengthX2 = length * 2;
367  auto lengthX3 = length * 3;
368 
369  auto strideX2 = stride * 2;
370  auto strideX3 = stride * 3;
371 
372  auto* twiddle1 = twiddleTable.getData();
373  auto* twiddle2 = twiddle1;
374  auto* twiddle3 = twiddle1;
375 
376  for (int i = length; --i >= 0;)
377  {
378  auto s0 = data[length] * *twiddle1;
379  auto s1 = data[lengthX2] * *twiddle2;
380  auto s2 = data[lengthX3] * *twiddle3;
381  auto s3 = s0; s3 += s2;
382  auto s4 = s0; s4 -= s2;
383  auto s5 = *data; s5 -= s1;
384 
385  *data += s1;
386  data[lengthX2] = *data;
387  data[lengthX2] -= s3;
388  twiddle1 += stride;
389  twiddle2 += strideX2;
390  twiddle3 += strideX3;
391  *data += s3;
392 
393  if (inverse)
394  {
395  data[length] = { s5.real() - s4.imag(),
396  s5.imag() + s4.real() };
397 
398  data[lengthX3] = { s5.real() + s4.imag(),
399  s5.imag() - s4.real() };
400  }
401  else
402  {
403  data[length] = { s5.real() + s4.imag(),
404  s5.imag() - s4.real() };
405 
406  data[lengthX3] = { s5.real() - s4.imag(),
407  s5.imag() + s4.real() };
408  }
409 
410  ++data;
411  }
412  }
413 
414  JUCE_DECLARE_NON_COPYABLE_WITH_LEAK_DETECTOR (FFTConfig)
415  };
416 
417  //==============================================================================
418  SpinLock processLock;
419  std::unique_ptr<FFTConfig> configForward, configInverse;
420  int size;
421 };
422 
423 FFT::EngineImpl<FFTFallback> fftFallback;
424 
425 //==============================================================================
426 //==============================================================================
427 #if (JUCE_MAC || JUCE_IOS) && JUCE_USE_VDSP_FRAMEWORK
428 struct AppleFFT : public FFT::Instance
429 {
430  static constexpr int priority = 5;
431 
432  static AppleFFT* create (int order)
433  {
434  return new AppleFFT (order);
435  }
436 
437  AppleFFT (int orderToUse)
438  : order (static_cast<vDSP_Length> (orderToUse)),
439  fftSetup (vDSP_create_fftsetup (order, 2)),
440  forwardNormalisation (0.5f),
441  inverseNormalisation (1.0f / static_cast<float> (1 << order))
442  {}
443 
444  ~AppleFFT() override
445  {
446  if (fftSetup != nullptr)
447  {
448  vDSP_destroy_fftsetup (fftSetup);
449  fftSetup = nullptr;
450  }
451  }
452 
453  void perform (const Complex<float>* input, Complex<float>* output, bool inverse) const noexcept override
454  {
455  auto size = (1 << order);
456 
457  DSPSplitComplex splitInput (toSplitComplex (const_cast<Complex<float>*> (input)));
458  DSPSplitComplex splitOutput (toSplitComplex (output));
459 
460  vDSP_fft_zop (fftSetup, &splitInput, 2, &splitOutput, 2,
461  order, inverse ? kFFTDirection_Inverse : kFFTDirection_Forward);
462 
463  float factor = (inverse ? inverseNormalisation : forwardNormalisation * 2.0f);
464  vDSP_vsmul ((float*) output, 1, &factor, (float*) output, 1, static_cast<size_t> (size << 1));
465  }
466 
467  void performRealOnlyForwardTransform (float* inoutData, bool ignoreNegativeFreqs) const noexcept override
468  {
469  auto size = (1 << order);
470  auto* inout = reinterpret_cast<Complex<float>*> (inoutData);
471  auto splitInOut (toSplitComplex (inout));
472 
473  inoutData[size] = 0.0f;
474  vDSP_fft_zrip (fftSetup, &splitInOut, 2, order, kFFTDirection_Forward);
475  vDSP_vsmul (inoutData, 1, &forwardNormalisation, inoutData, 1, static_cast<size_t> (size << 1));
476 
477  mirrorResult (inout, ignoreNegativeFreqs);
478  }
479 
480  void performRealOnlyInverseTransform (float* inoutData) const noexcept override
481  {
482  auto* inout = reinterpret_cast<Complex<float>*> (inoutData);
483  auto size = (1 << order);
484  auto splitInOut (toSplitComplex (inout));
485 
486  // Imaginary part of nyquist and DC frequencies are always zero
487  // so Apple uses the imaginary part of the DC frequency to store
488  // the real part of the nyquist frequency
489  if (size != 1)
490  inout[0] = Complex<float> (inout[0].real(), inout[size >> 1].real());
491 
492  vDSP_fft_zrip (fftSetup, &splitInOut, 2, order, kFFTDirection_Inverse);
493  vDSP_vsmul (inoutData, 1, &inverseNormalisation, inoutData, 1, static_cast<size_t> (size << 1));
494  vDSP_vclr (inoutData + size, 1, static_cast<size_t> (size));
495  }
496 
497 private:
498  //==============================================================================
499  void mirrorResult (Complex<float>* out, bool ignoreNegativeFreqs) const noexcept
500  {
501  auto size = (1 << order);
502  auto i = size >> 1;
503 
504  // Imaginary part of nyquist and DC frequencies are always zero
505  // so Apple uses the imaginary part of the DC frequency to store
506  // the real part of the nyquist frequency
507  out[i++] = { out[0].imag(), 0.0 };
508  out[0] = { out[0].real(), 0.0 };
509 
510  if (! ignoreNegativeFreqs)
511  for (; i < size; ++i)
512  out[i] = std::conj (out[size - i]);
513  }
514 
515  static DSPSplitComplex toSplitComplex (Complex<float>* data) noexcept
516  {
517  // this assumes that Complex interleaves real and imaginary parts
518  // and is tightly packed.
519  return { reinterpret_cast<float*> (data),
520  reinterpret_cast<float*> (data) + 1};
521  }
522 
523  //==============================================================================
524  vDSP_Length order;
525  FFTSetup fftSetup;
526  float forwardNormalisation, inverseNormalisation;
527 };
528 
529 FFT::EngineImpl<AppleFFT> appleFFT;
530 #endif
531 
532 //==============================================================================
533 //==============================================================================
534 #if JUCE_DSP_USE_SHARED_FFTW || JUCE_DSP_USE_STATIC_FFTW
535 
536 #if JUCE_DSP_USE_STATIC_FFTW
537 extern "C"
538 {
539  void* fftwf_plan_dft_1d (int, void*, void*, int, int);
540  void* fftwf_plan_dft_r2c_1d (int, void*, void*, int);
541  void* fftwf_plan_dft_c2r_1d (int, void*, void*, int);
542  void fftwf_destroy_plan (void*);
543  void fftwf_execute_dft (void*, void*, void*);
544  void fftwf_execute_dft_r2c (void*, void*, void*);
545  void fftwf_execute_dft_c2r (void*, void*, void*);
546 }
547 #endif
548 
549 struct FFTWImpl : public FFT::Instance
550 {
551  #if JUCE_DSP_USE_STATIC_FFTW
552  // if the JUCE developer has gone through the hassle of statically
553  // linking in fftw, they probably want to use it
554  static constexpr int priority = 10;
555  #else
556  static constexpr int priority = 3;
557  #endif
558 
559  struct FFTWPlan;
560  using FFTWPlanRef = FFTWPlan*;
561 
562  enum
563  {
564  measure = 0,
565  unaligned = (1 << 1),
566  estimate = (1 << 6)
567  };
568 
569  struct Symbols
570  {
571  FFTWPlanRef (*plan_dft_fftw) (unsigned, Complex<float>*, Complex<float>*, int, unsigned);
572  FFTWPlanRef (*plan_r2c_fftw) (unsigned, float*, Complex<float>*, unsigned);
573  FFTWPlanRef (*plan_c2r_fftw) (unsigned, Complex<float>*, float*, unsigned);
574  void (*destroy_fftw) (FFTWPlanRef);
575 
576  void (*execute_dft_fftw) (FFTWPlanRef, const Complex<float>*, Complex<float>*);
577  void (*execute_r2c_fftw) (FFTWPlanRef, float*, Complex<float>*);
578  void (*execute_c2r_fftw) (FFTWPlanRef, Complex<float>*, float*);
579 
580  #if JUCE_DSP_USE_STATIC_FFTW
581  template <typename FuncPtr, typename ActualSymbolType>
582  static bool symbol (FuncPtr& dst, ActualSymbolType sym)
583  {
584  dst = reinterpret_cast<FuncPtr> (sym);
585  return true;
586  }
587  #else
588  template <typename FuncPtr>
589  static bool symbol (DynamicLibrary& lib, FuncPtr& dst, const char* name)
590  {
591  dst = reinterpret_cast<FuncPtr> (lib.getFunction (name));
592  return (dst != nullptr);
593  }
594  #endif
595  };
596 
597  static FFTWImpl* create (int order)
598  {
599  DynamicLibrary lib;
600 
601  #if ! JUCE_DSP_USE_STATIC_FFTW
602  #if JUCE_MAC
603  auto libName = "libfftw3f.dylib";
604  #elif JUCE_WINDOWS
605  auto libName = "libfftw3f.dll";
606  #else
607  auto libName = "libfftw3f.so";
608  #endif
609 
610  if (lib.open (libName))
611  #endif
612  {
613  Symbols symbols;
614 
615  #if JUCE_DSP_USE_STATIC_FFTW
616  if (! Symbols::symbol (symbols.plan_dft_fftw, fftwf_plan_dft_1d)) return nullptr;
617  if (! Symbols::symbol (symbols.plan_r2c_fftw, fftwf_plan_dft_r2c_1d)) return nullptr;
618  if (! Symbols::symbol (symbols.plan_c2r_fftw, fftwf_plan_dft_c2r_1d)) return nullptr;
619  if (! Symbols::symbol (symbols.destroy_fftw, fftwf_destroy_plan)) return nullptr;
620 
621  if (! Symbols::symbol (symbols.execute_dft_fftw, fftwf_execute_dft)) return nullptr;
622  if (! Symbols::symbol (symbols.execute_r2c_fftw, fftwf_execute_dft_r2c)) return nullptr;
623  if (! Symbols::symbol (symbols.execute_c2r_fftw, fftwf_execute_dft_c2r)) return nullptr;
624  #else
625  if (! Symbols::symbol (lib, symbols.plan_dft_fftw, "fftwf_plan_dft_1d")) return nullptr;
626  if (! Symbols::symbol (lib, symbols.plan_r2c_fftw, "fftwf_plan_dft_r2c_1d")) return nullptr;
627  if (! Symbols::symbol (lib, symbols.plan_c2r_fftw, "fftwf_plan_dft_c2r_1d")) return nullptr;
628  if (! Symbols::symbol (lib, symbols.destroy_fftw, "fftwf_destroy_plan")) return nullptr;
629 
630  if (! Symbols::symbol (lib, symbols.execute_dft_fftw, "fftwf_execute_dft")) return nullptr;
631  if (! Symbols::symbol (lib, symbols.execute_r2c_fftw, "fftwf_execute_dft_r2c")) return nullptr;
632  if (! Symbols::symbol (lib, symbols.execute_c2r_fftw, "fftwf_execute_dft_c2r")) return nullptr;
633  #endif
634 
635  return new FFTWImpl (static_cast<size_t> (order), std::move (lib), symbols);
636  }
637 
638  return nullptr;
639  }
640 
641  FFTWImpl (size_t orderToUse, DynamicLibrary&& libraryToUse, const Symbols& symbols)
642  : fftwLibrary (std::move (libraryToUse)), fftw (symbols), order (static_cast<size_t> (orderToUse))
643  {
644  ScopedLock lock (getFFTWPlanLock());
645 
646  auto n = (1u << order);
647  HeapBlock<Complex<float>> in (n), out (n);
648 
649  c2cForward = fftw.plan_dft_fftw (n, in.getData(), out.getData(), -1, unaligned | estimate);
650  c2cInverse = fftw.plan_dft_fftw (n, in.getData(), out.getData(), +1, unaligned | estimate);
651 
652  r2c = fftw.plan_r2c_fftw (n, (float*) in.getData(), in.getData(), unaligned | estimate);
653  c2r = fftw.plan_c2r_fftw (n, in.getData(), (float*) in.getData(), unaligned | estimate);
654  }
655 
656  ~FFTWImpl() override
657  {
658  ScopedLock lock (getFFTWPlanLock());
659 
660  fftw.destroy_fftw (c2cForward);
661  fftw.destroy_fftw (c2cInverse);
662  fftw.destroy_fftw (r2c);
663  fftw.destroy_fftw (c2r);
664  }
665 
666  void perform (const Complex<float>* input, Complex<float>* output, bool inverse) const noexcept override
667  {
668  if (inverse)
669  {
670  auto n = (1u << order);
671  fftw.execute_dft_fftw (c2cInverse, input, output);
672  FloatVectorOperations::multiply ((float*) output, 1.0f / static_cast<float> (n), (int) n << 1);
673  }
674  else
675  {
676  fftw.execute_dft_fftw (c2cForward, input, output);
677  }
678  }
679 
680  void performRealOnlyForwardTransform (float* inputOutputData, bool ignoreNegativeFreqs) const noexcept override
681  {
682  if (order == 0)
683  return;
684 
685  auto* out = reinterpret_cast<Complex<float>*> (inputOutputData);
686 
687  fftw.execute_r2c_fftw (r2c, inputOutputData, out);
688 
689  auto size = (1 << order);
690 
691  if (! ignoreNegativeFreqs)
692  for (auto i = size >> 1; i < size; ++i)
693  out[i] = std::conj (out[size - i]);
694  }
695 
696  void performRealOnlyInverseTransform (float* inputOutputData) const noexcept override
697  {
698  auto n = (1u << order);
699 
700  fftw.execute_c2r_fftw (c2r, (Complex<float>*) inputOutputData, inputOutputData);
701  FloatVectorOperations::multiply ((float*) inputOutputData, 1.0f / static_cast<float> (n), (int) n);
702  }
703 
704  //==============================================================================
705  // fftw's plan_* and destroy_* methods are NOT thread safe. So we need to share
706  // a lock between all instances of FFTWImpl
707  static CriticalSection& getFFTWPlanLock() noexcept
708  {
709  static CriticalSection cs;
710  return cs;
711  }
712 
713  //==============================================================================
714  DynamicLibrary fftwLibrary;
715  Symbols fftw;
716  size_t order;
717 
718  FFTWPlanRef c2cForward, c2cInverse, r2c, c2r;
719 };
720 
721 FFT::EngineImpl<FFTWImpl> fftwEngine;
722 #endif
723 
724 //==============================================================================
725 //==============================================================================
726 #if JUCE_DSP_USE_INTEL_MKL
727 struct IntelFFT : public FFT::Instance
728 {
729  static constexpr int priority = 8;
730 
731  static bool succeeded (MKL_LONG status) noexcept { return status == 0; }
732 
733  static IntelFFT* create (int orderToUse)
734  {
735  DFTI_DESCRIPTOR_HANDLE mklc2c, mklc2r;
736 
737  if (DftiCreateDescriptor (&mklc2c, DFTI_SINGLE, DFTI_COMPLEX, 1, 1 << orderToUse) == 0)
738  {
739  if (succeeded (DftiSetValue (mklc2c, DFTI_PLACEMENT, DFTI_NOT_INPLACE))
740  && succeeded (DftiSetValue (mklc2c, DFTI_BACKWARD_SCALE, 1.0f / static_cast<float> (1 << orderToUse)))
741  && succeeded (DftiCommitDescriptor (mklc2c)))
742  {
743  if (succeeded (DftiCreateDescriptor (&mklc2r, DFTI_SINGLE, DFTI_REAL, 1, 1 << orderToUse)))
744  {
745  if (succeeded (DftiSetValue (mklc2r, DFTI_PLACEMENT, DFTI_INPLACE))
746  && succeeded (DftiSetValue (mklc2r, DFTI_BACKWARD_SCALE, 1.0f / static_cast<float> (1 << orderToUse)))
747  && succeeded (DftiCommitDescriptor (mklc2r)))
748  {
749  return new IntelFFT (static_cast<size_t> (orderToUse), mklc2c, mklc2r);
750  }
751 
752  DftiFreeDescriptor (&mklc2r);
753  }
754  }
755 
756  DftiFreeDescriptor (&mklc2c);
757  }
758 
759  return {};
760  }
761 
762  IntelFFT (size_t orderToUse, DFTI_DESCRIPTOR_HANDLE c2cToUse, DFTI_DESCRIPTOR_HANDLE cr2ToUse)
763  : order (orderToUse), c2c (c2cToUse), c2r (cr2ToUse)
764  {}
765 
766  ~IntelFFT()
767  {
768  DftiFreeDescriptor (&c2c);
769  DftiFreeDescriptor (&c2r);
770  }
771 
772  void perform (const Complex<float>* input, Complex<float>* output, bool inverse) const noexcept override
773  {
774  if (inverse)
775  DftiComputeBackward (c2c, (void*) input, output);
776  else
777  DftiComputeForward (c2c, (void*) input, output);
778  }
779 
780  void performRealOnlyForwardTransform (float* inputOutputData, bool ignoreNegativeFreqs) const noexcept override
781  {
782  if (order == 0)
783  return;
784 
785  DftiComputeForward (c2r, inputOutputData);
786 
787  auto* out = reinterpret_cast<Complex<float>*> (inputOutputData);
788  auto size = (1 << order);
789 
790  if (! ignoreNegativeFreqs)
791  for (auto i = size >> 1; i < size; ++i)
792  out[i] = std::conj (out[size - i]);
793  }
794 
795  void performRealOnlyInverseTransform (float* inputOutputData) const noexcept override
796  {
797  DftiComputeBackward (c2r, inputOutputData);
798  }
799 
800  size_t order;
801  DFTI_DESCRIPTOR_HANDLE c2c, c2r;
802 };
803 
804 FFT::EngineImpl<IntelFFT> fftwEngine;
805 #endif
806 
807 //==============================================================================
808 //==============================================================================
809 FFT::FFT (int order)
810  : engine (FFT::Engine::createBestEngineForPlatform (order)),
811  size (1 << order)
812 {
813 }
814 
816 
817 void FFT::perform (const Complex<float>* input, Complex<float>* output, bool inverse) const noexcept
818 {
819  if (engine != nullptr)
820  engine->perform (input, output, inverse);
821 }
822 
823 void FFT::performRealOnlyForwardTransform (float* inputOutputData, bool ignoreNeagtiveFreqs) const noexcept
824 {
825  if (engine != nullptr)
826  engine->performRealOnlyForwardTransform (inputOutputData, ignoreNeagtiveFreqs);
827 }
828 
829 void FFT::performRealOnlyInverseTransform (float* inputOutputData) const noexcept
830 {
831  if (engine != nullptr)
832  engine->performRealOnlyInverseTransform (inputOutputData);
833 }
834 
835 void FFT::performFrequencyOnlyForwardTransform (float* inputOutputData) const noexcept
836 {
837  if (size == 1)
838  return;
839 
840  performRealOnlyForwardTransform (inputOutputData);
841  auto* out = reinterpret_cast<Complex<float>*> (inputOutputData);
842 
843  for (auto i = 0; i < size; ++i)
844  inputOutputData[i] = std::abs (out[i]);
845 
846  zeromem (&inputOutputData[size], static_cast<size_t> (size) * sizeof (float));
847 }
848 
849 } // namespace dsp
850 } // namespace juce
void performFrequencyOnlyForwardTransform(float *inputOutputData) const noexcept
Definition: juce_FFT.cpp:835
void perform(const Complex< float > *input, Complex< float > *output, bool inverse) const noexcept
Definition: juce_FFT.cpp:817
FFT(int order)
Definition: juce_FFT.cpp:809
Definition: juce_Uuid.h:140
static void JUCE_CALLTYPE multiply(float *dest, const float *src, int numValues) noexcept
static const FloatType pi
GenericScopedLock< SpinLock > ScopedLockType
Definition: juce_SpinLock.h:73
void performRealOnlyInverseTransform(float *inputOutputData) const noexcept
Definition: juce_FFT.cpp:829
void performRealOnlyForwardTransform(float *inputOutputData, bool dontCalculateNegativeFrequencies=false) const noexcept
Definition: juce_FFT.cpp:823