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;
42 Engine (
int priorityToUse) : enginePriority (priorityToUse)
44 auto& list = getEngines();
46 std::sort (list.begin(), list.end(), [] (Engine* a, Engine* b) { return b->enginePriority < a->enginePriority; });
51 virtual FFT::Instance* create (
int order)
const = 0;
54 static FFT::Instance* createBestEngineForPlatform (
int order)
56 for (
auto* engine : getEngines())
57 if (auto* instance = engine->create (order))
65 static Array<Engine*>& getEngines()
67 static Array<Engine*> engines;
74template <
typename InstanceToUse>
75struct FFT::EngineImpl :
public FFT::Engine
77 EngineImpl() : FFT::Engine (InstanceToUse::priority) {}
78 FFT::Instance* create (
int order)
const override {
return InstanceToUse::create (order); }
83struct FFTFallback :
public FFT::Instance
86 static constexpr int priority = -1;
88 static FFTFallback* create (
int order)
90 return new FFTFallback (order);
93 FFTFallback (
int order)
95 configForward.reset (
new FFTConfig (1 << order,
false));
96 configInverse.reset (
new FFTConfig (1 << order,
true));
101 void perform (
const Complex<float>* input, Complex<float>* output,
bool inverse)
const noexcept override
111 jassert (configForward !=
nullptr);
115 configInverse->perform (input, output);
117 const float scaleFactor = 1.0f / size;
119 for (
int i = 0; i < size; ++i)
120 output[i] *= scaleFactor;
124 configForward->perform (input, output);
128 const size_t maxFFTScratchSpaceToAlloca = 256 * 1024;
130 void performRealOnlyForwardTransform (
float* d,
bool)
const noexcept override
135 const size_t scratchSize = 16 + (size_t) size *
sizeof (Complex<float>);
137 if (scratchSize < maxFFTScratchSpaceToAlloca)
139 performRealOnlyForwardTransform (
static_cast<Complex<float>*
> (alloca (scratchSize)), d);
143 HeapBlock<char> heapSpace (scratchSize);
144 performRealOnlyForwardTransform (
reinterpret_cast<Complex<float>*
> (heapSpace.getData()), d);
148 void performRealOnlyInverseTransform (
float* d)
const noexcept override
153 const size_t scratchSize = 16 + (size_t) size *
sizeof (Complex<float>);
155 if (scratchSize < maxFFTScratchSpaceToAlloca)
157 performRealOnlyInverseTransform (
static_cast<Complex<float>*
> (alloca (scratchSize)), d);
161 HeapBlock<char> heapSpace (scratchSize);
162 performRealOnlyInverseTransform (
reinterpret_cast<Complex<float>*
> (heapSpace.getData()), d);
166 void performRealOnlyForwardTransform (Complex<float>* scratch,
float* d)
const noexcept
168 for (
int i = 0; i < size; ++i)
169 scratch[i] = { d[i], 0 };
171 perform (scratch,
reinterpret_cast<Complex<float>*
> (d),
false);
174 void performRealOnlyInverseTransform (Complex<float>* scratch,
float* d)
const noexcept
176 auto* input =
reinterpret_cast<Complex<float>*
> (d);
178 for (
auto i = size >> 1; i < size; ++i)
179 input[i] = std::conj (input[size - i]);
181 perform (input, scratch,
true);
183 for (
int i = 0; i < size; ++i)
185 d[i] = scratch[i].real();
186 d[i + size] = scratch[i].imag();
193 FFTConfig (
int sizeOfFFT,
bool isInverse)
194 : fftSize (sizeOfFFT), inverse (isInverse), twiddleTable ((size_t) sizeOfFFT)
200 for (
int i = 0; i < fftSize; ++i)
202 auto phase = i * inverseFactor;
204 twiddleTable[i] = { (float) std::cos (phase),
205 (float) std::sin (phase) };
210 for (
int i = 0; i < fftSize / 4; ++i)
212 auto phase = i * inverseFactor;
214 twiddleTable[i] = { (float) std::cos (phase),
215 (float) std::sin (phase) };
218 for (
int i = fftSize / 4; i < fftSize / 2; ++i)
220 auto other = twiddleTable[i - fftSize / 4];
222 twiddleTable[i] = { inverse ? -other.imag() : other.imag(),
223 inverse ? other.real() : -other.real() };
226 twiddleTable[fftSize / 2].real (-1.0f);
227 twiddleTable[fftSize / 2].imag (0.0f);
229 for (
int i = fftSize / 2; i < fftSize; ++i)
231 auto index = fftSize / 2 - (i - fftSize / 2);
232 twiddleTable[i] = conj(twiddleTable[index]);
236 auto root = (int) std::sqrt ((
double) fftSize);
237 int divisor = 4, n = fftSize;
239 for (
int i = 0; i < numElementsInArray (factors); ++i)
241 while ((n % divisor) != 0)
243 if (divisor == 2) divisor = 3;
244 else if (divisor == 4) divisor = 2;
253 jassert (divisor == 1 || divisor == 2 || divisor == 4);
254 factors[i].radix = divisor;
255 factors[i].length = n;
259 void perform (
const Complex<float>* input, Complex<float>* output)
const noexcept
261 perform (input, output, 1, 1, factors);
267 struct Factor {
int radix, length; };
269 HeapBlock<Complex<float>> twiddleTable;
271 void perform (
const Complex<float>* input, Complex<float>* output,
int stride,
int strideIn,
const Factor* facs)
const noexcept
273 auto factor = *facs++;
274 auto* originalOutput = output;
275 auto* outputEnd = output + factor.radix * factor.length;
277 if (stride == 1 && factor.radix <= 5)
279 for (
int i = 0; i < factor.radix; ++i)
280 perform (input + stride * strideIn * i, output + i * factor.length, stride * factor.radix, strideIn, facs);
282 butterfly (factor, output, stride);
286 if (factor.length == 1)
291 input += stride * strideIn;
293 while (output < outputEnd);
299 perform (input, output, stride * factor.radix, strideIn, facs);
300 input += stride * strideIn;
301 output += factor.length;
303 while (output < outputEnd);
306 butterfly (factor, originalOutput, stride);
309 void butterfly (
const Factor factor, Complex<float>* data,
int stride)
const noexcept
311 switch (factor.radix)
314 case 2: butterfly2 (data, stride, factor.length);
return;
315 case 4: butterfly4 (data, stride, factor.length);
return;
316 default: jassertfalse;
break;
319 auto* scratch =
static_cast<Complex<float>*
> (alloca ((
size_t) factor.radix * sizeof (Complex<float>)));
321 for (
int i = 0; i < factor.length; ++i)
323 for (
int k = i, q1 = 0; q1 < factor.radix; ++q1)
325 scratch[q1] = data[k];
329 for (
int k = i, q1 = 0; q1 < factor.radix; ++q1)
331 int twiddleIndex = 0;
332 data[k] = scratch[0];
334 for (
int q = 1; q < factor.radix; ++q)
336 twiddleIndex += stride * k;
338 if (twiddleIndex >= fftSize)
339 twiddleIndex -= fftSize;
341 data[k] += scratch[q] * twiddleTable[twiddleIndex];
349 void butterfly2 (Complex<float>* data,
const int stride,
const int length)
const noexcept
351 auto* dataEnd = data + length;
352 auto* tw = twiddleTable.getData();
354 for (
int i = length; --i >= 0;)
359 *dataEnd++ = *data - s;
364 void butterfly4 (Complex<float>* data,
const int stride,
const int length)
const noexcept
366 auto lengthX2 = length * 2;
367 auto lengthX3 = length * 3;
369 auto strideX2 = stride * 2;
370 auto strideX3 = stride * 3;
372 auto* twiddle1 = twiddleTable.getData();
373 auto* twiddle2 = twiddle1;
374 auto* twiddle3 = twiddle1;
376 for (
int i = length; --i >= 0;)
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;
386 data[lengthX2] = *data;
387 data[lengthX2] -= s3;
389 twiddle2 += strideX2;
390 twiddle3 += strideX3;
395 data[length] = { s5.real() - s4.imag(),
396 s5.imag() + s4.real() };
398 data[lengthX3] = { s5.real() + s4.imag(),
399 s5.imag() - s4.real() };
403 data[length] = { s5.real() + s4.imag(),
404 s5.imag() - s4.real() };
406 data[lengthX3] = { s5.real() - s4.imag(),
407 s5.imag() + s4.real() };
414 JUCE_DECLARE_NON_COPYABLE_WITH_LEAK_DETECTOR (FFTConfig)
418 SpinLock processLock;
419 std::unique_ptr<FFTConfig> configForward, configInverse;
423FFT::EngineImpl<FFTFallback> fftFallback;
427#if (JUCE_MAC || JUCE_IOS) && JUCE_USE_VDSP_FRAMEWORK
428struct AppleFFT :
public FFT::Instance
430 static constexpr int priority = 5;
432 static AppleFFT* create (
int order)
434 return new AppleFFT (order);
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))
446 if (fftSetup !=
nullptr)
448 vDSP_destroy_fftsetup (fftSetup);
453 void perform (
const Complex<float>* input, Complex<float>* output,
bool inverse)
const noexcept override
455 auto size = (1 << order);
457 DSPSplitComplex splitInput (toSplitComplex (
const_cast<Complex<float>*
> (input)));
458 DSPSplitComplex splitOutput (toSplitComplex (output));
460 vDSP_fft_zop (fftSetup, &splitInput, 2, &splitOutput, 2,
461 order, inverse ? kFFTDirection_Inverse : kFFTDirection_Forward);
463 float factor = (inverse ? inverseNormalisation : forwardNormalisation * 2.0f);
464 vDSP_vsmul ((
float*) output, 1, &factor, (
float*) output, 1,
static_cast<size_t> (size << 1));
467 void performRealOnlyForwardTransform (
float* inoutData,
bool ignoreNegativeFreqs)
const noexcept override
469 auto size = (1 << order);
470 auto* inout =
reinterpret_cast<Complex<float>*
> (inoutData);
471 auto splitInOut (toSplitComplex (inout));
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));
477 mirrorResult (inout, ignoreNegativeFreqs);
480 void performRealOnlyInverseTransform (
float* inoutData)
const noexcept override
482 auto* inout =
reinterpret_cast<Complex<float>*
> (inoutData);
483 auto size = (1 << order);
484 auto splitInOut (toSplitComplex (inout));
490 inout[0] = Complex<float> (inout[0].real(), inout[size >> 1].real());
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));
499 void mirrorResult (Complex<float>* out,
bool ignoreNegativeFreqs)
const noexcept
501 auto size = (1 << order);
507 out[i++] = { out[0].imag(), 0.0 };
508 out[0] = { out[0].real(), 0.0 };
510 if (! ignoreNegativeFreqs)
511 for (; i < size; ++i)
512 out[i] = std::conj (out[size - i]);
515 static DSPSplitComplex toSplitComplex (Complex<float>* data)
noexcept
519 return {
reinterpret_cast<float*
> (data),
520 reinterpret_cast<float*
> (data) + 1};
526 float forwardNormalisation, inverseNormalisation;
529FFT::EngineImpl<AppleFFT> appleFFT;
534#if JUCE_DSP_USE_SHARED_FFTW || JUCE_DSP_USE_STATIC_FFTW
536#if JUCE_DSP_USE_STATIC_FFTW
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*);
549struct FFTWImpl :
public FFT::Instance
551 #if JUCE_DSP_USE_STATIC_FFTW
554 static constexpr int priority = 10;
556 static constexpr int priority = 3;
560 using FFTWPlanRef = FFTWPlan*;
565 unaligned = (1 << 1),
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);
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*);
580 #if JUCE_DSP_USE_STATIC_FFTW
581 template <
typename FuncPtr,
typename ActualSymbolType>
582 static bool symbol (FuncPtr& dst, ActualSymbolType sym)
584 dst =
reinterpret_cast<FuncPtr
> (sym);
588 template <
typename FuncPtr>
589 static bool symbol (DynamicLibrary& lib, FuncPtr& dst,
const char* name)
591 dst =
reinterpret_cast<FuncPtr
> (lib.getFunction (name));
592 return (dst !=
nullptr);
597 static FFTWImpl* create (
int order)
601 #if ! JUCE_DSP_USE_STATIC_FFTW
603 auto libName =
"libfftw3f.dylib";
605 auto libName =
"libfftw3f.dll";
607 auto libName =
"libfftw3f.so";
610 if (lib.open (libName))
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;
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;
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;
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;
635 return new FFTWImpl (
static_cast<size_t> (order), std::move (lib), symbols);
641 FFTWImpl (
size_t orderToUse, DynamicLibrary&& libraryToUse,
const Symbols& symbols)
642 : fftwLibrary (std::move (libraryToUse)), fftw (symbols), order (static_cast<size_t> (orderToUse))
644 ScopedLock lock (getFFTWPlanLock());
646 auto n = (1u << order);
647 HeapBlock<Complex<float>> in (n), out (n);
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);
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);
658 ScopedLock lock (getFFTWPlanLock());
660 fftw.destroy_fftw (c2cForward);
661 fftw.destroy_fftw (c2cInverse);
662 fftw.destroy_fftw (r2c);
663 fftw.destroy_fftw (c2r);
666 void perform (
const Complex<float>* input, Complex<float>* output,
bool inverse)
const noexcept override
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);
676 fftw.execute_dft_fftw (c2cForward, input, output);
680 void performRealOnlyForwardTransform (
float* inputOutputData,
bool ignoreNegativeFreqs)
const noexcept override
685 auto* out =
reinterpret_cast<Complex<float>*
> (inputOutputData);
687 fftw.execute_r2c_fftw (r2c, inputOutputData, out);
689 auto size = (1 << order);
691 if (! ignoreNegativeFreqs)
692 for (
auto i = size >> 1; i < size; ++i)
693 out[i] = std::conj (out[size - i]);
696 void performRealOnlyInverseTransform (
float* inputOutputData)
const noexcept override
698 auto n = (1u << order);
700 fftw.execute_c2r_fftw (c2r, (Complex<float>*) inputOutputData, inputOutputData);
701 FloatVectorOperations::multiply ((
float*) inputOutputData, 1.0f /
static_cast<float> (n), (
int) n);
707 static CriticalSection& getFFTWPlanLock() noexcept
709 static CriticalSection cs;
714 DynamicLibrary fftwLibrary;
718 FFTWPlanRef c2cForward, c2cInverse, r2c, c2r;
721FFT::EngineImpl<FFTWImpl> fftwEngine;
726#if JUCE_DSP_USE_INTEL_MKL
727struct IntelFFT :
public FFT::Instance
729 static constexpr int priority = 8;
731 static bool succeeded (MKL_LONG status)
noexcept {
return status == 0; }
733 static IntelFFT* create (
int orderToUse)
735 DFTI_DESCRIPTOR_HANDLE mklc2c, mklc2r;
737 if (DftiCreateDescriptor (&mklc2c, DFTI_SINGLE, DFTI_COMPLEX, 1, 1 << orderToUse) == 0)
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)))
743 if (succeeded (DftiCreateDescriptor (&mklc2r, DFTI_SINGLE, DFTI_REAL, 1, 1 << orderToUse)))
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)))
749 return new IntelFFT (
static_cast<size_t> (orderToUse), mklc2c, mklc2r);
752 DftiFreeDescriptor (&mklc2r);
756 DftiFreeDescriptor (&mklc2c);
762 IntelFFT (
size_t orderToUse, DFTI_DESCRIPTOR_HANDLE c2cToUse, DFTI_DESCRIPTOR_HANDLE cr2ToUse)
763 : order (orderToUse), c2c (c2cToUse), c2r (cr2ToUse)
768 DftiFreeDescriptor (&c2c);
769 DftiFreeDescriptor (&c2r);
772 void perform (
const Complex<float>* input, Complex<float>* output,
bool inverse)
const noexcept override
775 DftiComputeBackward (c2c, (
void*) input, output);
777 DftiComputeForward (c2c, (
void*) input, output);
780 void performRealOnlyForwardTransform (
float* inputOutputData,
bool ignoreNegativeFreqs)
const noexcept override
785 DftiComputeForward (c2r, inputOutputData);
787 auto* out =
reinterpret_cast<Complex<float>*
> (inputOutputData);
788 auto size = (1 << order);
790 if (! ignoreNegativeFreqs)
791 for (
auto i = size >> 1; i < size; ++i)
792 out[i] = std::conj (out[size - i]);
795 void performRealOnlyInverseTransform (
float* inputOutputData)
const noexcept override
797 DftiComputeBackward (c2r, inputOutputData);
801 DFTI_DESCRIPTOR_HANDLE c2c, c2r;
804FFT::EngineImpl<IntelFFT> fftwEngine;
810 : engine (
FFT::Engine::createBestEngineForPlatform (order)),
819 if (engine !=
nullptr)
820 engine->perform (input, output, inverse);
825 if (engine !=
nullptr)
831 if (engine !=
nullptr)
843 for (
auto i = 0; i < size; ++i)
846 zeromem (&
inputOutputData[size],
static_cast<size_t> (size) *
sizeof (
float));
GenericScopedLock< SpinLock > ScopedLockType
void performRealOnlyInverseTransform(float *inputOutputData) const noexcept
void performRealOnlyForwardTransform(float *inputOutputData, bool dontCalculateNegativeFrequencies=false) const noexcept
void perform(const Complex< float > *input, Complex< float > *output, bool inverse) const noexcept
void performFrequencyOnlyForwardTransform(float *inputOutputData) const noexcept
static const FloatType pi