blumia/juicysfplugin
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

fix macOS build (following Projucer changes made in Windows, which removed /Applications/JUCE/modules from its headers). move JUCE headers under source control, so that Windows and macOS can both build against same version of JUCE. remove AUv3 target (I think it's an iOS thing, so it will never work with this macOS fluidsynth dylib).

Browse Source
This commit is contained in:
Alex Birch
2018-06-17 13:34:53 +01:00
parent a2be47c887
commit dff4d13a1d
1563 changed files with 601601 additions and 3466 deletions
Download Patch File Download Diff File Expand all files Collapse all files

265
modules/juce_dsp/maths/juce_FastMathApproximations.h Normal file
View File

@ -0,0 +1,265 @@
/*
==============================================================================
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.
==============================================================================
*/
namespace juce
{
namespace dsp
{
/**
This class contains various fast mathematical function approximations.
@tags{DSP}
*/
struct FastMathApproximations
{
/** Provides a fast approximation of the function cosh(x) using a Pade approximant
continued fraction, calculated sample by sample.
Note : this is an approximation which works on a limited range. You are
advised to use input values only between -5 and +5 for limiting the error.
*/
template <typename FloatType>
static FloatType cosh (FloatType x) noexcept
{
auto x2 = x * x;
auto numerator = -(39251520 + x2 * (18471600 + x2 * (1075032 + 14615 * x2)));
auto denominator = -39251520 + x2 * (1154160 + x2 * (-16632 + 127 * x2));
return numerator / denominator;
}
/** Provides a fast approximation of the function cosh(x) using a Pade approximant
continued fraction, calculated on a whole buffer.
Note : this is an approximation which works on a limited range. You are
advised to use input values only between -5 and +5 for limiting the error.
*/
template <typename FloatType>
static void cosh (FloatType* values, size_t numValues) noexcept
{
for (size_t i = 0; i < numValues; ++i)
values[i] = FastMathApproximations::cosh (values[i]);
}
/** Provides a fast approximation of the function sinh(x) using a Pade approximant
continued fraction, calculated sample by sample.
Note : this is an approximation which works on a limited range. You are
advised to use input values only between -5 and +5 for limiting the error.
*/
template <typename FloatType>
static FloatType sinh (FloatType x) noexcept
{
auto x2 = x * x;
auto numerator = -x * (11511339840 + x2 * (1640635920 + x2 * (52785432 + x2 * 479249)));
auto denominator = -11511339840 + x2 * (277920720 + x2 * (-3177720 + x2 * 18361));
return numerator / denominator;
}
/** Provides a fast approximation of the function sinh(x) using a Pade approximant
continued fraction, calculated on a whole buffer.
Note : this is an approximation which works on a limited range. You are
advised to use input values only between -5 and +5 for limiting the error.
*/
template <typename FloatType>
static void sinh (FloatType* values, size_t numValues) noexcept
{
for (size_t i = 0; i < numValues; ++i)
values[i] = FastMathApproximations::sinh (values[i]);
}
/** Provides a fast approximation of the function tanh(x) using a Pade approximant
continued fraction, calculated sample by sample.
Note : this is an approximation which works on a limited range. You are
advised to use input values only between -5 and +5 for limiting the error.
*/
template <typename FloatType>
static FloatType tanh (FloatType x) noexcept
{
auto x2 = x * x;
auto numerator = x * (135135 + x2 * (17325 + x2 * (378 + x2)));
auto denominator = 135135 + x2 * (62370 + x2 * (3150 + 28 * x2));
return numerator / denominator;
}
/** Provides a fast approximation of the function tanh(x) using a Pade approximant
continued fraction, calculated on a whole buffer.
Note : this is an approximation which works on a limited range. You are
advised to use input values only between -5 and +5 for limiting the error.
*/
template <typename FloatType>
static void tanh (FloatType* values, size_t numValues) noexcept
{
for (size_t i = 0; i < numValues; ++i)
values[i] = FastMathApproximations::tanh (values[i]);
}
//==============================================================================
/** Provides a fast approximation of the function cos(x) using a Pade approximant
continued fraction, calculated sample by sample.
Note : this is an approximation which works on a limited range. You are
advised to use input values only between -pi and +pi for limiting the error.
*/
template <typename FloatType>
static FloatType cos (FloatType x) noexcept
{
auto x2 = x * x;
auto numerator = -(-39251520 + x2 * (18471600 + x2 * (-1075032 + 14615 * x2)));
auto denominator = 39251520 + x2 * (1154160 + x2 * (16632 + x2 * 127));
return numerator / denominator;
}
/** Provides a fast approximation of the function cos(x) using a Pade approximant
continued fraction, calculated on a whole buffer.
Note : this is an approximation which works on a limited range. You are
advised to use input values only between -pi and +pi for limiting the error.
*/
template <typename FloatType>
static void cos (FloatType* values, size_t numValues) noexcept
{
for (size_t i = 0; i < numValues; ++i)
values[i] = FastMathApproximations::cos (values[i]);
}
/** Provides a fast approximation of the function sin(x) using a Pade approximant
continued fraction, calculated sample by sample.
Note : this is an approximation which works on a limited range. You are
advised to use input values only between -pi and +pi for limiting the error.
*/
template <typename FloatType>
static FloatType sin (FloatType x) noexcept
{
auto x2 = x * x;
auto numerator = -x * (-11511339840 + x2 * (1640635920 + x2 * (-52785432 + x2 * 479249)));
auto denominator = 11511339840 + x2 * (277920720 + x2 * (3177720 + x2 * 18361));
return numerator / denominator;
}
/** Provides a fast approximation of the function sin(x) using a Pade approximant
continued fraction, calculated on a whole buffer.
Note : this is an approximation which works on a limited range. You are
advised to use input values only between -pi and +pi for limiting the error.
*/
template <typename FloatType>
static void sin (FloatType* values, size_t numValues) noexcept
{
for (size_t i = 0; i < numValues; ++i)
values[i] = FastMathApproximations::sin (values[i]);
}
/** Provides a fast approximation of the function tan(x) using a Pade approximant
continued fraction, calculated sample by sample.
Note : this is an approximation which works on a limited range. You are
advised to use input values only between -pi/2 and +pi/2 for limiting the error.
*/
template <typename FloatType>
static FloatType tan (FloatType x) noexcept
{
auto x2 = x * x;
auto numerator = x * (-135135 + x2 * (17325 + x2 * (-378 + x2)));
auto denominator = -135135 + x2 * (62370 + x2 * (-3150 + 28 * x2));
return numerator / denominator;
}
/** Provides a fast approximation of the function tan(x) using a Pade approximant
continued fraction, calculated on a whole buffer.
Note : this is an approximation which works on a limited range. You are
advised to use input values only between -pi/2 and +pi/2 for limiting the error.
*/
template <typename FloatType>
static void tan (FloatType* values, size_t numValues) noexcept
{
for (size_t i = 0; i < numValues; ++i)
values[i] = FastMathApproximations::tan (values[i]);
}
//==============================================================================
/** Provides a fast approximation of the function exp(x) using a Pade approximant
continued fraction, calculated sample by sample.
Note : this is an approximation which works on a limited range. You are
advised to use input values only between -6 and +4 for limiting the error.
*/
template <typename FloatType>
static FloatType exp (FloatType x) noexcept
{
auto numerator = 1680 + x * (840 + x * (180 + x * (20 + x)));
auto denominator = 1680 + x *(-840 + x * (180 + x * (-20 + x)));
return numerator / denominator;
}
/** Provides a fast approximation of the function exp(x) using a Pade approximant
continued fraction, calculated on a whole buffer.
Note : this is an approximation which works on a limited range. You are
advised to use input values only between -6 and +4 for limiting the error.
*/
template <typename FloatType>
static void exp (FloatType* values, size_t numValues) noexcept
{
for (size_t i = 0; i < numValues; ++i)
values[i] = FastMathApproximations::exp (values[i]);
}
/** Provides a fast approximation of the function log(x+1) using a Pade approximant
continued fraction, calculated sample by sample.
Note : this is an approximation which works on a limited range. You are
advised to use input values only between -0.8 and +5 for limiting the error.
*/
template <typename FloatType>
static FloatType logNPlusOne (FloatType x) noexcept
{
auto numerator = x * (7560 + x * (15120 + x * (9870 + x * (2310 + x * 137))));
auto denominator = 7560 + x * (18900 + x * (16800 + x * (6300 + x * (900 + 30 * x))));
return numerator / denominator;
}
/** Provides a fast approximation of the function log(x+1) using a Pade approximant
continued fraction, calculated on a whole buffer.
Note : this is an approximation which works on a limited range. You are
advised to use input values only between -0.8 and +5 for limiting the error.
*/
template <typename FloatType>
static void logNPlusOne (FloatType* values, size_t numValues) noexcept
{
for (size_t i = 0; i < numValues; ++i)
values[i] = FastMathApproximations::logNPlusOne (values[i]);
}
};
} // namespace dsp
} // namespace juce

157
modules/juce_dsp/maths/juce_LookupTable.cpp Normal file
View File

@ -0,0 +1,157 @@
/*
==============================================================================
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.
==============================================================================
*/
namespace juce
{
namespace dsp
{
template <typename FloatType>
LookupTable<FloatType>::LookupTable()
{
data.resize (1);
}
template <typename FloatType>
LookupTable<FloatType>::LookupTable (const std::function<FloatType (size_t)>& functionToApproximate,
size_t numPointsToUse)
{
initialise (functionToApproximate, numPointsToUse);
}
//==============================================================================
template <typename FloatType>
void LookupTable<FloatType>::initialise (const std::function<FloatType (size_t)>& functionToApproximate,
size_t numPointsToUse)
{
data.resize (static_cast<int> (getRequiredBufferSize (numPointsToUse)));
for (size_t i = 0; i < numPointsToUse; ++i)
{
auto value = functionToApproximate (i);
jassert (! std::isnan (value));
jassert (! std::isinf (value));
// Make sure functionToApproximate returns a sensible value for the entire specified range.
// E.g., this won't work for zero: [] (size_t i) { return 1.0f / i; }
data.getReference (static_cast<int> (i)) = value;
}
prepare();
}
template <typename FloatType>
void LookupTable<FloatType>::prepare() noexcept
{
auto guardIndex = static_cast<int> (getGuardIndex());
data.getReference (guardIndex) = data.getUnchecked (guardIndex - 1);
}
template <typename FloatType>
void LookupTableTransform<FloatType>::initialise (const std::function<FloatType (FloatType)>& functionToApproximate,
FloatType minInputValueToUse,
FloatType maxInputValueToUse,
size_t numPoints)
{
jassert (maxInputValueToUse > minInputValueToUse);
minInputValue = minInputValueToUse;
maxInputValue = maxInputValueToUse;
scaler = FloatType (numPoints - 1) / (maxInputValueToUse - minInputValueToUse);
offset = -minInputValueToUse * scaler;
const auto initFn = [functionToApproximate, minInputValueToUse, maxInputValueToUse, numPoints] (size_t i)
{
return functionToApproximate (
jlimit (
minInputValueToUse, maxInputValueToUse,
jmap (FloatType (i), FloatType (0), FloatType (numPoints - 1), minInputValueToUse, maxInputValueToUse))
);
};
lookupTable.initialise (initFn, numPoints);
}
//==============================================================================
template <typename FloatType>
double LookupTableTransform<FloatType>::calculateMaxRelativeError (const std::function<FloatType (FloatType)>& functionToApproximate,
FloatType minInputValue,
FloatType maxInputValue,
size_t numPoints,
size_t numTestPoints)
{
jassert (maxInputValue > minInputValue);
if (numTestPoints == 0)
numTestPoints = 100 * numPoints; // use default
LookupTableTransform transform (functionToApproximate, minInputValue, maxInputValue, numPoints);
double maxError = 0;
for (size_t i = 0; i < numTestPoints; ++i)
{
auto inputValue = jmap (FloatType (i), FloatType (0), FloatType (numTestPoints - 1), minInputValue, maxInputValue);
auto approximatedOutputValue = transform.processSample (inputValue);
auto referenceOutputValue = functionToApproximate (inputValue);
maxError = jmax (maxError, calculateRelativeDifference ((double) referenceOutputValue, (double) approximatedOutputValue));
}
return maxError;
}
//==============================================================================
template <typename FloatType>
double LookupTableTransform<FloatType>::calculateRelativeDifference (double x, double y) noexcept
{
static const auto eps = std::numeric_limits<double>::min();
auto absX = std::abs (x);
auto absY = std::abs (y);
auto absDiff = std::abs (x - y);
if (absX < eps)
{
if (absY >= eps)
return absDiff / absY;
return absDiff; // return the absolute error if both numbers are too close to zero
}
return absDiff / std::min (absX, absY);
}
//==============================================================================
template class LookupTable<float>;
template class LookupTable<double>;
template class LookupTableTransform<float>;
template class LookupTableTransform<double>;
} // namespace dsp
} // namespace juce

332
modules/juce_dsp/maths/juce_LookupTable.h Normal file
View File

@ -0,0 +1,332 @@
/*
==============================================================================
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.
==============================================================================
*/
namespace juce
{
namespace dsp
{
/**
Class for efficiently approximating expensive arithmetic operations.
The approximation is based on linear interpolation between pre-calculated values.
The approximated function should be passed as a callable object to the constructor
along with the number of data points to be pre-calculated. The accuracy of the
approximation can be increased by using more points at the cost of a larger memory
footprint.
Consider using LookupTableTransform as an easy-to-use alternative.
Example:
LookupTable<float> lut ([] (size_t i) { return std::sqrt ((float) i); }, 64);
auto outValue = lut[17];
@see LookupTableTransform
@tags{DSP}
*/
template <typename FloatType>
class LookupTable
{
public:
/** Creates an uninitialised LookupTable object.
You need to call initialise() before using the object. Prefer using the
non-default constructor instead.
@see initialise
*/
LookupTable();
/** Creates and initialises a LookupTable object.
@param functionToApproximate The function to be approximated. This should be a
mapping from the integer range [0, numPointsToUse - 1].
@param numPointsToUse The number of pre-calculated values stored.
*/
LookupTable (const std::function<FloatType (size_t)>& functionToApproximate, size_t numPointsToUse);
/** Initialises or changes the parameters of a LookupTable object.
This function can be used to change what function is approximated by an already
constructed LookupTable along with the number of data points used. If the function
to be approximated won't ever change, prefer using the non-default constructor.
@param functionToApproximate The function to be approximated. This should be a
mapping from the integer range [0, numPointsToUse - 1].
@param numPointsToUse The number of pre-calculated values stored.
*/
void initialise (const std::function<FloatType (size_t)>& functionToApproximate, size_t numPointsToUse);
//==============================================================================
/** Calculates the approximated value for the given index without range checking.
Use this if you can guarantee that the index is non-negative and less than numPoints.
Otherwise use get().
@param index The approximation is calculated for this non-integer index.
@return The approximated value at the given index.
@see get, operator[]
*/
FloatType getUnchecked (FloatType index) const noexcept
{
jassert (isInitialised()); // Use the non-default constructor or call initialise() before first use
jassert (isPositiveAndBelow (index, FloatType (getNumPoints())));
auto i = truncatePositiveToUnsignedInt (index);
auto f = index - FloatType (i);
jassert (isPositiveAndBelow (f, FloatType (1)));
auto x0 = data.getUnchecked (static_cast<int> (i));
auto x1 = data.getUnchecked (static_cast<int> (i + 1));
return jmap (f, x0, x1);
}
//==============================================================================
/** Calculates the approximated value for the given index with range checking.
This can be called with any input indices. If the provided index is out-of-range
either the bottom or the top element of the LookupTable is returned.
If the index is guaranteed to be in range use the faster getUnchecked() instead.
@param index The approximation is calculated for this non-integer index.
@return The approximated value at the given index.
@see getUnchecked, operator[]
*/
FloatType get (FloatType index) const noexcept
{
if (index >= getNumPoints())
index = static_cast<FloatType> (getGuardIndex());
else if (index < 0)
index = {};
return getUnchecked (index);
}
//==============================================================================
/** @see getUnchecked */
FloatType operator[] (FloatType index) const noexcept { return getUnchecked (index); }
/** Returns the size of the LookupTable, i.e., the number of pre-calculated data points. */
size_t getNumPoints() const noexcept { return static_cast<size_t> (data.size()) - 1; }
/** Returns true if the LookupTable is initialised and ready to be used. */
bool isInitialised() const noexcept { return data.size() > 1; }
private:
//==============================================================================
Array<FloatType> data;
void prepare() noexcept;
static size_t getRequiredBufferSize (size_t numPointsToUse) noexcept { return numPointsToUse + 1; }
size_t getGuardIndex() const noexcept { return getRequiredBufferSize (getNumPoints()) - 1; }
JUCE_DECLARE_NON_COPYABLE_WITH_LEAK_DETECTOR (LookupTable)
};
//==============================================================================
/** Class for approximating expensive arithmetic operations.
Once initialised, this class can be used just like the function it approximates
via operator().
Example:
LookupTableTransform<float> tanhApprox ([] (float x) { return std::tanh (x); }, -5.0f, 5.0f, 64);
auto outValue = tanhApprox (4.2f);
Note : if you try to call the function with an input outside the provided
range, it will return either the first or the last recorded LookupTable value.
@see LookupTable
@tags{DSP}
*/
template <typename FloatType>
class LookupTableTransform
{
public:
//==============================================================================
/** Creates an uninitialised LookupTableTransform object.
You need to call initialise() before using the object. Prefer using the
non-default constructor instead.
@see initialise
*/
LookupTableTransform()
{}
//==============================================================================
/** Creates and initialises a LookupTableTransform object.
@param functionToApproximate The function to be approximated. This should be a
mapping from a FloatType to FloatType.
@param minInputValueToUse The lowest input value used. The approximation will
fail for values lower than this.
@param maxInputValueToUse The highest input value used. The approximation will
fail for values higher than this.
@param numPoints The number of pre-calculated values stored.
*/
LookupTableTransform (const std::function<FloatType (FloatType)>& functionToApproximate,
FloatType minInputValueToUse,
FloatType maxInputValueToUse,
size_t numPoints)
{
initialise (functionToApproximate, minInputValueToUse, maxInputValueToUse, numPoints);
}
//==============================================================================
/** Initialises or changes the parameters of a LookupTableTransform object.
@param functionToApproximate The function to be approximated. This should be a
mapping from a FloatType to FloatType.
@param minInputValueToUse The lowest input value used. The approximation will
fail for values lower than this.
@param maxInputValueToUse The highest input value used. The approximation will
fail for values higher than this.
@param numPoints The number of pre-calculated values stored.
*/
void initialise (const std::function<FloatType (FloatType)>& functionToApproximate,
FloatType minInputValueToUse,
FloatType maxInputValueToUse,
size_t numPoints);
//==============================================================================
/** Calculates the approximated value for the given input value without range checking.
Use this if you can guarantee that the input value is within the range specified
in the constructor or initialise(), otherwise use processSample().
@param value The approximation is calculated for this input value.
@return The approximated value for the provided input value.
@see processSample, operator(), operator[]
*/
FloatType processSampleUnchecked (FloatType value) const noexcept
{
jassert (value >= minInputValue && value <= maxInputValue);
return lookupTable[scaler * value + offset];
}
//==============================================================================
/** Calculates the approximated value for the given input value with range checking.
This can be called with any input values. Out-of-range input values will be
clipped to the specified input range.
If the index is guaranteed to be in range use the faster processSampleUnchecked()
instead.
@param value The approximation is calculated for this input value.
@return The approximated value for the provided input value.
@see processSampleUnchecked, operator(), operator[]
*/
FloatType processSample (FloatType value) const noexcept
{
auto index = scaler * jlimit (minInputValue, maxInputValue, value) + offset;
jassert (isPositiveAndBelow (index, FloatType (lookupTable.getNumPoints())));
return lookupTable[index];
}
//==============================================================================
/** @see processSampleUnchecked */
FloatType operator[] (FloatType index) const noexcept { return processSampleUnchecked (index); }
/** @see processSample */
FloatType operator() (FloatType index) const noexcept { return processSample (index); }
//==============================================================================
/** Processes an array of input values without range checking
@see process
*/
void processUnchecked (const FloatType* input, FloatType* output, size_t numSamples) const noexcept
{
for (size_t i = 0; i < numSamples; ++i)
output[i] = processSampleUnchecked (input[i]);
}
//==============================================================================
/** Processes an array of input values with range checking
@see processUnchecked
*/
void process (const FloatType* input, FloatType* output, size_t numSamples) const noexcept
{
for (size_t i = 0; i < numSamples; ++i)
output[i] = processSample (input[i]);
}
//==============================================================================
/** Calculates the maximum relative error of the approximation for the specified
parameter set.
The closer the returned value is to zero the more accurate the approximation
is.
This function compares the approximated output of this class to the function
it approximates at a range of points and returns the maximum relative error.
This can be used to determine if the approximation is suitable for the given
problem. The accuracy of the approximation can generally be improved by
increasing numPoints.
@param functionToApproximate The approximated function. This should be a
mapping from a FloatType to FloatType.
@param minInputValue The lowest input value used.
@param maxInputValue The highest input value used.
@param numPoints The number of pre-calculated values stored.
@param numTestPoints The number of input values used for error
calculation. Higher numbers can increase the
accuracy of the error calculation. If it's zero
then 100 * numPoints will be used.
*/
static double calculateMaxRelativeError (const std::function<FloatType (FloatType)>& functionToApproximate,
FloatType minInputValue,
FloatType maxInputValue,
size_t numPoints,
size_t numTestPoints = 0);
private:
//==============================================================================
static double calculateRelativeDifference (double, double) noexcept;
//==============================================================================
LookupTable<FloatType> lookupTable;
FloatType minInputValue, maxInputValue;
FloatType scaler, offset;
JUCE_DECLARE_NON_COPYABLE_WITH_LEAK_DETECTOR (LookupTableTransform)
};
} // namespace dsp
} // namespace juce

318
modules/juce_dsp/maths/juce_Matrix.cpp Normal file
View File

@ -0,0 +1,318 @@
/*
==============================================================================
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.
==============================================================================
*/
namespace juce
{
namespace dsp
{
template <typename ElementType>
Matrix<ElementType> Matrix<ElementType>::identity (size_t size)
{
Matrix result (size, size);
for (size_t i = 0; i < size; ++i)
result(i, i) = 1;
return result;
}
template <typename ElementType>
Matrix<ElementType> Matrix<ElementType>::toeplitz (const Matrix& vector, size_t size)
{
jassert (vector.isOneColumnVector());
jassert (size <= vector.rows);
Matrix result (size, size);
for (size_t i = 0; i < size; ++i)
result (i, i) = vector (0, 0);
for (size_t i = 1; i < size; ++i)
for (size_t j = i; j < size; ++j)
result (j, j - i) = result (j - i, j) = vector (i, 0);
return result;
}
template <typename ElementType>
Matrix<ElementType> Matrix<ElementType>::hankel (const Matrix& vector, size_t size, size_t offset)
{
jassert(vector.isOneColumnVector());
jassert(vector.rows >= (2 * (size - 1) + 1));
Matrix result (size, size);
for (size_t i = 0; i < size; ++i)
result (i, i) = vector ((2 * i) + offset, 0);
for (size_t i = 1; i < size; ++i)
for (size_t j = i; j < size; ++j)
result (j, j - i) = result (j - i, j) = vector (i + 2 * (j - i) + offset, 0);
return result;
}
//==============================================================================
template <typename ElementType>
Matrix<ElementType>& Matrix<ElementType>::swapColumns (size_t columnOne, size_t columnTwo) noexcept
{
jassert (columnOne < columns && columnTwo < columns);
auto* p = data.getRawDataPointer();
for (size_t i = 0; i < rows; ++i)
{
auto offset = dataAcceleration.getUnchecked (static_cast<int> (i));
std::swap (p[offset + columnOne], p[offset + columnTwo]);
}
return *this;
}
template <typename ElementType>
Matrix<ElementType>& Matrix<ElementType>::swapRows (size_t rowOne, size_t rowTwo) noexcept
{
jassert (rowOne < rows && rowTwo < rows);
auto offset1 = rowOne * columns;
auto offset2 = rowTwo * columns;
auto* p = data.getRawDataPointer();
for (size_t i = 0; i < columns; ++i)
std::swap (p[offset1 + i], p[offset2 + i]);
return *this;
}
//==============================================================================
template <typename ElementType>
Matrix<ElementType> Matrix<ElementType>::operator* (const Matrix<ElementType>& other) const
{
auto n = getNumRows(), m = other.getNumColumns(), p = getNumColumns();
Matrix result (n, m);
jassert (p == other.getNumRows());
size_t offsetMat = 0, offsetlhs = 0;
auto* dst = result.getRawDataPointer();
auto* a = getRawDataPointer();
auto* b = other.getRawDataPointer();
for (size_t i = 0; i < n; ++i)
{
size_t offsetrhs = 0;
for (size_t k = 0; k < p; ++k)
{
auto ak = a[offsetlhs++];
for (size_t j = 0; j < m; ++j)
dst[offsetMat + j] += ak * b[offsetrhs + j];
offsetrhs += m;
}
offsetMat += m;
}
return result;
}
//==============================================================================
template <typename ElementType>
bool Matrix<ElementType>::compare (const Matrix& a, const Matrix& b, ElementType tolerance) noexcept
{
if (a.rows != b.rows || a.columns != b.columns)
return false;
tolerance = std::abs (tolerance);
auto* bPtr = b.begin();
for (auto aValue : a)
if (std::abs (aValue - *bPtr++) > tolerance)
return false;
return true;
}
//==============================================================================
template <typename ElementType>
bool Matrix<ElementType>::solve (Matrix& b) const noexcept
{
auto n = columns;
jassert (n == n && n == b.rows && b.isOneColumnVector());
auto* x = b.getRawDataPointer();
const auto& A = *this;
switch (n)
{
case 1:
{
auto denominator = A (0,0);
if (denominator == 0)
return false;
b (0, 0) /= denominator;
}
break;
case 2:
{
auto denominator = A (0, 0) * A (1, 1) - A (0, 1) * A (1, 0);
if (denominator == 0)
return false;
auto factor = (1 / denominator);
auto b0 = x[0], b1 = x[1];
x[0] = factor * (A (1, 1) * b0 - A (0, 1) * b1);
x[1] = factor * (A (0, 0) * b1 - A (1, 0) * b0);
}
break;
case 3:
{
auto denominator = A (0, 0) * (A (1, 1) * A (2, 2) - A (1, 2) * A (2, 1))
+ A (0, 1) * (A (1, 2) * A (2, 0) - A (1, 0) * A (2, 2))
+ A (0, 2) * (A (1, 0) * A (2, 1) - A (1, 1) * A (2, 0));
if (denominator == 0)
return false;
auto factor = 1 / denominator;
auto b0 = x[0], b1 = x[1], b2 = x[2];
x[0] = ( ( A (0, 1) * A (1, 2) - A (0, 2) * A (1, 1)) * b2
+ (-A (0, 1) * A (2, 2) + A (0, 2) * A (2, 1)) * b1
+ ( A (1, 1) * A (2, 2) - A (1, 2) * A (2, 1)) * b0) * factor;
x[1] = -( ( A (0, 0) * A (1, 2) - A (0, 2) * A (1, 0)) * b2
+ (-A (0, 0) * A (2, 2) + A (0, 2) * A (2, 0)) * b1
+ ( A (1, 0) * A (2, 2) - A (1, 2) * A (2, 0)) * b0) * factor;
x[2] = ( ( A (0, 0) * A (1, 1) - A (0, 1) * A (1, 0)) * b2
+ (-A (0, 0) * A (2, 1) + A (0, 1) * A (2, 0)) * b1
+ ( A (1, 0) * A (2, 1) - A (1, 1) * A (2, 0)) * b0) * factor;
}
break;
default:
{
Matrix<ElementType> M (A);
for (size_t j = 0; j < n; ++j)
{
if (M (j, j) == 0)
{
auto i = j;
while (i < n && M (i, j) == 0)
++i;
if (i == n)
return false;
for (size_t k = 0; k < n; ++k)
M (j, k) += M (i, k);
x[j] += x[i];
}
auto t = 1 / M (j, j);
for (size_t k = 0; k < n; ++k)
M (j, k) *= t;
x[j] *= t;
for (size_t k = j + 1; k < n; ++k)
{
auto u = -M (k, j);
for (size_t l = 0; l < n; ++l)
M (k, l) += u * M (j, l);
x[k] += u * x[j];
}
}
for (int k = static_cast<int> (n) - 2; k >= 0; --k)
for (size_t i = static_cast<size_t> (k) + 1; i < n; ++i)
x[k] -= M (static_cast<size_t> (k), i) * x[i];
}
}
return true;
}
//==============================================================================
template <typename ElementType>
String Matrix<ElementType>::toString() const
{
StringArray entries;
int sizeMax = 0;
auto* p = data.begin();
for (size_t i = 0; i < rows; ++i)
{
for (size_t j = 0; j < columns; ++j)
{
String entry (*p++, 4);
sizeMax = jmax (sizeMax, entry.length());
entries.add (entry);
}
}
sizeMax = ((sizeMax + 1) / 4 + 1) * 4;
MemoryOutputStream result;
auto n = static_cast<size_t> (entries.size());
for (size_t i = 0; i < n; ++i)
{
result << entries[(int) i].paddedRight (' ', sizeMax);
if (i % columns == (columns - 1))
result << newLine;
}
return result.toString();
}
template class Matrix<float>;
template class Matrix<double>;
} // namespace dsp
} // namespace juce

255
modules/juce_dsp/maths/juce_Matrix.h Normal file
View File

@ -0,0 +1,255 @@
/*
==============================================================================
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.
==============================================================================
*/
namespace juce
{
namespace dsp
{
/**
General matrix and vectors class, meant for classic math manipulation such as
additions, multiplications, and linear systems of equations solving.
@see LinearAlgebra
@tags{DSP}
*/
template<typename ElementType>
class Matrix
{
public:
//==============================================================================
/** Creates a new matrix with a given number of rows and columns. */
Matrix (size_t numRows, size_t numColumns)
: rows (numRows), columns (numColumns)
{
resize();
clear();
}
/** Creates a new matrix with a given number of rows and columns, with initial
data coming from an array, stored in row-major order.
*/
Matrix (size_t numRows, size_t numColumns, const ElementType* dataPointer)
: rows (numRows), columns (numColumns)
{
resize();
memcpy (data.getRawDataPointer(), dataPointer, rows * columns * sizeof (ElementType));
}
/** Creates a copy of another matrix. */
Matrix (const Matrix&) = default;
/** Moves a copy of another matrix. */
Matrix (Matrix&&) noexcept = default;
/** Creates a copy of another matrix. */
Matrix& operator= (const Matrix&) = default;
/** Moves another matrix into this one */
Matrix& operator= (Matrix&&) noexcept = default;
//==============================================================================
/** Creates the identity matrix */
static Matrix identity (size_t size);
/** Creates a Toeplitz Matrix from a vector with a given squared size */
static Matrix toeplitz (const Matrix& vector, size_t size);
/** Creates a squared size x size Hankel Matrix from a vector with an optional offset.
@param vector The vector from which the Hankel matrix should be generated.
Its number of rows should be at least 2 * (size - 1) + 1
@param size The size of resulting square matrix.
@param offset An optional offset into the given vector.
*/
static Matrix hankel (const Matrix& vector, size_t size, size_t offset = 0);
//==============================================================================
/** Returns the number of rows in the matrix. */
size_t getNumRows() const noexcept { return rows; }
/** Returns the number of columns in the matrix. */
size_t getNumColumns() const noexcept { return columns; }
/** Returns an Array of 2 integers with the number of rows and columns in the
matrix.
*/
Array<size_t> getSize() const noexcept { return { rows, columns }; }
/** Fills the contents of the matrix with zeroes. */
void clear() noexcept { zeromem (data.begin(), sizeof (ElementType) * (size_t) data.size()); }
//==============================================================================
/** Swaps the contents of two rows in the matrix and returns a reference to itself. */
Matrix& swapRows (size_t rowOne, size_t rowTwo) noexcept;
/** Swaps the contents of two columns in the matrix and returns a reference to itself. */
Matrix& swapColumns (size_t columnOne, size_t columnTwo) noexcept;
//==============================================================================
/** Returns the value of the matrix at a given row and column (for reading). */
inline ElementType operator() (size_t row, size_t column) const noexcept
{
jassert (row < rows && column < columns);
return data.getReference (static_cast<int> (dataAcceleration.getReference (static_cast<int> (row))) + static_cast<int> (column));
}
/** Returns the value of the matrix at a given row and column (for modifying). */
inline ElementType& operator() (size_t row, size_t column) noexcept
{
jassert (row < rows && column < columns);
return data.getReference (static_cast<int> (dataAcceleration.getReference (static_cast<int> (row))) + static_cast<int> (column));
}
/** Returns a pointer to the raw data of the matrix object, ordered in row-major
order (for modifying).
*/
inline ElementType* getRawDataPointer() noexcept { return data.getRawDataPointer(); }
/** Returns a pointer to the raw data of the matrix object, ordered in row-major
order (for reading).
*/
inline const ElementType* getRawDataPointer() const noexcept { return data.begin(); }
//==============================================================================
/** Addition of two matrices */
inline Matrix& operator+= (const Matrix& other) noexcept { return apply (other, [] (ElementType a, ElementType b) { return a + b; } ); }
/** Subtraction of two matrices */
inline Matrix& operator-= (const Matrix& other) noexcept { return apply (other, [] (ElementType a, ElementType b) { return a - b; } ); }
/** Scalar multiplication */
inline Matrix& operator*= (ElementType scalar) noexcept
{
std::for_each (begin(), end(), [scalar] (ElementType& x) { x *= scalar; });
return *this;
}
/** Addition of two matrices */
inline Matrix operator+ (const Matrix& other) const { Matrix result (*this); result += other; return result; }
/** Addition of two matrices */
inline Matrix operator- (const Matrix& other) const { Matrix result (*this); result -= other; return result; }
/** Scalar multiplication */
inline Matrix operator* (ElementType scalar) const { Matrix result (*this); result *= scalar; return result; }
/** Matrix multiplication */
Matrix operator* (const Matrix& other) const;
/** Does a hadarmard product with the receiver and other and stores the result in the receiver */
inline Matrix& hadarmard (const Matrix& other) noexcept { return apply (other, [] (ElementType a, ElementType b) { return a * b; } ); }
/** Does a hadarmard product with a and b returns the result. */
static inline Matrix hadarmard (const Matrix& a, const Matrix& b) { Matrix result (a); result.hadarmard (b); return result; }
//==============================================================================
/** Compare to matrices with a given tolerance */
static bool compare (const Matrix& a, const Matrix& b, ElementType tolerance = 0) noexcept;
/* Comparison operator */
inline bool operator== (const Matrix& other) const noexcept { return compare (*this, other); }
//==============================================================================
/** Tells if the matrix is a square matrix */
bool isSquare() const noexcept { return rows == columns; }
/** Tells if the matrix is a vector */
bool isVector() const noexcept { return isOneColumnVector() || isOneRowVector(); }
/** Tells if the matrix is a one column vector */
bool isOneColumnVector() const noexcept { return columns == 1; }
/** Tells if the matrix is a one row vector */
bool isOneRowVector() const noexcept { return rows == 1; }
/** Tells if the matrix is a null matrix */
bool isNullMatrix() const noexcept { return rows == 0 || columns == 0; }
//==============================================================================
/** Solves a linear system of equations represented by this object and the argument b,
using various algorithms depending on the size of the arguments.
The matrix must be a square matrix N times N, and b must be a vector N times 1,
with the coefficients of b. After the execution of the algorithm,
the vector b will contain the solution.
Returns true if the linear system of euqations was successfully solved.
*/
bool solve (Matrix& b) const noexcept;
//==============================================================================
/** Returns a String displaying in a convenient way the matrix contents. */
String toString() const;
//==============================================================================
ElementType* begin() noexcept { return data.begin(); }
ElementType* end() noexcept { return data.end(); }
const ElementType* begin() const noexcept { return &data.getReference (0); }
const ElementType* end() const noexcept { return begin() + data.size(); }
private:
//==============================================================================
/** Resizes the matrix. */
void resize()
{
data.resize (static_cast<int> (columns * rows));
dataAcceleration.resize (static_cast<int> (rows));
for (size_t i = 0; i < rows; ++i)
dataAcceleration.setUnchecked (static_cast<int> (i), i * columns);
}
template <typename BinaryOperation>
Matrix& apply (const Matrix& other, BinaryOperation binaryOp)
{
jassert (rows == other.rows && columns == other.columns);
auto* dst = getRawDataPointer();
for (auto src : other)
{
*dst = binaryOp (*dst, src);
++dst;
}
return *this;
}
//==============================================================================
Array<ElementType> data;
Array<size_t> dataAcceleration;
size_t rows, columns;
//==============================================================================
JUCE_LEAK_DETECTOR (Matrix)
};
} // namespace dsp
} // namespace juce

172
modules/juce_dsp/maths/juce_Matrix_test.cpp Normal file
View File

@ -0,0 +1,172 @@
/*
==============================================================================
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.
==============================================================================
*/
namespace juce
{
namespace dsp
{
struct LinearAlgebraUnitTest : public UnitTest
{
LinearAlgebraUnitTest() : UnitTest ("Linear Algebra UnitTests", "DSP") {}
struct AdditionTest
{
template <typename ElementType>
static void run (LinearAlgebraUnitTest& u)
{
const ElementType data1[] = { 1, 2, 3, 4, 5, 6, 7, 8 };
const ElementType data2[] = { 1, -1, 3, -1, 5, -1, 7, -1 };
const ElementType data3[] = { 2, 1, 6, 3, 10, 5, 14, 7 };
Matrix<ElementType> mat1 (2, 4, data1);
Matrix<ElementType> mat2 (2, 4, data2);
Matrix<ElementType> mat3 (2, 4, data3);
u.expect((mat1 + mat2) == mat3);
}
};
struct DifferenceTest
{
template <typename ElementType>
static void run (LinearAlgebraUnitTest& u)
{
const ElementType data1[] = { 1, 2, 3, 4, 5, 6, 7, 8 };
const ElementType data2[] = { 1, -1, 3, -1, 5, -1, 7, -1 };
const ElementType data3[] = { 0, 3, 0, 5, 0, 7, 0, 9 };
Matrix<ElementType> mat1 (2, 4, data1);
Matrix<ElementType> mat2 (2, 4, data2);
Matrix<ElementType> mat3 (2, 4, data3);
u.expect((mat1 - mat2) == mat3);
}
};
struct ScalarMultiplicationTest
{
template <typename ElementType>
static void run (LinearAlgebraUnitTest& u)
{
const ElementType data1[] = { 1, 2, 3, 4, 5, 6, 7, 8 };
const ElementType scalar = 2.0;
const ElementType data2[] = { 2, 4, 6, 8, 10, 12, 14, 16 };
Matrix<ElementType> x (2, 4, data1);
Matrix<ElementType> expected (2, 4, data2);
u.expect ((x * scalar) == expected);
}
};
struct HadamardProductTest
{
template <typename ElementType>
static void run (LinearAlgebraUnitTest& u)
{
const ElementType data1[] = { 1, 2, 3, 4, 5, 6, 7, 8 };
const ElementType data2[] = { 1, -1, 3, -1, 5, -1, 7, -1 };
const ElementType data3[] = { 1, -2, 9, -4, 25, -6, 49, -8 };
Matrix<ElementType> mat1 (2, 4, data1);
Matrix<ElementType> mat2 (2, 4, data2);
Matrix<ElementType> mat3 (2, 4, data3);
u.expect (Matrix<ElementType>::hadarmard (mat1, mat2) == mat3);
}
};
struct MultiplicationTest
{
template <typename ElementType>
static void run (LinearAlgebraUnitTest& u)
{
const ElementType data1[] = { 1, 2, 3, 4, 5, 6, 7, 8 };
const ElementType data2[] = { 1, -1, 3, -1, 5, -1, 7, -1 };
const ElementType data3[] = { 50, -10, 114, -26 };
Matrix<ElementType> mat1 (2, 4, data1);
Matrix<ElementType> mat2 (4, 2, data2);
Matrix<ElementType> mat3 (2, 2, data3);
u.expect((mat1 * mat2) == mat3);
}
};
struct IdentityMatrixTest
{
template <typename ElementType>
static void run (LinearAlgebraUnitTest& u)
{
const ElementType data1[] = { 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1};
u.expect (Matrix<ElementType>::identity (4) == Matrix<ElementType> (4, 4, data1));
}
};
struct SolvingTest
{
template <typename ElementType>
static void run (LinearAlgebraUnitTest& u)
{
const ElementType data1[] = { 1, -1, 2, -2 };
const ElementType data2[] = { -1, 0, -1, -7 };
const ElementType data3[] = { 1, 4, 2, 1, -1, 1, 4, 3, -2, -1, 1, 1, -1, 0, 1, 4 };
Matrix<ElementType> X (4, 1, data1);
Matrix<ElementType> B (4, 1, data2);
Matrix<ElementType> A (4, 4, data3);
u.expect (A.solve (B));
u.expect (Matrix<ElementType>::compare (X, B, (ElementType) 1e-4));
}
};
template <class TheTest>
void runTestForAllTypes (const char* unitTestName)
{
beginTest (unitTestName);
TheTest::template run<float> (*this);
TheTest::template run<double> (*this);
}
void runTest() override
{
runTestForAllTypes<AdditionTest> ("AdditionTest");
runTestForAllTypes<DifferenceTest> ("DifferenceTest");
runTestForAllTypes<ScalarMultiplicationTest> ("ScalarMultiplication");
runTestForAllTypes<HadamardProductTest> ("HadamardProductTest");
runTestForAllTypes<MultiplicationTest> ("MultiplicationTest");
runTestForAllTypes<IdentityMatrixTest> ("IdentityMatrixTest");
runTestForAllTypes<SolvingTest> ("SolvingTest");
}
};
static LinearAlgebraUnitTest linearAlgebraUnitTest;
} // namespace dsp
} // namespace juce

68
modules/juce_dsp/maths/juce_Phase.h Normal file
View File

@ -0,0 +1,68 @@
/*
==============================================================================
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.
==============================================================================
*/
namespace juce
{
namespace dsp
{
/**
Represents an increasing phase value between 0 and 2*pi.
This represents a value which can be incremented, and which wraps back to 0 when it
goes past 2 * pi.
@tags{DSP}
*/
template <typename Type>
struct Phase
{
/** Resets the phase to 0. */
void reset() noexcept { phase = 0; }
/** Returns the current value, and increments the phase by the given increment.
The increment must be a positive value, it can't go backwards!
The new value of the phase after calling this function will be (phase + increment) % (2 * pi).
*/
Type advance (Type increment) noexcept
{
jassert (increment >= 0); // cannot run this value backwards!
auto last = phase;
auto next = last + increment;
while (next >= MathConstants<Type>::twoPi)
next -= MathConstants<Type>::twoPi;
phase = next;
return last;
}
Type phase = 0;
};
} // namespace dsp
} // namespace juce

169
modules/juce_dsp/maths/juce_Polynomial.h Normal file
View File

@ -0,0 +1,169 @@
/*
==============================================================================
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.
==============================================================================
*/
namespace juce
{
namespace dsp
{
/**
A class representing a polynomial
@tags{DSP}
*/
template <typename FloatingType>
class Polynomial
{
public:
//==============================================================================
/** Creates a new polynomial which will always evaluate to zero. */
Polynomial()
{
coeffs.add (0);
}
/** Creates a new polynomial with given coefficients.
@param numCoefficients The number of coefficients stored in coefficients.
This is also the order of the returned polynomial.
@param coefficients The coefficients which will be used by the newly
created polynomial. The Polynomial class will keep
a private copy of the coefficients.
*/
Polynomial (const FloatingType* coefficients, int numCoefficients)
: coeffs (coefficients, numCoefficients)
{
jassert (! coeffs.isEmpty());
}
/** Creates a copy of another polynomial. */
Polynomial (const Polynomial&) = default;
/** Creates a copy of another polynomial. */
Polynomial (Polynomial&&) = default;
/** Creates a copy of another polynomial. */
Polynomial& operator= (const Polynomial&) = default;
/** Creates a copy of another polynomial. */
Polynomial& operator= (Polynomial&&) = default;
/** Creates a new polynomial with coefficients by a C++11 initializer list.
This function can be used in the following way:
Polynomial<float> p ({0.5f, -0.3f, 0.2f});
*/
template <typename... Values>
Polynomial (Values... items) : coeffs (items...)
{
jassert (! coeffs.isEmpty());
}
//==============================================================================
/** Returns a single coefficient of the receiver for reading */
FloatingType operator[] (int index) const noexcept { return coeffs.getUnchecked (index); }
/** Returns a single coefficient of the receiver for modifying. */
FloatingType& operator[] (int index) noexcept { return coeffs.getReference (index); }
/** Evaluates the value of the polynomial at a single point x. */
FloatingType operator() (FloatingType x) const noexcept
{
// Horner's method
FloatingType y (0);
for (int i = coeffs.size(); --i >= 0;)
y = (x * y) + coeffs.getUnchecked(i);
return y;
}
/** Returns the order of the polynomial. */
int getOrder() noexcept
{
return coeffs.size() - 1;
}
//==============================================================================
/** Returns the polynomial with all its coefficients multiplied with a gain factor */
Polynomial<FloatingType> withGain (double gain) const
{
auto result = *this;
for (auto& c : result.coeffs)
c *= gain;
return result;
}
/** Returns the sum of this polynomial with another */
Polynomial<FloatingType> getSumWith (const Polynomial<FloatingType>& other) const
{
if (coeffs.size() < other.coeffs.size())
return other.getSumWith (*this);
auto result = *this;
for (int i = 0; i < other.coeffs.size(); ++i)
result[i] += other[i];
return result;
}
/** computes the product of two polynomials and return the result */
Polynomial<FloatingType> getProductWith (const Polynomial<FloatingType>& other) const
{
Polynomial<FloatingType> result;
result.coeffs.clearQuick();
auto N1 = coeffs.size();
auto N2 = other.coeffs.size();
auto Nmax = jmax (N1, N2);
auto N = N1 + N2 - 1;
for (int i = 0; i < N; ++i)
{
FloatingType value (0);
for (int j = 0; j < Nmax; ++j)
if (j >= 0 && j < N1 && i - j >= 0 && i - j < N2)
value = value + (*this)[j] * other[i - j];
result.coeffs.add (value);
}
return result;
}
private:
//==============================================================================
Array<FloatingType> coeffs;
JUCE_LEAK_DETECTOR (Polynomial)
};
} // namespace dsp
} // namespace juce

144
modules/juce_dsp/maths/juce_SpecialFunctions.cpp Normal file
View File

@ -0,0 +1,144 @@
/*
==============================================================================
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.
==============================================================================
*/
namespace juce
{
namespace dsp
{
double SpecialFunctions::besselI0 (double x) noexcept
{
auto ax = std::abs (x);
if (ax < 3.75)
{
auto y = x / 3.75;
y *= y;
return 1.0 + y * (3.5156229 + y * (3.0899424 + y * (1.2067492
+ y * (0.2659732 + y * (0.360768e-1 + y * 0.45813e-2)))));
}
auto y = 3.75 / ax;
return (std::exp (ax) / std::sqrt (ax))
* (0.39894228 + y * (0.1328592e-1 + y * (0.225319e-2 + y * (-0.157565e-2 + y * (0.916281e-2
+ y * (-0.2057706e-1 + y * (0.2635537e-1 + y * (-0.1647633e-1 + y * 0.392377e-2))))))));
}
void SpecialFunctions::ellipticIntegralK (double k, double& K, double& Kp) noexcept
{
constexpr int M = 4;
K = MathConstants<double>::halfPi;
auto lastK = k;
for (int i = 0; i < M; ++i)
{
lastK = std::pow (lastK / (1 + std::sqrt (1 - std::pow (lastK, 2.0))), 2.0);
K *= 1 + lastK;
}
Kp = MathConstants<double>::halfPi;
auto last = std::sqrt (1 - k * k);
for (int i = 0; i < M; ++i)
{
last = std::pow (last / (1.0 + std::sqrt (1.0 - std::pow (last, 2.0))), 2.0);
Kp *= 1 + last;
}
}
Complex<double> SpecialFunctions::cde (Complex<double> u, double k) noexcept
{
constexpr int M = 4;
double ke[M + 1];
double* kei = ke;
*kei = k;
for (int i = 0; i < M; ++i)
{
auto next = std::pow (*kei / (1.0 + std::sqrt (1.0 - std::pow (*kei, 2.0))), 2.0);
*++kei = next;
}
// NB: the spurious cast to double here is a workaround for a very odd link-time failure
std::complex<double> last = std::cos (u * (double) MathConstants<double>::halfPi);
for (int i = M - 1; i >= 0; --i)
last = (1.0 + ke[i + 1]) / (1.0 / last + ke[i + 1] * last);
return last;
}
Complex<double> SpecialFunctions::sne (Complex<double> u, double k) noexcept
{
constexpr int M = 4;
double ke[M + 1];
double* kei = ke;
*kei = k;
for (int i = 0; i < M; ++i)
{
auto next = std::pow (*kei / (1 + std::sqrt (1 - std::pow (*kei, 2.0))), 2.0);
*++kei = next;
}
// NB: the spurious cast to double here is a workaround for a very odd link-time failure
std::complex<double> last = std::sin (u * (double) MathConstants<double>::halfPi);
for (int i = M - 1; i >= 0; --i)
last = (1.0 + ke[i + 1]) / (1.0 / last + ke[i + 1] * last);
return last;
}
Complex<double> SpecialFunctions::asne (Complex<double> w, double k) noexcept
{
constexpr int M = 4;
double ke[M + 1];
double* kei = ke;
*kei = k;
for (int i = 0; i < M; ++i)
{
auto next = std::pow (*kei / (1.0 + std::sqrt (1.0 - std::pow (*kei, 2.0))), 2.0);
*++kei = next;
}
std::complex<double> last = w;
for (int i = 1; i <= M; ++i)
last = 2.0 * last / ((1.0 + ke[i]) * (1.0 + std::sqrt (1.0 - std::pow (ke[i - 1] * last, 2.0))));
return 2.0 / MathConstants<double>::pi * std::asin (last);
}
} // namespace dsp
} // namespace juce

68
modules/juce_dsp/maths/juce_SpecialFunctions.h Normal file
View File

@ -0,0 +1,68 @@
/*
==============================================================================
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.
==============================================================================
*/
namespace juce
{
namespace dsp
{
/**
Contains miscellaneous filter design and windowing functions.
@tags{DSP}
*/
struct SpecialFunctions
{
/** Computes the modified Bessel function of the first kind I0 for a
given double value x. Modified Bessel functions are useful to solve
various mathematical problems involving differential equations.
*/
static double besselI0 (double x) noexcept;
/** Computes the complete elliptic integral of the first kind K for a
given double value k, and the associated complete elliptic integral
of the first kind Kp for the complementary modulus of k.
*/
static void ellipticIntegralK (double k, double& K, double& Kp) noexcept;
/** Computes the Jacobian elliptic function cd for the elliptic
modulus k and the quarter-period units u.
*/
static Complex<double> cde (Complex<double> u, double k) noexcept;
/** Computes the Jacobian elliptic function sn for the elliptic
modulus k and the quarter-period units u.
*/
static Complex<double> sne (Complex<double> u, double k) noexcept;
/** Computes the inverse of the Jacobian elliptic function sn
for the elliptic modulus k and the quarter-period units u.
*/
static Complex<double> asne (Complex<double> w, double k) noexcept;
};
} // namespace dsp
} // namespace juce