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