diff options
Diffstat (limited to 'audio/softsynth/mt32/srchelper/srctools/src')
5 files changed, 672 insertions, 0 deletions
diff --git a/audio/softsynth/mt32/srchelper/srctools/src/FIRResampler.cpp b/audio/softsynth/mt32/srchelper/srctools/src/FIRResampler.cpp new file mode 100644 index 0000000000..2cded0c3d7 --- /dev/null +++ b/audio/softsynth/mt32/srchelper/srctools/src/FIRResampler.cpp @@ -0,0 +1,107 @@ +/* Copyright (C) 2015-2017 Sergey V. Mikayev + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 2.1 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <cmath> +#include <cstring> + +#include "FIRResampler.h" + +using namespace SRCTools; + +FIRResampler::Constants::Constants(const unsigned int upsampleFactor, const double downsampleFactor, const FIRCoefficient kernel[], const unsigned int kernelLength) { + usePhaseInterpolation = downsampleFactor != floor(downsampleFactor); + FIRCoefficient *kernelCopy = new FIRCoefficient[kernelLength]; + memcpy(kernelCopy, kernel, kernelLength * sizeof(FIRCoefficient)); + taps = kernelCopy; + numberOfTaps = kernelLength; + numberOfPhases = upsampleFactor; + phaseIncrement = downsampleFactor; + unsigned int minDelayLineLength = static_cast<unsigned int>(ceil(double(kernelLength) / upsampleFactor)); + unsigned int delayLineLength = 2; + while (delayLineLength < minDelayLineLength) delayLineLength <<= 1; + delayLineMask = delayLineLength - 1; + ringBuffer = new FloatSample[delayLineLength][FIR_INTERPOLATOR_CHANNEL_COUNT]; + FloatSample *s = *ringBuffer; + FloatSample *e = ringBuffer[delayLineLength]; + while (s < e) *(s++) = 0; +} + +FIRResampler::FIRResampler(const unsigned int upsampleFactor, const double downsampleFactor, const FIRCoefficient kernel[], const unsigned int kernelLength) : + constants(upsampleFactor, downsampleFactor, kernel, kernelLength), + ringBufferPosition(0), + phase(constants.numberOfPhases) +{} + +FIRResampler::~FIRResampler() { + delete[] constants.ringBuffer; + delete[] constants.taps; +} + +void FIRResampler::process(const FloatSample *&inSamples, unsigned int &inLength, FloatSample *&outSamples, unsigned int &outLength) { + while (outLength > 0) { + while (needNextInSample()) { + if (inLength == 0) return; + addInSamples(inSamples); + --inLength; + } + getOutSamplesStereo(outSamples); + --outLength; + } +} + +unsigned int FIRResampler::estimateInLength(const unsigned int outLength) const { + return static_cast<unsigned int>((outLength * constants.phaseIncrement + phase) / constants.numberOfPhases); +} + +bool FIRResampler::needNextInSample() const { + return constants.numberOfPhases <= phase; +} + +void FIRResampler::addInSamples(const FloatSample *&inSamples) { + ringBufferPosition = (ringBufferPosition - 1) & constants.delayLineMask; + for (unsigned int i = 0; i < FIR_INTERPOLATOR_CHANNEL_COUNT; i++) { + constants.ringBuffer[ringBufferPosition][i] = *(inSamples++); + } + phase -= constants.numberOfPhases; +} + +// Optimised for processing stereo interleaved streams +void FIRResampler::getOutSamplesStereo(FloatSample *&outSamples) { + FloatSample leftSample = 0.0; + FloatSample rightSample = 0.0; + unsigned int delaySampleIx = ringBufferPosition; + if (constants.usePhaseInterpolation) { + double phaseFraction = phase - floor(phase); + unsigned int maxTapIx = phaseFraction == 0 ? constants.numberOfTaps : constants.numberOfTaps - 1; + for (unsigned int tapIx = static_cast<unsigned int>(phase); tapIx < maxTapIx; tapIx += constants.numberOfPhases) { + FIRCoefficient tap = FIRCoefficient(constants.taps[tapIx] + (constants.taps[tapIx + 1] - constants.taps[tapIx]) * phaseFraction); + leftSample += tap * constants.ringBuffer[delaySampleIx][0]; + rightSample += tap * constants.ringBuffer[delaySampleIx][1]; + delaySampleIx = (delaySampleIx + 1) & constants.delayLineMask; + } + } else { + // Optimised for rational resampling ratios when phase is always integer + for (unsigned int tapIx = static_cast<unsigned int>(phase); tapIx < constants.numberOfTaps; tapIx += constants.numberOfPhases) { + FIRCoefficient tap = constants.taps[tapIx]; + leftSample += tap * constants.ringBuffer[delaySampleIx][0]; + rightSample += tap * constants.ringBuffer[delaySampleIx][1]; + delaySampleIx = (delaySampleIx + 1) & constants.delayLineMask; + } + } + *(outSamples++) = leftSample; + *(outSamples++) = rightSample; + phase += constants.phaseIncrement; +} diff --git a/audio/softsynth/mt32/srchelper/srctools/src/IIR2xResampler.cpp b/audio/softsynth/mt32/srchelper/srctools/src/IIR2xResampler.cpp new file mode 100644 index 0000000000..061006a1e7 --- /dev/null +++ b/audio/softsynth/mt32/srchelper/srctools/src/IIR2xResampler.cpp @@ -0,0 +1,229 @@ +/* Copyright (C) 2015-2017 Sergey V. Mikayev + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 2.1 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <cstddef> + +#include "IIR2xResampler.h" + +namespace SRCTools { + + // Avoid denormals degrading performance, using biased input + static const BufferedSample BIAS = 1e-35f; + + // Sharp elliptic filter with symmetric ripple: N=18, Ap=As=-106 dB, fp=0.238, fs = 0.25 (in terms of sample rate) + static const IIRCoefficient FIR_BEST = 0.0014313792470984f; + static const IIRSection SECTIONS_BEST[] = { + { 2.85800356692148000f,-0.2607342682253230f,-0.602478421807085f, 0.109823442522145f }, + { -4.39519408383016000f, 1.4651975326003500f,-0.533817668127954f, 0.226045921792036f }, + { 0.86638550740991800f,-2.1053851417898500f,-0.429134968401065f, 0.403512574222174f }, + { 1.67161485530774000f, 0.7963595880494520f,-0.324989203363446f, 0.580756666711889f }, + { -1.19962759276471000f, 0.5873595178851540f,-0.241486447489019f, 0.724264899930934f }, + { 0.01631779946479250f,-0.6282334739461620f,-0.182766025706656f, 0.827774001858882f }, + { 0.28404415859352400f, 0.1038619997715160f,-0.145276649558926f, 0.898510501923554f }, + { -0.08105788424234910f, 0.0781551578108934f,-0.123965846623366f, 0.947105257601873f }, + { -0.00872608905948005f,-0.0222098231712466f,-0.115056854360748f, 0.983542001125711f } + }; + + // Average elliptic filter with symmetric ripple: N=12, Ap=As=-106 dB, fp=0.193, fs = 0.25 (in terms of sample rate) + static const IIRCoefficient FIR_GOOD = 0.000891054570268146f; + static const IIRSection SECTIONS_GOOD[] = { + { 2.2650157226725700f,-0.4034180565140230f,-0.750061486095301f, 0.157801404511953f }, + { -3.2788261989161700f, 1.3952152147542600f,-0.705854270206788f, 0.265564985645774f }, + { 0.4397975114813240f,-1.3957634748753100f,-0.639718853965265f, 0.435324134360315f }, + { 0.9827040216680520f, 0.1837182774040940f,-0.578569965618418f, 0.615205557837542f }, + { -0.3759752818621670f, 0.3266073609399490f,-0.540913588637109f, 0.778264420176574f }, + { -0.0253548089519618f,-0.0925779221603846f,-0.537704370375240f, 0.925800083252964f } + }; + + // Fast elliptic filter with symmetric ripple: N=8, Ap=As=-99 dB, fp=0.125, fs = 0.25 (in terms of sample rate) + static const IIRCoefficient FIR_FAST = 0.000882837778745889f; + static const IIRSection SECTIONS_FAST[] = { + { 1.215377077431620f,-0.35864455030878000f,-0.972220718789242f, 0.252934735930620f }, + { -1.525654419254140f, 0.86784918631245500f,-0.977713689358124f, 0.376580616703668f }, + { 0.136094441564220f,-0.50414116798010400f,-1.007004471865290f, 0.584048854845331f }, + { 0.180604082285806f,-0.00467624342403851f,-1.093486919012100f, 0.844904524843996f } + }; + + static inline BufferedSample calcNumerator(const IIRSection §ion, const BufferedSample buffer1, const BufferedSample buffer2) { + return section.num1 * buffer1 + section.num2 * buffer2; + } + + static inline BufferedSample calcDenominator(const IIRSection §ion, const BufferedSample input, const BufferedSample buffer1, const BufferedSample buffer2) { + return input - section.den1 * buffer1 - section.den2 * buffer2; + } + +} // namespace SRCTools + +using namespace SRCTools; + +double IIRResampler::getPassbandFractionForQuality(Quality quality) { + switch (quality) { + case FAST: + return 0.5; + case GOOD: + return 0.7708; + case BEST: + return 0.9524; + default: + return 0; + } +} + +IIRResampler::Constants::Constants(const unsigned int useSectionsCount, const IIRCoefficient useFIR, const IIRSection useSections[], const Quality quality) { + if (quality == CUSTOM) { + sectionsCount = useSectionsCount; + fir = useFIR; + sections = useSections; + } else { + unsigned int sectionsSize; + switch (quality) { + case FAST: + fir = FIR_FAST; + sections = SECTIONS_FAST; + sectionsSize = sizeof(SECTIONS_FAST); + break; + case GOOD: + fir = FIR_GOOD; + sections = SECTIONS_GOOD; + sectionsSize = sizeof(SECTIONS_GOOD); + break; + case BEST: + fir = FIR_BEST; + sections = SECTIONS_BEST; + sectionsSize = sizeof(SECTIONS_BEST); + break; + default: + sectionsSize = 0; + break; + } + sectionsCount = (sectionsSize / sizeof(IIRSection)); + } + const unsigned int delayLineSize = IIR_RESAMPER_CHANNEL_COUNT * sectionsCount; + buffer = new SectionBuffer[delayLineSize]; + BufferedSample *s = buffer[0]; + BufferedSample *e = buffer[delayLineSize]; + while (s < e) *(s++) = 0; +} + +IIRResampler::IIRResampler(const Quality quality) : + constants(0, 0.0f, NULL, quality) +{} + +IIRResampler::IIRResampler(const unsigned int useSectionsCount, const IIRCoefficient useFIR, const IIRSection useSections[]) : + constants(useSectionsCount, useFIR, useSections, IIRResampler::CUSTOM) +{} + +IIRResampler::~IIRResampler() { + delete[] constants.buffer; +} + +IIR2xInterpolator::IIR2xInterpolator(const Quality quality) : + IIRResampler(quality), + phase() +{ + for (unsigned int chIx = 0; chIx < IIR_RESAMPER_CHANNEL_COUNT; ++chIx) { + lastInputSamples[chIx] = 0; + } +} + +IIR2xInterpolator::IIR2xInterpolator(const unsigned int useSectionsCount, const IIRCoefficient useFIR, const IIRSection useSections[]) : + IIRResampler(useSectionsCount, useFIR, useSections), + phase() +{ + for (unsigned int chIx = 0; chIx < IIR_RESAMPER_CHANNEL_COUNT; ++chIx) { + lastInputSamples[chIx] = 0; + } +} + +void IIR2xInterpolator::process(const FloatSample *&inSamples, unsigned int &inLength, FloatSample *&outSamples, unsigned int &outLength) { + static const IIRCoefficient INTERPOLATOR_AMP = 2.0; + + while (outLength > 0 && inLength > 0) { + SectionBuffer *bufferp = constants.buffer; + for (unsigned int chIx = 0; chIx < IIR_RESAMPER_CHANNEL_COUNT; ++chIx) { + const FloatSample lastInputSample = lastInputSamples[chIx]; + const FloatSample inSample = inSamples[chIx]; + BufferedSample tmpOut = phase == 0 ? 0 : inSample * constants.fir; + for (unsigned int i = 0; i < constants.sectionsCount; ++i) { + const IIRSection §ion = constants.sections[i]; + SectionBuffer &buffer = *bufferp; + // For 2x interpolation, calculation of the numerator reduces to a single multiplication depending on the phase. + if (phase == 0) { + const BufferedSample numOutSample = section.num1 * lastInputSample; + const BufferedSample denOutSample = calcDenominator(section, BIAS + numOutSample, buffer[0], buffer[1]); + buffer[1] = denOutSample; + tmpOut += denOutSample; + } else { + const BufferedSample numOutSample = section.num2 * lastInputSample; + const BufferedSample denOutSample = calcDenominator(section, BIAS + numOutSample, buffer[1], buffer[0]); + buffer[0] = denOutSample; + tmpOut += denOutSample; + } + bufferp++; + } + *(outSamples++) = FloatSample(INTERPOLATOR_AMP * tmpOut); + if (phase > 0) { + lastInputSamples[chIx] = inSample; + } + } + outLength--; + if (phase > 0) { + inSamples += IIR_RESAMPER_CHANNEL_COUNT; + inLength--; + phase = 0; + } else { + phase = 1; + } + } +} + +unsigned int IIR2xInterpolator::estimateInLength(const unsigned int outLength) const { + return outLength >> 1; +} + +IIR2xDecimator::IIR2xDecimator(const Quality quality) : + IIRResampler(quality) +{} + +IIR2xDecimator::IIR2xDecimator(const unsigned int useSectionsCount, const IIRCoefficient useFIR, const IIRSection useSections[]) : + IIRResampler(useSectionsCount, useFIR, useSections) +{} + +void IIR2xDecimator::process(const FloatSample *&inSamples, unsigned int &inLength, FloatSample *&outSamples, unsigned int &outLength) { + while (outLength > 0 && inLength > 1) { + SectionBuffer *bufferp = constants.buffer; + for (unsigned int chIx = 0; chIx < IIR_RESAMPER_CHANNEL_COUNT; ++chIx) { + BufferedSample tmpOut = inSamples[chIx] * constants.fir; + for (unsigned int i = 0; i < constants.sectionsCount; ++i) { + const IIRSection §ion = constants.sections[i]; + SectionBuffer &buffer = *bufferp; + // For 2x decimation, calculation of the numerator is not performed for odd output samples which are to be omitted. + tmpOut += calcNumerator(section, buffer[0], buffer[1]); + buffer[1] = calcDenominator(section, BIAS + inSamples[chIx], buffer[0], buffer[1]); + buffer[0] = calcDenominator(section, BIAS + inSamples[chIx + IIR_RESAMPER_CHANNEL_COUNT], buffer[1], buffer[0]); + bufferp++; + } + *(outSamples++) = FloatSample(tmpOut); + } + outLength--; + inLength -= 2; + inSamples += 2 * IIR_RESAMPER_CHANNEL_COUNT; + } +} + +unsigned int IIR2xDecimator::estimateInLength(const unsigned int outLength) const { + return outLength << 1; +} diff --git a/audio/softsynth/mt32/srchelper/srctools/src/LinearResampler.cpp b/audio/softsynth/mt32/srchelper/srctools/src/LinearResampler.cpp new file mode 100644 index 0000000000..98b9c77c73 --- /dev/null +++ b/audio/softsynth/mt32/srchelper/srctools/src/LinearResampler.cpp @@ -0,0 +1,47 @@ +/* Copyright (C) 2015-2017 Sergey V. Mikayev + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 2.1 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include "LinearResampler.h" + +using namespace SRCTools; + +LinearResampler::LinearResampler(double sourceSampleRate, double targetSampleRate) : + inputToOutputRatio(sourceSampleRate / targetSampleRate), + position(1.0) // Preload delay line which effectively makes resampler zero phase +{} + +void LinearResampler::process(const FloatSample *&inSamples, unsigned int &inLength, FloatSample *&outSamples, unsigned int &outLength) { + if (inLength == 0) return; + while (outLength > 0) { + while (1.0 <= position) { + position--; + inLength--; + for (unsigned int chIx = 0; chIx < LINEAR_RESAMPER_CHANNEL_COUNT; ++chIx) { + lastInputSamples[chIx] = *(inSamples++); + } + if (inLength == 0) return; + } + for (unsigned int chIx = 0; chIx < LINEAR_RESAMPER_CHANNEL_COUNT; chIx++) { + *(outSamples++) = FloatSample(lastInputSamples[chIx] + position * (inSamples[chIx] - lastInputSamples[chIx])); + } + outLength--; + position += inputToOutputRatio; + } +} + +unsigned int LinearResampler::estimateInLength(const unsigned int outLength) const { + return static_cast<unsigned int>(outLength * inputToOutputRatio); +} diff --git a/audio/softsynth/mt32/srchelper/srctools/src/ResamplerModel.cpp b/audio/softsynth/mt32/srchelper/srctools/src/ResamplerModel.cpp new file mode 100644 index 0000000000..4d2d930837 --- /dev/null +++ b/audio/softsynth/mt32/srchelper/srctools/src/ResamplerModel.cpp @@ -0,0 +1,153 @@ +/* Copyright (C) 2015-2017 Sergey V. Mikayev + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 2.1 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <cmath> +#include <cstddef> + +#include "ResamplerModel.h" + +#include "ResamplerStage.h" +#include "SincResampler.h" +#include "IIR2xResampler.h" +#include "LinearResampler.h" + +namespace SRCTools { + +namespace ResamplerModel { + +static const unsigned int CHANNEL_COUNT = 2; +static const unsigned int MAX_SAMPLES_PER_RUN = 4096; + +class CascadeStage : public FloatSampleProvider { +friend void freeResamplerModel(FloatSampleProvider &model, FloatSampleProvider &source); +public: + CascadeStage(FloatSampleProvider &source, ResamplerStage &resamplerStage); + + void getOutputSamples(FloatSample *outBuffer, unsigned int size); + +protected: + ResamplerStage &resamplerStage; + +private: + FloatSampleProvider &source; + FloatSample buffer[CHANNEL_COUNT * MAX_SAMPLES_PER_RUN]; + const FloatSample *bufferPtr; + unsigned int size; +}; + +class InternalResamplerCascadeStage : public CascadeStage { +public: + InternalResamplerCascadeStage(FloatSampleProvider &useSource, ResamplerStage &useResamplerStage) : + CascadeStage(useSource, useResamplerStage) + {} + + ~InternalResamplerCascadeStage() { + delete &resamplerStage; + } +}; + +} // namespace ResamplerModel + +} // namespace SRCTools + +using namespace SRCTools; + +FloatSampleProvider &ResamplerModel::createResamplerModel(FloatSampleProvider &source, double sourceSampleRate, double targetSampleRate, Quality quality) { + if (sourceSampleRate == targetSampleRate) { + return source; + } + if (quality == FASTEST) { + return *new InternalResamplerCascadeStage(source, *new LinearResampler(sourceSampleRate, targetSampleRate)); + } + const IIRResampler::Quality iirQuality = static_cast<IIRResampler::Quality>(quality); + const double iirPassbandFraction = IIRResampler::getPassbandFractionForQuality(iirQuality); + if (sourceSampleRate < targetSampleRate) { + ResamplerStage *iir2xInterpolator = new IIR2xInterpolator(iirQuality); + FloatSampleProvider &iir2xInterpolatorStage = *new InternalResamplerCascadeStage(source, *iir2xInterpolator); + + if (2.0 * sourceSampleRate == targetSampleRate) { + return iir2xInterpolatorStage; + } + + double passband = 0.5 * sourceSampleRate * iirPassbandFraction; + double stopband = 1.5 * sourceSampleRate; + ResamplerStage *sincResampler = SincResampler::createSincResampler(2.0 * sourceSampleRate, targetSampleRate, passband, stopband, DEFAULT_DB_SNR, DEFAULT_WINDOWED_SINC_MAX_UPSAMPLE_FACTOR); + return *new InternalResamplerCascadeStage(iir2xInterpolatorStage, *sincResampler); + } + + if (sourceSampleRate == 2.0 * targetSampleRate) { + ResamplerStage *iir2xDecimator = new IIR2xDecimator(iirQuality); + return *new InternalResamplerCascadeStage(source, *iir2xDecimator); + } + + double passband = 0.5 * targetSampleRate * iirPassbandFraction; + double stopband = 1.5 * targetSampleRate; + double sincOutSampleRate = 2.0 * targetSampleRate; + const unsigned int maxUpsampleFactor = static_cast<unsigned int>(ceil(DEFAULT_WINDOWED_SINC_MAX_DOWNSAMPLE_FACTOR * sincOutSampleRate / sourceSampleRate)); + ResamplerStage *sincResampler = SincResampler::createSincResampler(sourceSampleRate, sincOutSampleRate, passband, stopband, DEFAULT_DB_SNR, maxUpsampleFactor); + FloatSampleProvider &sincResamplerStage = *new InternalResamplerCascadeStage(source, *sincResampler); + + ResamplerStage *iir2xDecimator = new IIR2xDecimator(iirQuality); + return *new InternalResamplerCascadeStage(sincResamplerStage, *iir2xDecimator); +} + +FloatSampleProvider &ResamplerModel::createResamplerModel(FloatSampleProvider &source, ResamplerStage **resamplerStages, unsigned int stageCount) { + FloatSampleProvider *prevStage = &source; + for (unsigned int i = 0; i < stageCount; i++) { + prevStage = new CascadeStage(*prevStage, *(resamplerStages[i])); + } + return *prevStage; +} + +FloatSampleProvider &ResamplerModel::createResamplerModel(FloatSampleProvider &source, ResamplerStage &stage) { + return *new CascadeStage(source, stage); +} + +void ResamplerModel::freeResamplerModel(FloatSampleProvider &model, FloatSampleProvider &source) { + FloatSampleProvider *currentStage = &model; + while (currentStage != &source) { + CascadeStage *cascadeStage = dynamic_cast<CascadeStage *>(currentStage); + if (cascadeStage == NULL) return; + FloatSampleProvider &prevStage = cascadeStage->source; + delete currentStage; + currentStage = &prevStage; + } +} + +using namespace ResamplerModel; + +CascadeStage::CascadeStage(FloatSampleProvider &useSource, ResamplerStage &useResamplerStage) : + resamplerStage(useResamplerStage), + source(useSource), + bufferPtr(buffer), + size() +{} + +void CascadeStage::getOutputSamples(FloatSample *outBuffer, unsigned int length) { + while (length > 0) { + if (size == 0) { + size = resamplerStage.estimateInLength(length); + if (size < 1) { + size = 1; + } else if (MAX_SAMPLES_PER_RUN < size) { + size = MAX_SAMPLES_PER_RUN; + } + source.getOutputSamples(buffer, size); + bufferPtr = buffer; + } + resamplerStage.process(bufferPtr, size, outBuffer, length); + } +} diff --git a/audio/softsynth/mt32/srchelper/srctools/src/SincResampler.cpp b/audio/softsynth/mt32/srchelper/srctools/src/SincResampler.cpp new file mode 100644 index 0000000000..3ed028d261 --- /dev/null +++ b/audio/softsynth/mt32/srchelper/srctools/src/SincResampler.cpp @@ -0,0 +1,136 @@ +/* Copyright (C) 2015-2017 Sergey V. Mikayev + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation, either version 2.1 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + */ + +#include <cmath> + +#ifdef SINC_RESAMPLER_DEBUG_LOG +#include <iostream> +#endif + +#include "SincResampler.h" + +#ifndef M_PI +static const double M_PI = 3.1415926535897932; +#endif + +using namespace SRCTools; + +using namespace SincResampler; + +using namespace Utils; + +void Utils::computeResampleFactors(unsigned int &upsampleFactor, double &downsampleFactor, const double inputFrequency, const double outputFrequency, const unsigned int maxUpsampleFactor) { + static const double RATIONAL_RATIO_ACCURACY_FACTOR = 1E15; + + upsampleFactor = static_cast<unsigned int>(outputFrequency); + unsigned int downsampleFactorInt = static_cast<unsigned int>(inputFrequency); + if ((upsampleFactor == outputFrequency) && (downsampleFactorInt == inputFrequency)) { + // Input and output frequencies are integers, try to reduce them + const unsigned int gcd = greatestCommonDivisor(upsampleFactor, downsampleFactorInt); + if (gcd > 1) { + upsampleFactor /= gcd; + downsampleFactor = downsampleFactorInt / gcd; + } else { + downsampleFactor = downsampleFactorInt; + } + if (upsampleFactor <= maxUpsampleFactor) return; + } else { + // Try to recover rational resample ratio by brute force + double inputToOutputRatio = inputFrequency / outputFrequency; + for (unsigned int i = 1; i <= maxUpsampleFactor; ++i) { + double testFactor = inputToOutputRatio * i; + if (floor(RATIONAL_RATIO_ACCURACY_FACTOR * testFactor + 0.5) == RATIONAL_RATIO_ACCURACY_FACTOR * floor(testFactor + 0.5)) { + // inputToOutputRatio found to be rational within the accuracy + upsampleFactor = i; + downsampleFactor = floor(testFactor + 0.5); + return; + } + } + } + // Use interpolation of FIR taps as the last resort + upsampleFactor = maxUpsampleFactor; + downsampleFactor = maxUpsampleFactor * inputFrequency / outputFrequency; +} + +unsigned int Utils::greatestCommonDivisor(unsigned int a, unsigned int b) { + while (0 < b) { + unsigned int r = a % b; + a = b; + b = r; + } + return a; +} + +double KaizerWindow::estimateBeta(double dbRipple) { + return 0.1102 * (dbRipple - 8.7); +} + +unsigned int KaizerWindow::estimateOrder(double dbRipple, double fp, double fs) { + const double transBW = (fs - fp); + return static_cast<unsigned int>(ceil((dbRipple - 8) / (2.285 * 2 * M_PI * transBW))); +} + +double KaizerWindow::bessel(const double x) { + static const double EPS = 1.11E-16; + + double sum = 0.0; + double f = 1.0; + for (unsigned int i = 1;; ++i) { + f *= (0.5 * x / i); + double f2 = f * f; + if (f2 <= sum * EPS) break; + sum += f2; + } + return 1.0 + sum; +} + +void KaizerWindow::windowedSinc(FIRCoefficient kernel[], const unsigned int order, const double fc, const double beta, const double amp) { + const double fc_pi = M_PI * fc; + const double recipOrder = 1.0 / order; + const double mult = 2.0 * fc * amp / bessel(beta); + for (int i = order, j = 0; 0 <= i; i -= 2, ++j) { + double xw = i * recipOrder; + double win = bessel(beta * sqrt(fabs(1.0 - xw * xw))); + double xs = i * fc_pi; + double sinc = (i == 0) ? 1.0 : sin(xs) / xs; + FIRCoefficient imp = FIRCoefficient(mult * sinc * win); + kernel[j] = imp; + kernel[order - j] = imp; + } +} + +ResamplerStage *SincResampler::createSincResampler(const double inputFrequency, const double outputFrequency, const double passbandFrequency, const double stopbandFrequency, const double dbSNR, const unsigned int maxUpsampleFactor) { + unsigned int upsampleFactor; + double downsampleFactor; + computeResampleFactors(upsampleFactor, downsampleFactor, inputFrequency, outputFrequency, maxUpsampleFactor); + double baseSamplePeriod = 1.0 / (inputFrequency * upsampleFactor); + double fp = passbandFrequency * baseSamplePeriod; + double fs = stopbandFrequency * baseSamplePeriod; + double fc = 0.5 * (fp + fs); + double beta = KaizerWindow::estimateBeta(dbSNR); + unsigned int order = KaizerWindow::estimateOrder(dbSNR, fp, fs); + const unsigned int kernelLength = order + 1; + +#ifdef SINC_RESAMPLER_DEBUG_LOG + std::clog << "FIR: " << upsampleFactor << "/" << downsampleFactor << ", N=" << kernelLength << ", NPh=" << kernelLength / double(upsampleFactor) << ", C=" << 0.5 / fc << ", fp=" << fp << ", fs=" << fs << ", M=" << maxUpsampleFactor << std::endl; +#endif + + FIRCoefficient *windowedSincKernel = new FIRCoefficient[kernelLength]; + KaizerWindow::windowedSinc(windowedSincKernel, order, fc, beta, upsampleFactor); + ResamplerStage *windowedSincStage = new FIRResampler(upsampleFactor, downsampleFactor, windowedSincKernel, kernelLength); + delete[] windowedSincKernel; + return windowedSincStage; +} |