diff options
Diffstat (limited to 'libs/qm-dsp/dsp/rateconversion')
| -rw-r--r-- | libs/qm-dsp/dsp/rateconversion/Decimator.cpp | 22 | ||||
| -rw-r--r-- | libs/qm-dsp/dsp/rateconversion/Decimator.h | 34 | ||||
| -rw-r--r-- | libs/qm-dsp/dsp/rateconversion/DecimatorB.cpp | 160 | ||||
| -rw-r--r-- | libs/qm-dsp/dsp/rateconversion/DecimatorB.h | 64 | ||||
| -rw-r--r-- | libs/qm-dsp/dsp/rateconversion/Resampler.cpp | 416 | ||||
| -rw-r--r-- | libs/qm-dsp/dsp/rateconversion/Resampler.h | 102 |
6 files changed, 785 insertions, 13 deletions
diff --git a/libs/qm-dsp/dsp/rateconversion/Decimator.cpp b/libs/qm-dsp/dsp/rateconversion/Decimator.cpp index c150ee0e11..593474f7be 100644 --- a/libs/qm-dsp/dsp/rateconversion/Decimator.cpp +++ b/libs/qm-dsp/dsp/rateconversion/Decimator.cpp @@ -199,10 +199,15 @@ void Decimator::doAntiAlias(const float *src, double *dst, unsigned int length) void Decimator::process(const double *src, double *dst) { - if( m_decFactor != 1 ) - { - doAntiAlias( src, decBuffer, m_inputLength ); + if (m_decFactor == 1) { + for( unsigned int i = 0; i < m_outputLength; i++ ) { + dst[i] = src[i]; + } + return; } + + doAntiAlias( src, decBuffer, m_inputLength ); + unsigned idx = 0; for( unsigned int i = 0; i < m_outputLength; i++ ) @@ -213,10 +218,15 @@ void Decimator::process(const double *src, double *dst) void Decimator::process(const float *src, float *dst) { - if( m_decFactor != 1 ) - { - doAntiAlias( src, decBuffer, m_inputLength ); + if (m_decFactor == 1) { + for( unsigned int i = 0; i < m_outputLength; i++ ) { + dst[i] = src[i]; + } + return; } + + doAntiAlias( src, decBuffer, m_inputLength ); + unsigned idx = 0; for( unsigned int i = 0; i < m_outputLength; i++ ) diff --git a/libs/qm-dsp/dsp/rateconversion/Decimator.h b/libs/qm-dsp/dsp/rateconversion/Decimator.h index f8a113a0db..e3913f8722 100644 --- a/libs/qm-dsp/dsp/rateconversion/Decimator.h +++ b/libs/qm-dsp/dsp/rateconversion/Decimator.h @@ -15,12 +15,15 @@ #ifndef DECIMATOR_H #define DECIMATOR_H -class Decimator +/** + * Decimator carries out a fast downsample by a power-of-two + * factor. Only a limited number of factors are supported, from two to + * whatever getHighestSupportedFactor() returns. This is much faster + * than Resampler but has a worse signal-noise ratio. + */ +class Decimator { public: - void process( const double* src, double* dst ); - void process( const float* src, float* dst ); - /** * Construct a Decimator to operate on input blocks of length * inLength, with decimation factor decFactor. inLength should be @@ -34,11 +37,28 @@ public: Decimator( unsigned int inLength, unsigned int decFactor ); virtual ~Decimator(); + /** + * Process inLength samples (as supplied to constructor) from src + * and write inLength / decFactor samples to dst. Note that src + * and dst may be the same or overlap (an intermediate buffer is + * used). + */ + void process( const double* src, double* dst ); + + /** + * Process inLength samples (as supplied to constructor) from src + * and write inLength / decFactor samples to dst. Note that src + * and dst may be the same or overlap (an intermediate buffer is + * used). + */ + void process( const float* src, float* dst ); + int getFactor() const { return m_decFactor; } static int getHighestSupportedFactor() { return 8; } -private: void resetFilter(); + +private: void deInitialise(); void initialise( unsigned int inLength, unsigned int decFactor ); void doAntiAlias( const double* src, double* dst, unsigned int length ); @@ -55,8 +75,8 @@ private: double a[ 9 ]; double b[ 9 ]; - + double* decBuffer; }; -#endif // +#endif // diff --git a/libs/qm-dsp/dsp/rateconversion/DecimatorB.cpp b/libs/qm-dsp/dsp/rateconversion/DecimatorB.cpp new file mode 100644 index 0000000000..55df1e9fc0 --- /dev/null +++ b/libs/qm-dsp/dsp/rateconversion/DecimatorB.cpp @@ -0,0 +1,160 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ + +/* + QM DSP Library + + Centre for Digital Music, Queen Mary, University of London. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. See the file + COPYING included with this distribution for more information. +*/ + +#include "DecimatorB.h" + +#include "maths/MathUtilities.h" + +#include <iostream> + +using std::vector; + +DecimatorB::DecimatorB(int inLength, int decFactor) +{ + m_inputLength = 0; + m_outputLength = 0; + m_decFactor = 1; + m_aaBuffer = 0; + m_tmpBuffer = 0; + + initialise(inLength, decFactor); +} + +DecimatorB::~DecimatorB() +{ + deInitialise(); +} + +void DecimatorB::initialise(int inLength, int decFactor) +{ + m_inputLength = inLength; + m_decFactor = decFactor; + m_outputLength = m_inputLength / m_decFactor; + + if (m_decFactor < 2 || !MathUtilities::isPowerOfTwo(m_decFactor)) { + std::cerr << "ERROR: DecimatorB::initialise: Decimation factor must be a power of 2 and at least 2 (was: " << m_decFactor << ")" << std::endl; + m_decFactor = 0; + return; + } + + if (m_inputLength % m_decFactor != 0) { + std::cerr << "ERROR: DecimatorB::initialise: inLength must be a multiple of decimation factor (was: " << m_inputLength << ", factor is " << m_decFactor << ")" << std::endl; + m_decFactor = 0; + return; + } + + m_aaBuffer = new double[m_inputLength]; + m_tmpBuffer = new double[m_inputLength]; + + // Order 6 Butterworth lowpass filter + // Calculated using e.g. MATLAB butter(6, 0.5, 'low') + + m_b[0] = 0.029588223638661; + m_b[1] = 0.177529341831965; + m_b[2] = 0.443823354579912; + m_b[3] = 0.591764472773216; + m_b[4] = 0.443823354579912; + m_b[5] = 0.177529341831965; + m_b[6] = 0.029588223638661; + + m_a[0] = 1.000000000000000; + m_a[1] = 0.000000000000000; + m_a[2] = 0.777695961855673; + m_a[3] = 0.000000000000000; + m_a[4] = 0.114199425062434; + m_a[5] = 0.000000000000000; + m_a[6] = 0.001750925956183; + + for (int factor = m_decFactor; factor > 1; factor /= 2) { + m_o.push_back(vector<double>(6, 0.0)); + } +} + +void DecimatorB::deInitialise() +{ + delete [] m_aaBuffer; + delete [] m_tmpBuffer; +} + +void DecimatorB::doAntiAlias(const double *src, double *dst, int length, + int filteridx) +{ + vector<double> &o = m_o[filteridx]; + + for (int i = 0; i < length; i++) { + + double input = src[i]; + double output = input * m_b[0] + o[0]; + + o[0] = input * m_b[1] - output * m_a[1] + o[1]; + o[1] = input * m_b[2] - output * m_a[2] + o[2]; + o[2] = input * m_b[3] - output * m_a[3] + o[3]; + o[3] = input * m_b[4] - output * m_a[4] + o[4]; + o[4] = input * m_b[5] - output * m_a[5] + o[5]; + o[5] = input * m_b[6] - output * m_a[6]; + + dst[i] = output; + } +} + +void DecimatorB::doProcess() +{ + int filteridx = 0; + int factorDone = 1; + int factorRemaining = m_decFactor; + + while (factorDone < m_decFactor) { + + doAntiAlias(m_tmpBuffer, m_aaBuffer, + m_inputLength / factorDone, + filteridx); + + filteridx ++; + factorDone *= 2; + + for (int i = 0; i < m_inputLength / factorDone; ++i) { + m_tmpBuffer[i] = m_aaBuffer[i * 2]; + } + } +} + +void DecimatorB::process(const double *src, double *dst) +{ + if (m_decFactor == 0) return; + + for (int i = 0; i < m_inputLength; ++i) { + m_tmpBuffer[i] = src[i]; + } + + doProcess(); + + for (int i = 0; i < m_outputLength; ++i) { + dst[i] = m_tmpBuffer[i]; + } +} + +void DecimatorB::process(const float *src, float *dst) +{ + if (m_decFactor == 0) return; + + for (int i = 0; i < m_inputLength; ++i) { + m_tmpBuffer[i] = src[i]; + } + + doProcess(); + + for (int i = 0; i < m_outputLength; ++i) { + dst[i] = m_tmpBuffer[i]; + } +} diff --git a/libs/qm-dsp/dsp/rateconversion/DecimatorB.h b/libs/qm-dsp/dsp/rateconversion/DecimatorB.h new file mode 100644 index 0000000000..8458e61061 --- /dev/null +++ b/libs/qm-dsp/dsp/rateconversion/DecimatorB.h @@ -0,0 +1,64 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + QM DSP Library + + Centre for Digital Music, Queen Mary, University of London. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. See the file + COPYING included with this distribution for more information. +*/ + +#ifndef DECIMATORB_H +#define DECIMATORB_H + +#include <vector> + +/** + * DecimatorB carries out a fast downsample by a power-of-two + * factor. It only knows how to decimate by a factor of 2, and will + * use repeated decimation for higher factors. A Butterworth filter of + * order 6 is used for the lowpass filter. + */ +class DecimatorB +{ +public: + void process( const double* src, double* dst ); + void process( const float* src, float* dst ); + + /** + * Construct a DecimatorB to operate on input blocks of length + * inLength, with decimation factor decFactor. inLength should be + * a multiple of decFactor. Output blocks will be of length + * inLength / decFactor. + * + * decFactor must be a power of two. + */ + DecimatorB(int inLength, int decFactor); + virtual ~DecimatorB(); + + int getFactor() const { return m_decFactor; } + +private: + void deInitialise(); + void initialise(int inLength, int decFactor); + void doAntiAlias(const double* src, double* dst, int length, int filteridx); + void doProcess(); + + int m_inputLength; + int m_outputLength; + int m_decFactor; + + std::vector<std::vector<double> > m_o; + + double m_a[7]; + double m_b[7]; + + double *m_aaBuffer; + double *m_tmpBuffer; +}; + +#endif + diff --git a/libs/qm-dsp/dsp/rateconversion/Resampler.cpp b/libs/qm-dsp/dsp/rateconversion/Resampler.cpp new file mode 100644 index 0000000000..f0598cab2c --- /dev/null +++ b/libs/qm-dsp/dsp/rateconversion/Resampler.cpp @@ -0,0 +1,416 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + QM DSP Library + + Centre for Digital Music, Queen Mary, University of London. + This file by Chris Cannam. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. See the file + COPYING included with this distribution for more information. +*/ + +#include "Resampler.h" + +#include "maths/MathUtilities.h" +#include "base/KaiserWindow.h" +#include "base/SincWindow.h" +#include "thread/Thread.h" + +#include <iostream> +#include <vector> +#include <map> +#include <cassert> + +using std::vector; +using std::map; +using std::cerr; +using std::endl; + +//#define DEBUG_RESAMPLER 1 +//#define DEBUG_RESAMPLER_VERBOSE 1 + +Resampler::Resampler(int sourceRate, int targetRate) : + m_sourceRate(sourceRate), + m_targetRate(targetRate) +{ + initialise(100, 0.02); +} + +Resampler::Resampler(int sourceRate, int targetRate, + double snr, double bandwidth) : + m_sourceRate(sourceRate), + m_targetRate(targetRate) +{ + initialise(snr, bandwidth); +} + +Resampler::~Resampler() +{ + delete[] m_phaseData; +} + +// peakToPole -> length -> beta -> window +static map<double, map<int, map<double, vector<double> > > > +knownFilters; + +static Mutex +knownFilterMutex; + +void +Resampler::initialise(double snr, double bandwidth) +{ + int higher = std::max(m_sourceRate, m_targetRate); + int lower = std::min(m_sourceRate, m_targetRate); + + m_gcd = MathUtilities::gcd(lower, higher); + m_peakToPole = higher / m_gcd; + + if (m_targetRate < m_sourceRate) { + // antialiasing filter, should be slightly below nyquist + m_peakToPole = m_peakToPole / (1.0 - bandwidth/2.0); + } + + KaiserWindow::Parameters params = + KaiserWindow::parametersForBandwidth(snr, bandwidth, higher / m_gcd); + + params.length = + (params.length % 2 == 0 ? params.length + 1 : params.length); + + params.length = + (params.length > 200001 ? 200001 : params.length); + + m_filterLength = params.length; + + vector<double> filter; + knownFilterMutex.lock(); + + if (knownFilters[m_peakToPole][m_filterLength].find(params.beta) == + knownFilters[m_peakToPole][m_filterLength].end()) { + + KaiserWindow kw(params); + SincWindow sw(m_filterLength, m_peakToPole * 2); + + filter = vector<double>(m_filterLength, 0.0); + for (int i = 0; i < m_filterLength; ++i) filter[i] = 1.0; + sw.cut(filter.data()); + kw.cut(filter.data()); + + knownFilters[m_peakToPole][m_filterLength][params.beta] = filter; + } + + filter = knownFilters[m_peakToPole][m_filterLength][params.beta]; + knownFilterMutex.unlock(); + + int inputSpacing = m_targetRate / m_gcd; + int outputSpacing = m_sourceRate / m_gcd; + +#ifdef DEBUG_RESAMPLER + cerr << "resample " << m_sourceRate << " -> " << m_targetRate + << ": inputSpacing " << inputSpacing << ", outputSpacing " + << outputSpacing << ": filter length " << m_filterLength + << endl; +#endif + + // Now we have a filter of (odd) length flen in which the lower + // sample rate corresponds to every n'th point and the higher rate + // to every m'th where n and m are higher and lower rates divided + // by their gcd respectively. So if x coordinates are on the same + // scale as our filter resolution, then source sample i is at i * + // (targetRate / gcd) and target sample j is at j * (sourceRate / + // gcd). + + // To reconstruct a single target sample, we want a buffer (real + // or virtual) of flen values formed of source samples spaced at + // intervals of (targetRate / gcd), in our example case 3. This + // is initially formed with the first sample at the filter peak. + // + // 0 0 0 0 a 0 0 b 0 + // + // and of course we have our filter + // + // f1 f2 f3 f4 f5 f6 f7 f8 f9 + // + // We take the sum of products of non-zero values from this buffer + // with corresponding values in the filter + // + // a * f5 + b * f8 + // + // Then we drop (sourceRate / gcd) values, in our example case 4, + // from the start of the buffer and fill until it has flen values + // again + // + // a 0 0 b 0 0 c 0 0 + // + // repeat to reconstruct the next target sample + // + // a * f1 + b * f4 + c * f7 + // + // and so on. + // + // Above I said the buffer could be "real or virtual" -- ours is + // virtual. We don't actually store all the zero spacing values, + // except for padding at the start; normally we store only the + // values that actually came from the source stream, along with a + // phase value that tells us how many virtual zeroes there are at + // the start of the virtual buffer. So the two examples above are + // + // 0 a b [ with phase 1 ] + // a b c [ with phase 0 ] + // + // Having thus broken down the buffer so that only the elements we + // need to multiply are present, we can also unzip the filter into + // every-nth-element subsets at each phase, allowing us to do the + // filter multiplication as a simply vector multiply. That is, rather + // than store + // + // f1 f2 f3 f4 f5 f6 f7 f8 f9 + // + // we store separately + // + // f1 f4 f7 + // f2 f5 f8 + // f3 f6 f9 + // + // Each time we complete a multiply-and-sum, we need to work out + // how many (real) samples to drop from the start of our buffer, + // and how many to add at the end of it for the next multiply. We + // know we want to drop enough real samples to move along by one + // computed output sample, which is our outputSpacing number of + // virtual buffer samples. Depending on the relationship between + // input and output spacings, this may mean dropping several real + // samples, one real sample, or none at all (and simply moving to + // a different "phase"). + + m_phaseData = new Phase[inputSpacing]; + + for (int phase = 0; phase < inputSpacing; ++phase) { + + Phase p; + + p.nextPhase = phase - outputSpacing; + while (p.nextPhase < 0) p.nextPhase += inputSpacing; + p.nextPhase %= inputSpacing; + + p.drop = int(ceil(std::max(0.0, double(outputSpacing - phase)) + / inputSpacing)); + + int filtZipLength = int(ceil(double(m_filterLength - phase) + / inputSpacing)); + + for (int i = 0; i < filtZipLength; ++i) { + p.filter.push_back(filter[i * inputSpacing + phase]); + } + + m_phaseData[phase] = p; + } + +#ifdef DEBUG_RESAMPLER + int cp = 0; + int totDrop = 0; + for (int i = 0; i < inputSpacing; ++i) { + cerr << "phase = " << cp << ", drop = " << m_phaseData[cp].drop + << ", filter length = " << m_phaseData[cp].filter.size() + << ", next phase = " << m_phaseData[cp].nextPhase << endl; + totDrop += m_phaseData[cp].drop; + cp = m_phaseData[cp].nextPhase; + } + cerr << "total drop = " << totDrop << endl; +#endif + + // The May implementation of this uses a pull model -- we ask the + // resampler for a certain number of output samples, and it asks + // its source stream for as many as it needs to calculate + // those. This means (among other things) that the source stream + // can be asked for enough samples up-front to fill the buffer + // before the first output sample is generated. + // + // In this implementation we're using a push model in which a + // certain number of source samples is provided and we're asked + // for as many output samples as that makes available. But we + // can't return any samples from the beginning until half the + // filter length has been provided as input. This means we must + // either return a very variable number of samples (none at all + // until the filter fills, then half the filter length at once) or + // else have a lengthy declared latency on the output. We do the + // latter. (What do other implementations do?) + // + // We want to make sure the first "real" sample will eventually be + // aligned with the centre sample in the filter (it's tidier, and + // easier to do diagnostic calculations that way). So we need to + // pick the initial phase and buffer fill accordingly. + // + // Example: if the inputSpacing is 2, outputSpacing is 3, and + // filter length is 7, + // + // x x x x a b c ... input samples + // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 ... + // i j k l ... output samples + // [--------|--------] <- filter with centre mark + // + // Let h be the index of the centre mark, here 3 (generally + // int(filterLength/2) for odd-length filters). + // + // The smallest n such that h + n * outputSpacing > filterLength + // is 2 (that is, ceil((filterLength - h) / outputSpacing)), and + // (h + 2 * outputSpacing) % inputSpacing == 1, so the initial + // phase is 1. + // + // To achieve our n, we need to pre-fill the "virtual" buffer with + // 4 zero samples: the x's above. This is int((h + n * + // outputSpacing) / inputSpacing). It's the phase that makes this + // buffer get dealt with in such a way as to give us an effective + // index for sample a of 9 rather than 8 or 10 or whatever. + // + // This gives us output latency of 2 (== n), i.e. output samples i + // and j will appear before the one in which input sample a is at + // the centre of the filter. + + int h = int(m_filterLength / 2); + int n = ceil(double(m_filterLength - h) / outputSpacing); + + m_phase = (h + n * outputSpacing) % inputSpacing; + + int fill = (h + n * outputSpacing) / inputSpacing; + + m_latency = n; + + m_buffer = vector<double>(fill, 0); + m_bufferOrigin = 0; + +#ifdef DEBUG_RESAMPLER + cerr << "initial phase " << m_phase << " (as " << (m_filterLength/2) << " % " << inputSpacing << ")" + << ", latency " << m_latency << endl; +#endif +} + +double +Resampler::reconstructOne() +{ + Phase &pd = m_phaseData[m_phase]; + double v = 0.0; + int n = pd.filter.size(); + + assert(n + m_bufferOrigin <= (int)m_buffer.size()); + + const double *const __restrict__ buf = m_buffer.data() + m_bufferOrigin; + const double *const __restrict__ filt = pd.filter.data(); + + for (int i = 0; i < n; ++i) { + // NB gcc can only vectorize this with -ffast-math + v += buf[i] * filt[i]; + } + + m_bufferOrigin += pd.drop; + m_phase = pd.nextPhase; + return v; +} + +int +Resampler::process(const double *src, double *dst, int n) +{ + for (int i = 0; i < n; ++i) { + m_buffer.push_back(src[i]); + } + + int maxout = int(ceil(double(n) * m_targetRate / m_sourceRate)); + int outidx = 0; + +#ifdef DEBUG_RESAMPLER + cerr << "process: buf siz " << m_buffer.size() << " filt siz for phase " << m_phase << " " << m_phaseData[m_phase].filter.size() << endl; +#endif + + double scaleFactor = (double(m_targetRate) / m_gcd) / m_peakToPole; + + while (outidx < maxout && + m_buffer.size() >= m_phaseData[m_phase].filter.size() + m_bufferOrigin) { + dst[outidx] = scaleFactor * reconstructOne(); + outidx++; + } + + m_buffer = vector<double>(m_buffer.begin() + m_bufferOrigin, m_buffer.end()); + m_bufferOrigin = 0; + + return outidx; +} + +vector<double> +Resampler::process(const double *src, int n) +{ + int maxout = int(ceil(double(n) * m_targetRate / m_sourceRate)); + vector<double> out(maxout, 0.0); + int got = process(src, out.data(), n); + assert(got <= maxout); + if (got < maxout) out.resize(got); + return out; +} + +vector<double> +Resampler::resample(int sourceRate, int targetRate, const double *data, int n) +{ + Resampler r(sourceRate, targetRate); + + int latency = r.getLatency(); + + // latency is the output latency. We need to provide enough + // padding input samples at the end of input to guarantee at + // *least* the latency's worth of output samples. that is, + + int inputPad = int(ceil((double(latency) * sourceRate) / targetRate)); + + // that means we are providing this much input in total: + + int n1 = n + inputPad; + + // and obtaining this much output in total: + + int m1 = int(ceil((double(n1) * targetRate) / sourceRate)); + + // in order to return this much output to the user: + + int m = int(ceil((double(n) * targetRate) / sourceRate)); + +#ifdef DEBUG_RESAMPLER + cerr << "n = " << n << ", sourceRate = " << sourceRate << ", targetRate = " << targetRate << ", m = " << m << ", latency = " << latency << ", inputPad = " << inputPad << ", m1 = " << m1 << ", n1 = " << n1 << ", n1 - n = " << n1 - n << endl; +#endif + + vector<double> pad(n1 - n, 0.0); + vector<double> out(m1 + 1, 0.0); + + int gotData = r.process(data, out.data(), n); + int gotPad = r.process(pad.data(), out.data() + gotData, pad.size()); + int got = gotData + gotPad; + +#ifdef DEBUG_RESAMPLER + cerr << "resample: " << n << " in, " << pad.size() << " padding, " << got << " out (" << gotData << " data, " << gotPad << " padding, latency = " << latency << ")" << endl; +#endif +#ifdef DEBUG_RESAMPLER_VERBOSE + int printN = 50; + cerr << "first " << printN << " in:" << endl; + for (int i = 0; i < printN && i < n; ++i) { + if (i % 5 == 0) cerr << endl << i << "... "; + cerr << data[i] << " "; + } + cerr << endl; +#endif + + int toReturn = got - latency; + if (toReturn > m) toReturn = m; + + vector<double> sliced(out.begin() + latency, + out.begin() + latency + toReturn); + +#ifdef DEBUG_RESAMPLER_VERBOSE + cerr << "first " << printN << " out (after latency compensation), length " << sliced.size() << ":"; + for (int i = 0; i < printN && i < sliced.size(); ++i) { + if (i % 5 == 0) cerr << endl << i << "... "; + cerr << sliced[i] << " "; + } + cerr << endl; +#endif + + return sliced; +} + diff --git a/libs/qm-dsp/dsp/rateconversion/Resampler.h b/libs/qm-dsp/dsp/rateconversion/Resampler.h new file mode 100644 index 0000000000..92c0169ba0 --- /dev/null +++ b/libs/qm-dsp/dsp/rateconversion/Resampler.h @@ -0,0 +1,102 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +/* + QM DSP Library + + Centre for Digital Music, Queen Mary, University of London. + This file by Chris Cannam. + + This program is free software; you can redistribute it and/or + modify it under the terms of the GNU General Public License as + published by the Free Software Foundation; either version 2 of the + License, or (at your option) any later version. See the file + COPYING included with this distribution for more information. +*/ + +#ifndef RESAMPLER_H +#define RESAMPLER_H + +#include <vector> + +/** + * Resampler resamples a stream from one integer sample rate to + * another (arbitrary) rate, using a kaiser-windowed sinc filter. The + * results and performance are pretty similar to libraries such as + * libsamplerate, though this implementation does not support + * time-varying ratios (the ratio is fixed on construction). + * + * See also Decimator, which is faster and rougher but supports only + * power-of-two downsampling factors. + */ +class Resampler +{ +public: + /** + * Construct a Resampler to resample from sourceRate to + * targetRate. + */ + Resampler(int sourceRate, int targetRate); + + /** + * Construct a Resampler to resample from sourceRate to + * targetRate, using the given filter parameters. + */ + Resampler(int sourceRate, int targetRate, + double snr, double bandwidth); + + virtual ~Resampler(); + + /** + * Read n input samples from src and write resampled data to + * dst. The return value is the number of samples written, which + * will be no more than ceil((n * targetRate) / sourceRate). The + * caller must ensure the dst buffer has enough space for the + * samples returned. + */ + int process(const double *src, double *dst, int n); + + /** + * Read n input samples from src and return resampled data by + * value. + */ + std::vector<double> process(const double *src, int n); + + /** + * Return the number of samples of latency at the output due by + * the filter. (That is, the output will be delayed by this number + * of samples relative to the input.) + */ + int getLatency() const { return m_latency; } + + /** + * Carry out a one-off resample of a single block of n + * samples. The output is latency-compensated. + */ + static std::vector<double> resample + (int sourceRate, int targetRate, const double *data, int n); + +private: + int m_sourceRate; + int m_targetRate; + int m_gcd; + int m_filterLength; + int m_bufferLength; + int m_latency; + double m_peakToPole; + + struct Phase { + int nextPhase; + std::vector<double> filter; + int drop; + }; + + Phase *m_phaseData; + int m_phase; + std::vector<double> m_buffer; + int m_bufferOrigin; + + void initialise(double, double); + double reconstructOne(); +}; + +#endif + |