mirror/JUCE
1
0
Fork 0
mirror of https://github.com/juce-framework/JUCE.git synced 2026-01-10 23:44:24 +00:00
Code Activity

Added the JUCE DSP module

Browse source
This commit is contained in:
hogliux 2017-07-27 12:39:07 +01:00
parent 281c2fe2af
commit 244a944857
212 changed files with 37051 additions and 1301 deletions
Download patch file Download diff file Expand all files Collapse all files

1110
modules/juce_dsp/frequency/juce_Convolution.cpp Normal file
View file

File diff suppressed because it is too large Load diff

135
modules/juce_dsp/frequency/juce_Convolution.h Normal file
View file

@ -0,0 +1,135 @@
/*
==============================================================================
This file is part of the JUCE library.
Copyright (c) 2017 - ROLI Ltd.
JUCE is an open source library subject to commercial or open-source
licensing.
By using JUCE, you agree to the terms of both the JUCE 5 End-User License
Agreement and JUCE 5 Privacy Policy (both updated and effective as of the
27th April 2017).
End User License Agreement: www.juce.com/juce-5-licence
Privacy Policy: www.juce.com/juce-5-privacy-policy
Or: You may also use this code under the terms of the GPL v3 (see
www.gnu.org/licenses).
JUCE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY, AND ALL WARRANTIES, WHETHER
EXPRESSED OR IMPLIED, INCLUDING MERCHANTABILITY AND FITNESS FOR PURPOSE, ARE
DISCLAIMED.
==============================================================================
*/
/**
Performs stereo uniform-partitioned convolution of an input signal with an
impulse response in the frequency domain, using the juce FFT class.
It provides some thread-safe functions to load impulse responses as well,
from audio files or memory on the fly without any noticeable artefacts,
performing resampling and trimming if necessary.
The processing is equivalent to the time domain convolution done in the
class FIRFilter, with a FIRFilter::Coefficients object having as
coefficients the samples of the impulse response. However, it is more
efficient in general to do frequency domain convolution when the size of
the impulse response is higher than 64 samples.
see @FIRFilter, @FIRFilter::Coefficients, @FFT
*/
class JUCE_API Convolution
{
public:
//==============================================================================
/** Initialises an object for performing convolution in the frequency domain. */
Convolution();
/** Destructor. */
~Convolution();
//==============================================================================
/** Must be called before loading any impulse response, to provide to the
convolution the maximumBufferSize to handle, and the sample rate useful for
optional resampling.
*/
void prepare (const ProcessSpec&);
/** Resets the processing pipeline, ready to start a new stream of data. */
void reset() noexcept;
/** Performs the filter operation on the given set of samples, with optional
stereo processing.
*/
template <typename ProcessContext>
void process (const ProcessContext& context) noexcept
{
static_assert (std::is_same<typename ProcessContext::SampleType, float>::value,
"Convolution engine only supports single precision floating point data");
processSamples (context.getInputBlock(), context.getOutputBlock(), context.isBypassed);
}
//==============================================================================
/** This function loads an impulse response audio file from memory, added in a
JUCE project with the Projucer as binary data. It can load any of the audio
formats registered in JUCE, and performs some resampling and pre-processing
as well if needed.
Note : obviously, don't try to use this function on float samples, since the
data is supposed to be an audio file in its binary format, and be sure that
the original data is not going to move at all its memory location during the
process !!
@param sourceData the block of data to use as the stream's source
@param sourceDataSize the number of bytes in the source data block
@param wantsStereo requests to load both stereo channels or only one mono channel
@param size the expected size for the impulse response after loading
*/
void loadImpulseResponse (const void* sourceData, size_t sourceDataSize,
bool wantsStereo, size_t size);
/** This function loads an impulse response from an audio file on any drive. It
can load any of the audio formats registered in JUCE, and performs some
resampling and pre-processing as well if needed.
@param fileImpulseResponse the location of the audio file
@param wantsStereo requests to load both stereo channels or only one mono channel
@param size the expected size for the impulse response after loading
*/
void loadImpulseResponse (const File& fileImpulseResponse,
bool wantsStereo, size_t size);
/** This function loads an impulse response from an audio buffer, which is
copied before doing anything else. Performs some resampling and
pre-processing as well if needed.
@param buffer the AudioBuffer to use
@param bufferSampleRate the sampleRate of the data in the AudioBuffer
@param wantsStereo requests to load both stereo channels or only one mono channel
@param size the expected size for the impulse response after loading
*/
void copyAndLoadImpulseResponseFromBuffer (const AudioBuffer<float>& buffer, double bufferSampleRate,
bool wantsStereo, size_t size);
private:
//==============================================================================
struct Pimpl;
ScopedPointer<Pimpl> pimpl;
//==============================================================================
void processSamples (const AudioBlock<float>&, AudioBlock<float>&, bool isBypassed) noexcept;
//==============================================================================
double sampleRate;
bool currentIsBypassed = false;
LinearSmoothedValue<float> volumeDry[2], volumeWet[2];
AudioBlock<float> dryBuffer;
HeapBlock<char> dryBufferStorage;
//==============================================================================
JUCE_DECLARE_NON_COPYABLE_WITH_LEAK_DETECTOR (Convolution)
};

815
modules/juce_dsp/frequency/juce_FFT.cpp Normal file
View file

@ -0,0 +1,815 @@
/*
==============================================================================
This file is part of the JUCE library.
Copyright (c) 2017 - ROLI Ltd.
JUCE is an open source library subject to commercial or open-source
licensing.
By using JUCE, you agree to the terms of both the JUCE 5 End-User License
Agreement and JUCE 5 Privacy Policy (both updated and effective as of the
27th April 2017).
End User License Agreement: www.juce.com/juce-5-licence
Privacy Policy: www.juce.com/juce-5-privacy-policy
Or: You may also use this code under the terms of the GPL v3 (see
www.gnu.org/licenses).
JUCE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY, AND ALL WARRANTIES, WHETHER
EXPRESSED OR IMPLIED, INCLUDING MERCHANTABILITY AND FITNESS FOR PURPOSE, ARE
DISCLAIMED.
==============================================================================
*/
struct FFT::Instance
{
virtual ~Instance() {}
virtual void perform (const Complex<float>* input, Complex<float>* output, bool inverse) const noexcept = 0;
virtual void performRealOnlyForwardTransform (float*) const noexcept = 0;
virtual void performRealOnlyInverseTransform (float*) const noexcept = 0;
};
struct FFT::Engine
{
Engine (int priorityToUse) : enginePriority (priorityToUse)
{
EnginePriorityComparator comparator;
getEngines().addSorted (comparator, this);
}
virtual ~Engine() {}
virtual FFT::Instance* create (int order) const = 0;
//==============================================================================
static FFT::Instance* createBestEngineForPlatform (int order)
{
for (auto* engine : getEngines())
if (auto* instance = engine->create (order))
return instance;
jassertfalse; // This should never happen as the fallback engine should always work!
return nullptr;
}
private:
struct EnginePriorityComparator
{
static int compareElements (Engine* first, Engine* second) noexcept
{
// sort in reverse order
return DefaultElementComparator<int>::compareElements (second->enginePriority, first->enginePriority);
}
};
static Array<Engine*>& getEngines()
{
static Array<Engine*> engines;
return engines;
}
int enginePriority; // used so that faster engines have priority over slower ones
};
template <typename InstanceToUse>
struct FFT::EngineImpl : public FFT::Engine
{
EngineImpl() : FFT::Engine (InstanceToUse::priority) {}
FFT::Instance* create (int order) const override { return InstanceToUse::create (order); }
};
//==============================================================================
//==============================================================================
struct FFTFallback : public FFT::Instance
{
// this should have the least priority of all engines
static constexpr int priority = -1;
static FFTFallback* create (int order)
{
return new FFTFallback (order);
}
FFTFallback (int order)
{
configForward = new FFTConfig (1 << order, false);
configInverse = new FFTConfig (1 << order, true);
size = 1 << order;
}
void perform (const Complex<float>* input, Complex<float>* output, bool inverse) const noexcept override
{
if (size == 1)
{
*output = *input;
return;
}
const SpinLock::ScopedLockType sl(processLock);
jassert (configForward != nullptr);
if (inverse)
{
configInverse->perform (input, output);
const float scaleFactor = 1.0f / size;
for (int i = 0; i < size; ++i)
output[i] *= scaleFactor;
}
else
{
configForward->perform (input, output);
}
}
const size_t maxFFTScratchSpaceToAlloca = 256 * 1024;
void performRealOnlyForwardTransform (float* d) const noexcept override
{
if (size == 1)
return;
const size_t scratchSize = 16 + sizeof (Complex<float>) * (size_t) size;
if (scratchSize < maxFFTScratchSpaceToAlloca)
{
performRealOnlyForwardTransform (static_cast<Complex<float>*> (alloca (scratchSize)), d);
}
else
{
HeapBlock<char> heapSpace (scratchSize);
performRealOnlyForwardTransform (reinterpret_cast<Complex<float>*> (heapSpace.getData()), d);
}
}
void performRealOnlyInverseTransform (float* d) const noexcept override
{
if (size == 1)
return;
const size_t scratchSize = 16 + sizeof (Complex<float>) * (size_t) size;
if (scratchSize < maxFFTScratchSpaceToAlloca)
{
performRealOnlyInverseTransform (static_cast<Complex<float>*> (alloca (scratchSize)), d);
}
else
{
HeapBlock<char> heapSpace (scratchSize);
performRealOnlyInverseTransform (reinterpret_cast<Complex<float>*> (heapSpace.getData()), d);
}
}
void performRealOnlyForwardTransform (Complex<float>* scratch, float* d) const noexcept
{
for (int i = 0; i < size; ++i)
{
scratch[i].real (d[i]);
scratch[i].imag (0);
}
perform (scratch, reinterpret_cast<Complex<float>*> (d), false);
}
void performRealOnlyInverseTransform (Complex<float>* scratch, float* d) const noexcept
{
perform (reinterpret_cast<const Complex<float>*> (d), scratch, true);
for (int i = 0; i < size; ++i)
{
d[i] = scratch[i].real();
d[i + size] = scratch[i].imag();
}
}
//==============================================================================
struct FFTConfig
{
FFTConfig (int sizeOfFFT, bool isInverse)
: fftSize (sizeOfFFT), inverse (isInverse), twiddleTable ((size_t) sizeOfFFT)
{
const double inverseFactor = (inverse ? 2.0 : -2.0) * double_Pi / (double) fftSize;
if (fftSize <= 4)
{
for (int i = 0; i < fftSize; ++i)
{
const double phase = i * inverseFactor;
twiddleTable[i].real ((float) std::cos (phase));
twiddleTable[i].imag ((float) std::sin (phase));
}
}
else
{
for (int i = 0; i < fftSize / 4; ++i)
{
const double phase = i * inverseFactor;
twiddleTable[i].real ((float) std::cos (phase));
twiddleTable[i].imag ((float) std::sin (phase));
}
for (int i = fftSize / 4; i < fftSize / 2; ++i)
{
const int index = i - fftSize / 4;
twiddleTable[i].real (inverse ? -twiddleTable[index].imag() : twiddleTable[index].imag());
twiddleTable[i].imag (inverse ? twiddleTable[index].real() : -twiddleTable[index].real());
}
twiddleTable[fftSize / 2].real (-1.0f);
twiddleTable[fftSize / 2].imag (0.0f);
for (int i = fftSize / 2; i < fftSize; ++i)
{
const int index = fftSize / 2 - (i - fftSize / 2);
twiddleTable[i] = conj(twiddleTable[index]);
}
}
const int root = (int) std::sqrt ((double) fftSize);
int divisor = 4, n = fftSize;
for (int i = 0; i < numElementsInArray (factors); ++i)
{
while ((n % divisor) != 0)
{
if (divisor == 2) divisor = 3;
else if (divisor == 4) divisor = 2;
else divisor += 2;
if (divisor > root)
divisor = n;
}
n /= divisor;
jassert (divisor == 1 || divisor == 2 || divisor == 4);
factors[i].radix = divisor;
factors[i].length = n;
}
}
void perform (const Complex<float>* input, Complex<float>* output) const noexcept
{
perform (input, output, 1, 1, factors);
}
const int fftSize;
const bool inverse;
struct Factor { int radix, length; };
Factor factors[32];
HeapBlock<Complex<float>> twiddleTable;
void perform (const Complex<float>* input, Complex<float>* output, int stride, int strideIn, const Factor* facs) const noexcept
{
auto factor = *facs++;
auto* originalOutput = output;
auto* outputEnd = output + factor.radix * factor.length;
if (stride == 1 && factor.radix <= 5)
{
for (int i = 0; i < factor.radix; ++i)
perform (input + stride * strideIn * i, output + i * factor.length, stride * factor.radix, strideIn, facs);
butterfly (factor, output, stride);
return;
}
if (factor.length == 1)
{
do
{
*output++ = *input;
input += stride * strideIn;
}
while (output < outputEnd);
}
else
{
do
{
perform (input, output, stride * factor.radix, strideIn, facs);
input += stride * strideIn;
output += factor.length;
}
while (output < outputEnd);
}
butterfly (factor, originalOutput, stride);
}
void butterfly (const Factor factor, Complex<float>* data, int stride) const noexcept
{
switch (factor.radix)
{
case 1: break;
case 2: butterfly2 (data, stride, factor.length); return;
case 4: butterfly4 (data, stride, factor.length); return;
default: jassertfalse; break;
}
auto* scratch = static_cast<Complex<float>*> (alloca (sizeof (Complex<float>) * (size_t) factor.radix));
for (int i = 0; i < factor.length; ++i)
{
for (int k = i, q1 = 0; q1 < factor.radix; ++q1)
{
scratch[q1] = data[k];
k += factor.length;
}
for (int k = i, q1 = 0; q1 < factor.radix; ++q1)
{
int twiddleIndex = 0;
data[k] = scratch[0];
for (int q = 1; q < factor.radix; ++q)
{
twiddleIndex += stride * k;
if (twiddleIndex >= fftSize)
twiddleIndex -= fftSize;
data[k] += scratch[q] * twiddleTable[twiddleIndex];
}
k += factor.length;
}
}
}
void butterfly2 (Complex<float>* data, const int stride, const int length) const noexcept
{
auto* dataEnd = data + length;
auto* tw = twiddleTable.getData();
for (int i = length; --i >= 0;)
{
auto s = *dataEnd;
s *= (*tw);
tw += stride;
*dataEnd++ = *data - s;
*data++ += s;
}
}
void butterfly4 (Complex<float>* data, const int stride, const int length) const noexcept
{
auto lengthX2 = length * 2;
auto lengthX3 = length * 3;
auto strideX2 = stride * 2;
auto strideX3 = stride * 3;
auto* twiddle1 = twiddleTable.getData();
auto* twiddle2 = twiddle1;
auto* twiddle3 = twiddle1;
for (int i = length; --i >= 0;)
{
auto s0 = data[length] * *twiddle1;
auto s1 = data[lengthX2] * *twiddle2;
auto s2 = data[lengthX3] * *twiddle3;
auto s3 = s0; s3 += s2;
auto s4 = s0; s4 -= s2;
auto s5 = *data; s5 -= s1;
*data += s1;
data[lengthX2] = *data;
data[lengthX2] -= s3;
twiddle1 += stride;
twiddle2 += strideX2;
twiddle3 += strideX3;
*data += s3;
if (inverse)
{
data[length].real (s5.real() - s4.imag());
data[length].imag (s5.imag() + s4.real());
data[lengthX3].real (s5.real() + s4.imag());
data[lengthX3].imag (s5.imag() - s4.real());
}
else
{
data[length].real (s5.real() + s4.imag());
data[length].imag (s5.imag() - s4.real());
data[lengthX3].real (s5.real() - s4.imag());
data[lengthX3].imag (s5.imag() + s4.real());
}
++data;
}
}
JUCE_DECLARE_NON_COPYABLE_WITH_LEAK_DETECTOR (FFTConfig)
};
//==============================================================================
SpinLock processLock;
ScopedPointer<FFTConfig> configForward, configInverse;
int size;
};
FFT::EngineImpl<FFTFallback> fftFallback;
//==============================================================================
//==============================================================================
#if (JUCE_MAC || JUCE_IOS) && JUCE_USE_VDSP_FRAMEWORK
struct AppleFFT : public FFT::Instance
{
static constexpr int priority = 5;
static AppleFFT* create (int order)
{
return new AppleFFT (order);
}
AppleFFT (int orderToUse)
: order (static_cast<vDSP_Length> (orderToUse)),
fftSetup (vDSP_create_fftsetup (order, 2)),
forwardNormalisation (.5f),
inverseNormalisation (1.0f / static_cast<float> (1 << order))
{}
~AppleFFT() override
{
if (fftSetup != nullptr)
{
vDSP_destroy_fftsetup (fftSetup);
fftSetup = nullptr;
}
}
void perform (const Complex<float>* input, Complex<float>* output, bool inverse) const noexcept override
{
auto size = (1 << order);
DSPSplitComplex splitInput (toSplitComplex (const_cast<Complex<float>*> (input)));
DSPSplitComplex splitOutput (toSplitComplex (output));
vDSP_fft_zop (fftSetup, &splitInput, 2, &splitOutput, 2,
order, inverse ? kFFTDirection_Inverse : kFFTDirection_Forward);
float factor = (inverse ? inverseNormalisation : forwardNormalisation * 2.0f);
vDSP_vsmul ((float*) output, 1, &factor, (float*) output, 1, static_cast<size_t> (size << 1));
}
void performRealOnlyForwardTransform (float* inoutData) const noexcept override
{
auto size = (1 << order);
auto* inout = reinterpret_cast<Complex<float>*> (inoutData);
auto splitInOut (toSplitComplex (inout));
inoutData[size] = 0.0f;
vDSP_fft_zrip (fftSetup, &splitInOut, 2, order, kFFTDirection_Forward);
vDSP_vsmul (inoutData, 1, &forwardNormalisation, inoutData, 1, static_cast<size_t> (size << 1));
mirrorResult (inout);
}
void performRealOnlyInverseTransform (float* inoutData) const noexcept override
{
auto* inout = reinterpret_cast<Complex<float>*> (inoutData);
auto size = (1 << order);
auto splitInOut (toSplitComplex (inout));
// Imaginary part of nyquist and DC frequencies are always zero
// so Apple uses the imaginary part of the DC frequency to store
// the real part of the nyquist frequency
if (size != 1)
inout[0] = Complex<float> (inout[0].real(), inout[size >> 1].real());
vDSP_fft_zrip (fftSetup, &splitInOut, 2, order, kFFTDirection_Inverse);
vDSP_vsmul (inoutData, 1, &inverseNormalisation, inoutData, 1, static_cast<size_t> (size << 1));
vDSP_vclr (inoutData + size, 1, static_cast<size_t> (size));
}
private:
//==============================================================================
void mirrorResult (Complex<float>* out) const noexcept
{
auto size = (1 << order);
auto i = size >> 1;
// Imaginary part of nyquist and DC frequencies are always zero
// so Apple uses the imaginary part of the DC frequency to store
// the real part of the nyquist frequency
out[i++] = { out[0].imag(), 0.0 };
out[0] = { out[0].real(), 0.0 };
for (; i < size; ++i)
out[i] = std::conj (out[size - i]);
}
static DSPSplitComplex toSplitComplex (Complex<float>* data) noexcept
{
// this assumes that Complex interleaves real and imaginary parts
// and is tightly packed.
return { reinterpret_cast<float*> (data),
reinterpret_cast<float*> (data) + 1};
}
//==============================================================================
vDSP_Length order;
FFTSetup fftSetup;
float forwardNormalisation, inverseNormalisation;
};
FFT::EngineImpl<AppleFFT> appleFFT;
#endif
//==============================================================================
//==============================================================================
#if JUCE_DSP_USE_SHARED_FFTW || JUCE_DSP_USE_STATIC_FFTW
struct FFTWImpl : public FFT::Instance
{
#if JUCE_DSP_USE_STATIC_FFTW
// if the JUCE developer has gone through the hassle of statically
// linking in fftw, they probably want to use it
static constexpr int priority = 10;
#else
static constexpr int priority = 3;
#endif
struct FFTWPlan;
using FFTWPlanRef = FFTWPlan*;
enum
{
measure = 0,
unaligned = (1 << 1),
estimate = (1 << 6)
};
struct Symbols
{
FFTWPlanRef (*plan_dft_fftw) (unsigned, Complex<float>*, Complex<float>*, int, unsigned);
FFTWPlanRef (*plan_r2c_fftw) (unsigned, float*, Complex<float>*, unsigned);
FFTWPlanRef (*plan_c2r_fftw) (unsigned, Complex<float>*, float*, unsigned);
void (*destroy_fftw) (FFTWPlanRef);
void (*execute_dft_fftw) (FFTWPlanRef, const Complex<float>*, Complex<float>*);
void (*execute_r2c_fftw) (FFTWPlanRef, float*, Complex<float>*);
void (*execute_c2r_fftw) (FFTWPlanRef, Complex<float>*, float*);
#if JUCE_DSP_USE_STATIC_FFTW
template <typename FuncPtr, typename ActualSymbolType>
static bool symbol (FuncPtr& dst, ActualSymbolType sym)
{
dst = reinterpret_cast<FuncPtr> (sym);
return true;
}
#else
template <typename FuncPtr>
static bool symbol (DynamicLibrary& lib, FuncPtr& dst, const char* name)
{
dst = reinterpret_cast<FuncPtr> (lib.getFunction (name));
return (dst != nullptr);
}
#endif
};
static FFTWImpl* create (int order)
{
DynamicLibrary lib;
#if ! JUCE_DSP_USE_STATIC_FFTW
#if JUCE_MAC
const char* libsuffix = "dylib";
#elif JUCE_WINDOWS
const char* libsuffix = "dll";
#else
const char* libsuffix = "so";
#endif
if (lib.open (String ("libfftw3f.") + libsuffix))
#endif
{
Symbols symbols;
#if JUCE_DSP_USE_STATIC_FFTW
if (! Symbols::symbol (symbols.plan_dft_fftw, fftwf_plan_dft_1d)) return nullptr;
if (! Symbols::symbol (symbols.plan_r2c_fftw, fftwf_plan_dft_r2c_1d)) return nullptr;
if (! Symbols::symbol (symbols.plan_c2r_fftw, fftwf_plan_dft_c2r_1d)) return nullptr;
if (! Symbols::symbol (symbols.destroy_fftw, fftwf_destroy_plan)) return nullptr;
if (! Symbols::symbol (symbols.execute_dft_fftw, fftwf_execute_dft)) return nullptr;
if (! Symbols::symbol (symbols.execute_r2c_fftw, fftwf_execute_dft_r2c)) return nullptr;
if (! Symbols::symbol (symbols.execute_c2r_fftw, fftwf_execute_dft_c2r)) return nullptr;
#else
if (! Symbols::symbol (lib, symbols.plan_dft_fftw, "fftwf_plan_dft_1d")) return nullptr;
if (! Symbols::symbol (lib, symbols.plan_r2c_fftw, "fftwf_plan_dft_r2c_1d")) return nullptr;
if (! Symbols::symbol (lib, symbols.plan_c2r_fftw, "fftwf_plan_dft_c2r_1d")) return nullptr;
if (! Symbols::symbol (lib, symbols.destroy_fftw, "fftwf_destroy_plan")) return nullptr;
if (! Symbols::symbol (lib, symbols.execute_dft_fftw, "fftwf_execute_dft")) return nullptr;
if (! Symbols::symbol (lib, symbols.execute_r2c_fftw, "fftwf_execute_dft_r2c")) return nullptr;
if (! Symbols::symbol (lib, symbols.execute_c2r_fftw, "fftwf_execute_dft_c2r")) return nullptr;
#endif
return new FFTWImpl (static_cast<size_t> (order), static_cast<DynamicLibrary&&> (lib), symbols);
}
return nullptr;
}
FFTWImpl (size_t orderToUse, DynamicLibrary&& libraryToUse, const Symbols& symbols)
: fftwLibrary (std::move (libraryToUse)), fftw (symbols), order (static_cast<size_t> (orderToUse))
{
auto n = (1u << order);
HeapBlock<Complex<float>> in (n), out (n);
c2cForward = fftw.plan_dft_fftw (n, in.getData(), out.getData(), -1, unaligned | estimate);
c2cInverse = fftw.plan_dft_fftw (n, in.getData(), out.getData(), +1, unaligned | estimate);
r2c = fftw.plan_r2c_fftw (n, (float*) in.getData(), in.getData(), unaligned | estimate);
c2r = fftw.plan_c2r_fftw (n, in.getData(), (float*) in.getData(), unaligned | estimate);
}
~FFTWImpl() override
{
fftw.destroy_fftw (c2cForward);
fftw.destroy_fftw (c2cInverse);
fftw.destroy_fftw (r2c);
fftw.destroy_fftw (c2r);
}
void perform (const Complex<float>* input, Complex<float>* output, bool inverse) const noexcept override
{
if (inverse)
{
auto n = (1u << order);
fftw.execute_dft_fftw (c2cInverse, input, output);
FloatVectorOperations::multiply ((float*) output, 1.0f / static_cast<float> (n), (int) n << 1);
}
else
{
fftw.execute_dft_fftw (c2cForward, input, output);
}
}
void performRealOnlyForwardTransform (float* inputOutputData) const noexcept override
{
if (order == 0)
return;
auto* out = reinterpret_cast<Complex<float>*> (inputOutputData);
fftw.execute_r2c_fftw (r2c, inputOutputData, out);
auto size = (1 << order);
for (auto i = size >> 1; i < size; ++i)
out[i] = std::conj (out[size - i]);
}
void performRealOnlyInverseTransform (float* inputOutputData) const noexcept override
{
auto n = (1u << order);
fftw.execute_c2r_fftw (c2r, (Complex<float>*) inputOutputData, inputOutputData);
FloatVectorOperations::multiply ((float*) inputOutputData, 1.0f / static_cast<float> (n), (int) n);
}
//==============================================================================
DynamicLibrary fftwLibrary;
Symbols fftw;
size_t order;
FFTWPlanRef c2cForward, c2cInverse, r2c, c2r;
};
FFT::EngineImpl<FFTWImpl> fftwEngine;
#endif
//==============================================================================
//==============================================================================
#if JUCE_DSP_USE_INTEL_MKL
struct IntelFFT : public FFT::Instance
{
static constexpr int priority = 8;
static bool succeeded (MKL_LONG status) noexcept { return status == 0; }
static IntelFFT* create (int orderToUse)
{
DFTI_DESCRIPTOR_HANDLE mklc2c, mklc2r;
if (DftiCreateDescriptor (&mklc2c, DFTI_SINGLE, DFTI_COMPLEX, 1, 1 << orderToUse) == 0)
{
if (succeeded (DftiSetValue (mklc2c, DFTI_PLACEMENT, DFTI_NOT_INPLACE))
&& succeeded (DftiSetValue (mklc2c, DFTI_BACKWARD_SCALE, 1.0f / static_cast<float> (1 << orderToUse)))
&& succeeded (DftiCommitDescriptor (mklc2c)))
{
if (succeeded (DftiCreateDescriptor (&mklc2r, DFTI_SINGLE, DFTI_REAL, 1, 1 << orderToUse)))
{
if (succeeded (DftiSetValue (mklc2r, DFTI_PLACEMENT, DFTI_INPLACE))
&& succeeded (DftiSetValue (mklc2r, DFTI_BACKWARD_SCALE, 1.0f / static_cast<float> (1 << orderToUse)))
&& succeeded (DftiCommitDescriptor (mklc2r)))
{
return new IntelFFT (static_cast<size_t> (orderToUse), mklc2c, mklc2r);
}
DftiFreeDescriptor (&mklc2r);
}
}
DftiFreeDescriptor (&mklc2c);
}
return {};
}
IntelFFT (size_t orderToUse, DFTI_DESCRIPTOR_HANDLE c2cToUse, DFTI_DESCRIPTOR_HANDLE cr2ToUse)
: order (orderToUse), c2c (c2cToUse), c2r (cr2ToUse)
{}
~IntelFFT()
{
DftiFreeDescriptor (&c2c);
DftiFreeDescriptor (&c2r);
}
void perform (const Complex<float>* input, Complex<float>* output, bool inverse) const noexcept override
{
if (inverse)
DftiComputeBackward (c2c, (void*) input, output);
else
DftiComputeForward (c2c, (void*) input, output);
}
void performRealOnlyForwardTransform (float* inputOutputData) const noexcept override
{
if (order == 0)
return;
DftiComputeForward (c2r, inputOutputData);
auto* out = reinterpret_cast<Complex<float>*> (inputOutputData);
auto size = (1 << order);
for (auto i = size >> 1; i < size; ++i)
out[i] = std::conj (out[size - i]);
}
void performRealOnlyInverseTransform (float* inputOutputData) const noexcept override
{
DftiComputeBackward (c2r, inputOutputData);
}
size_t order;
DFTI_DESCRIPTOR_HANDLE c2c, c2r;
};
FFT::EngineImpl<IntelFFT> fftwEngine;
#endif
//==============================================================================
//==============================================================================
FFT::FFT (int order)
: engine (FFT::Engine::createBestEngineForPlatform (order)),
size (1 << order)
{}
FFT::~FFT() {}
void FFT::perform (const Complex<float>* input, Complex<float>* output, bool inverse) const noexcept
{
if (engine != nullptr)
engine->perform (input, output, inverse);
}
void FFT::performRealOnlyForwardTransform (float* inputOutputData) const noexcept
{
if (engine != nullptr)
engine->performRealOnlyForwardTransform (inputOutputData);
}
void FFT::performRealOnlyInverseTransform (float* inputOutputData) const noexcept
{
if (engine != nullptr)
engine->performRealOnlyInverseTransform (inputOutputData);
}
void FFT::performFrequencyOnlyForwardTransform (float* inputOutputData) const noexcept
{
if (size == 1)
return;
performRealOnlyForwardTransform (inputOutputData);
auto* out = reinterpret_cast<Complex<float>*> (inputOutputData);
for (auto i = 0; i < size; ++i)
inputOutputData[i] = std::abs (out[i]);
zeromem (&inputOutputData[size], sizeof (float) * static_cast<size_t> (size));
}

99
modules/juce_dsp/frequency/juce_FFT.h Normal file
View file

@ -0,0 +1,99 @@
/*
==============================================================================
This file is part of the JUCE library.
Copyright (c) 2017 - ROLI Ltd.
JUCE is an open source library subject to commercial or open-source
licensing.
By using JUCE, you agree to the terms of both the JUCE 5 End-User License
Agreement and JUCE 5 Privacy Policy (both updated and effective as of the
27th April 2017).
End User License Agreement: www.juce.com/juce-5-licence
Privacy Policy: www.juce.com/juce-5-privacy-policy
Or: You may also use this code under the terms of the GPL v3 (see
www.gnu.org/licenses).
JUCE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY, AND ALL WARRANTIES, WHETHER
EXPRESSED OR IMPLIED, INCLUDING MERCHANTABILITY AND FITNESS FOR PURPOSE, ARE
DISCLAIMED.
==============================================================================
*/
/**
Performs a fast fourier transform.
This is only a simple low-footprint implementation and isn't tuned for speed - it may
be useful for simple applications where one of the more complex FFT libraries would be
overkill. (But in the future it may end up becoming optimised of course...)
The FFT class itself contains lookup tables, so there's some overhead in creating
one, you should create and cache an FFT object for each size/direction of transform
that you need, and re-use them to perform the actual operation.
*/
class JUCE_API FFT
{
public:
//==============================================================================
/** Initialises an object for performing forward and inverse FFT with the given size.
The number of points the FFT will operate on will be 2 ^ order.
*/
FFT (int order);
/** Destructor. */
~FFT();
//==============================================================================
/** Performs an out-of-place FFT, either forward or inverse.
The arrays must contain at least getSize() elements.
*/
void perform (const Complex<float> *input, Complex<float> * output, bool inverse) const noexcept;
/** Performs an in-place forward transform on a block of real data.
The size of the array passed in must be 2 * getSize(), and the first half
should contain your raw input sample data. On return, the array will contain
size complex real + imaginary parts data interleaved, and can be passed to
performRealOnlyInverseTransform() in order to convert it back to reals.
*/
void performRealOnlyForwardTransform (float* inputOutputData) const noexcept;
/** Performs a reverse operation to data created in performRealOnlyForwardTransform().
The size of the array passed in must be 2 * getSize(), containing size complex
real and imaginary parts interleaved numbers. On return, the first half of the
array will contain the reconstituted samples.
*/
void performRealOnlyInverseTransform (float* inputOutputData) const noexcept;
/** Takes an array and simply transforms it to the magnitude frequency response
spectrum. This may be handy for things like frequency displays or analysis.
The size of the array passed in must be 2 * getSize().
*/
void performFrequencyOnlyForwardTransform (float* inputOutputData) const noexcept;
/** Returns the number of data points that this FFT was created to work with. */
int getSize() const noexcept { return size; }
//==============================================================================
#ifndef DOXYGEN
/* internal */
struct Instance;
template <typename> struct EngineImpl;
#endif
private:
//==============================================================================
struct Engine;
ScopedPointer<Instance> engine;
int size;
//==============================================================================
JUCE_DECLARE_NON_COPYABLE_WITH_LEAK_DETECTOR (FFT)
};

198
modules/juce_dsp/frequency/juce_FFT_test.cpp Normal file
View file

@ -0,0 +1,198 @@
/*
==============================================================================
This file is part of the JUCE library.
Copyright (c) 2017 - ROLI Ltd.
JUCE is an open source library subject to commercial or open-source
licensing.
By using JUCE, you agree to the terms of both the JUCE 5 End-User License
Agreement and JUCE 5 Privacy Policy (both updated and effective as of the
27th April 2017).
End User License Agreement: www.juce.com/juce-5-licence
Privacy Policy: www.juce.com/juce-5-privacy-policy
Or: You may also use this code under the terms of the GPL v3 (see
www.gnu.org/licenses).
JUCE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY, AND ALL WARRANTIES, WHETHER
EXPRESSED OR IMPLIED, INCLUDING MERCHANTABILITY AND FITNESS FOR PURPOSE, ARE
DISCLAIMED.
==============================================================================
*/
struct FFTUnitTest : public UnitTest
{
FFTUnitTest() : UnitTest("FFT") {}
static void fillRandom (Random& random, Complex<float>* buffer, size_t n)
{
for (size_t i = 0; i < n; ++i)
buffer[i] = Complex<float> ((2.0f * random.nextFloat()) - 1.0f,
(2.0f * random.nextFloat()) - 1.0f);
}
static void fillRandom (Random& random, float* buffer, size_t n)
{
for (size_t i = 0; i < n; ++i)
buffer[i] = (2.0f * random.nextFloat()) - 1.0f;
}
static Complex<float> freqConvolution (const Complex<float>* in, float freq, size_t n)
{
Complex<float> sum (0.0, 0.0);
for (size_t i = 0; i < n; ++i)
sum += in[i] * exp (Complex<float> (0, static_cast<float> (i) * freq));
return sum;
}
static void performReferenceFourier (const Complex<float>* in, Complex<float>* out,
size_t n, bool reverve)
{
float base_freq = static_cast<float>(((reverve ? 1.0 : -1.0) * 2.0 * double_Pi)
/ static_cast<float> (n));
for (size_t i = 0; i < n; ++i)
out[i] = freqConvolution (in, static_cast<float>(i) * base_freq, n);
}
static void performReferenceFourier (const float* in, Complex<float>* out,
size_t n, bool reverve)
{
HeapBlock<Complex<float>> buffer (n);
for (size_t i = 0; i < n; ++i)
buffer.getData()[i] = Complex<float> (in[i], 0.0f);
float base_freq = static_cast<float>(((reverve ? 1.0 : -1.0) * 2.0 * double_Pi)
/ static_cast<float> (n));
for (size_t i = 0; i < n; ++i)
out[i] = freqConvolution (buffer.getData(), static_cast<float>(i) * base_freq, n);
}
//==============================================================================
template <typename Type>
static bool checkArrayIsSimilar (Type* a, Type* b, size_t n) noexcept
{
for (size_t i = 0; i < n; ++i)
if (std::abs (a[i] - b[i]) > 1e-3f)
return false;
return true;
}
struct RealTest
{
static void run (FFTUnitTest& u)
{
Random random (378272);
for (size_t order = 0; order <= 8; ++order)
{
auto n = (1u << order);
FFT fft ((int) order);
HeapBlock<float> input (n);
HeapBlock<Complex<float>> reference (n), output (n);
fillRandom (random, input.getData(), n);
performReferenceFourier (input.getData(), reference.getData(), n, false);
// fill only first half with real numbers
zeromem (output.getData(), n * sizeof (Complex<float>));
memcpy (reinterpret_cast<float*> (output.getData()), input.getData(), n * sizeof (float));
fft.performRealOnlyForwardTransform ((float*) output.getData());
u.expect (checkArrayIsSimilar (reference.getData(), output.getData(), n));
memcpy (output.getData(), reference.getData(), n * sizeof (Complex<float>));
fft.performRealOnlyInverseTransform ((float*) output.getData());
u.expect (checkArrayIsSimilar ((float*) output.getData(), input.getData(), n));
}
}
};
struct FrequencyOnlyTest
{
static void run(FFTUnitTest& u)
{
Random random (378272);
for (size_t order = 0; order <= 8; ++order)
{
auto n = (1u << order);
FFT fft ((int) order);
HeapBlock<float> inout (n << 1), reference (n << 1);
HeapBlock<Complex<float> > frequency (n);
fillRandom (random, inout.getData(), n);
zeromem (reference.getData(), sizeof (float) * (n << 1));
performReferenceFourier (inout.getData(), frequency.getData(), n, false);
for (size_t i = 0; i < n; ++i)
reference.getData()[i] = std::abs (frequency.getData()[i]);
fft.performFrequencyOnlyForwardTransform (inout.getData());
u.expect (checkArrayIsSimilar (inout.getData(), reference.getData(), n));
}
}
};
struct ComplexTest
{
static void run(FFTUnitTest& u)
{
Random random (378272);
for (size_t order = 0; order <= 7; ++order)
{
auto n = (1u << order);
FFT fft ((int) order);
HeapBlock<Complex<float> > input (n), buffer (n), output (n), reference (n);
fillRandom (random, input.getData(), n);
performReferenceFourier (input.getData(), reference.getData(), n, false);
memcpy (buffer.getData(), input.getData(), sizeof (Complex<float>) * n);
fft.perform (buffer.getData(), output.getData(), false);
u.expect (checkArrayIsSimilar (output.getData(), reference.getData(), n));
memcpy (buffer.getData(), reference.getData(), sizeof (Complex<float>) * n);
fft.perform (buffer.getData(), output.getData(), true);
u.expect (checkArrayIsSimilar (output.getData(), input.getData(), n));
}
}
};
template <class TheTest>
void runTestForAllTypes (const char* unitTestName)
{
beginTest (unitTestName);
TheTest::run (*this);
}
void runTest() override
{
runTestForAllTypes<RealTest> ("Real input numbers Test");
runTestForAllTypes<FrequencyOnlyTest> ("Frequency only Test");
runTestForAllTypes<ComplexTest> ("Complex input numbers Test");
}
};
static FFTUnitTest fftUnitTest;

186
modules/juce_dsp/frequency/juce_Windowing.cpp Normal file
View file

@ -0,0 +1,186 @@
/*
==============================================================================
This file is part of the JUCE library.
Copyright (c) 2017 - ROLI Ltd.
JUCE is an open source library subject to commercial or open-source
licensing.
By using JUCE, you agree to the terms of both the JUCE 5 End-User License
Agreement and JUCE 5 Privacy Policy (both updated and effective as of the
27th April 2017).
End User License Agreement: www.juce.com/juce-5-licence
Privacy Policy: www.juce.com/juce-5-privacy-policy
Or: You may also use this code under the terms of the GPL v3 (see
www.gnu.org/licenses).
JUCE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY, AND ALL WARRANTIES, WHETHER
EXPRESSED OR IMPLIED, INCLUDING MERCHANTABILITY AND FITNESS FOR PURPOSE, ARE
DISCLAIMED.
==============================================================================
*/
template <typename FloatType>
static inline FloatType ncos (size_t order, size_t i, size_t size) noexcept
{
return std::cos (static_cast<FloatType> (order * i)
* MathConstants<FloatType>::pi / static_cast<FloatType> (size - 1));
}
template <typename FloatType>
WindowingFunction<FloatType>::WindowingFunction (size_t size, WindowingMethod type, bool normalize, FloatType beta)
{
fillWindowingTables (size, type, normalize, beta);
}
template <typename FloatType>
void WindowingFunction<FloatType>::fillWindowingTables (size_t size, WindowingMethod type,
bool normalize, FloatType beta) noexcept
{
windowTable.resize (static_cast<int> (size));
fillWindowingTables (windowTable.getRawDataPointer(), size, type, normalize, beta);
}
template <typename FloatType>
void WindowingFunction<FloatType>::fillWindowingTables (FloatType* samples, size_t size,
WindowingMethod type, bool normalize,
FloatType beta) noexcept
{
switch (type)
{
case rectangular:
{
for (size_t i = 0; i < size; ++i)
samples[i] = static_cast<FloatType> (1);
}
break;
case triangular:
{
auto halfSlots = static_cast<FloatType> (0.5) * static_cast<FloatType> (size - 1);
for (size_t i = 0; i < size; ++i)
samples[i] = static_cast<FloatType> (1.0) - std::fabs ((static_cast<FloatType> (i) - halfSlots) / halfSlots);
}
break;
case hann:
{
for (size_t i = 0; i < size; ++i)
{
auto cos2 = ncos<FloatType> (2, i, size);
samples[i] = static_cast<FloatType> (0.5 - 0.5 * cos2);
}
}
break;
case hamming:
{
for (size_t i = 0; i < size; ++i)
{
auto cos2 = ncos<FloatType> (2, i, size);
samples[i] = static_cast<FloatType> (0.54 - 0.46 * cos2);
}
}
break;
case blackmann:
{
constexpr FloatType alpha = 0.16f;
for (size_t i = 0; i < size; ++i)
{
auto cos2 = ncos<FloatType> (2, i, size);
auto cos4 = ncos<FloatType> (4, i, size);
samples[i] = static_cast<FloatType> (0.5 * (1 - alpha) - 0.5 * cos2 + 0.5 * alpha * cos4);
}
}
break;
case blackmannHarris:
{
for (size_t i = 0; i < size; ++i)
{
auto cos2 = ncos<FloatType> (2, i, size);
auto cos4 = ncos<FloatType> (4, i, size);
auto cos6 = ncos<FloatType> (6, i, size);
samples[i] = static_cast<FloatType> (0.35875 - 0.48829 * cos2 + 0.14128 * cos4 - 0.01168 * cos6);
}
}
break;
case flatTop:
{
for (size_t i = 0; i < size; ++i)
{
auto cos2 = ncos<FloatType> (2, i, size);
auto cos4 = ncos<FloatType> (4, i, size);
auto cos6 = ncos<FloatType> (6, i, size);
auto cos8 = ncos<FloatType> (8, i, size);
samples[i] = static_cast<FloatType> (1.0 - 1.93 * cos2 + 1.29 * cos4 - 0.388 * cos6 + 0.028 * cos8);
}
}
break;
case kaiser:
{
const double factor = 1.0 / SpecialFunctions::besselI0 (beta);
for (size_t i = 0; i < size; ++i)
samples[i] = static_cast<FloatType> (SpecialFunctions::besselI0 (beta * std::sqrt (1.0 - std::pow ((i - 0.5 * (size - 1.0))
/ ( 0.5 * (size - 1.0)), 2.0)))
* factor);
}
break;
default:
jassertfalse;
break;
}
// DC frequency amplitude must be one
if (normalize)
{
FloatType sum = {};
for (size_t i = 0; i < size; ++i)
sum += samples[i];
auto factor = static_cast<FloatType> (size) / sum;
FloatVectorOperations::multiply (samples, factor, static_cast<int> (size));
}
}
template <typename FloatType>
void WindowingFunction<FloatType>::multiplyWithWindowingTable (FloatType* samples, size_t size) noexcept
{
FloatVectorOperations::multiply (samples, windowTable.getRawDataPointer(), jmin (static_cast<int> (size), windowTable.size()));
}
template <typename FloatType>
const char* WindowingFunction<FloatType>::getWindowingMethodName (WindowingMethod type) noexcept
{
switch (type)
{
case rectangular: return "Rectangular";
case triangular: return "Triangular";
case hann: return "Hann";
case hamming: return "Hamming";
case blackmann: return "Blackmann";
case blackmannHarris: return "Blackmann-Harris";
case flatTop: return "FlatTop";
case kaiser: return "Kaiser";
default: jassertfalse; return "";
}
}
template struct WindowingFunction<float>;
template struct WindowingFunction<double>;

73
modules/juce_dsp/frequency/juce_Windowing.h Normal file
View file

@ -0,0 +1,73 @@
/*
==============================================================================
This file is part of the JUCE library.
Copyright (c) 2017 - ROLI Ltd.
JUCE is an open source library subject to commercial or open-source
licensing.
By using JUCE, you agree to the terms of both the JUCE 5 End-User License
Agreement and JUCE 5 Privacy Policy (both updated and effective as of the
27th April 2017).
End User License Agreement: www.juce.com/juce-5-licence
Privacy Policy: www.juce.com/juce-5-privacy-policy
Or: You may also use this code under the terms of the GPL v3 (see
www.gnu.org/licenses).
JUCE IS PROVIDED "AS IS" WITHOUT ANY WARRANTY, AND ALL WARRANTIES, WHETHER
EXPRESSED OR IMPLIED, INCLUDING MERCHANTABILITY AND FITNESS FOR PURPOSE, ARE
DISCLAIMED.
==============================================================================
*/
/**
A class which provides multiple windowing functions useful for filter design
and spectrum analyzers
*/
template <typename FloatType>
struct WindowingFunction
{
enum WindowingMethod
{
rectangular = 0,
triangular,
hann,
hamming,
blackmann,
blackmannHarris,
flatTop,
kaiser,
numWindowingMethods
};
//==============================================================================
WindowingFunction (size_t size, WindowingMethod,
bool normalize = true, FloatType beta = 0);
//==============================================================================
/** Fills the content of an array with a given windowing method table */
void fillWindowingTables (size_t size, WindowingMethod type,
bool normalize = true, FloatType beta = 0) noexcept;
/** Fills the content of an array with a given windowing method table */
static void fillWindowingTables (FloatType* samples, size_t size, WindowingMethod,
bool normalize = true, FloatType beta = 0) noexcept;
/** Multiply the content of a buffer with the given window */
void multiplyWithWindowingTable (FloatType* samples, size_t size) noexcept;
/** Returns the name of a given windowing method */
static const char* getWindowingMethodName (WindowingMethod) noexcept;
private:
//==============================================================================
Array<FloatType> windowTable;
JUCE_DECLARE_NON_COPYABLE_WITH_LEAK_DETECTOR (WindowingFunction)
};