diff --git a/libs/qm-dsp/README.txt b/libs/qm-dsp/README.txt index 557ac23294..5856500f49 100644 --- a/libs/qm-dsp/README.txt +++ b/libs/qm-dsp/README.txt @@ -1,7 +1,6 @@ This is a stripped down version of the Queen Mary DSP library -Version 1.7.1 from Sept 2015. -https://code.soundsoftware.ac.uk/attachments/download/1582/qm-dsp-1.7.1.tar.gz +https://github.com/c4dm/qm-dsp -- see gitrev.txt for version --- diff --git a/libs/qm-dsp/base/KaiserWindow.h b/libs/qm-dsp/base/KaiserWindow.h index 0253d6d4c4..f16a4b6c16 100644 --- a/libs/qm-dsp/base/KaiserWindow.h +++ b/libs/qm-dsp/base/KaiserWindow.h @@ -81,7 +81,7 @@ public: } const double *getWindow() const { - return &m_window[0]; + return m_window.data(); } void cut(double *src) const { diff --git a/libs/qm-dsp/base/SincWindow.h b/libs/qm-dsp/base/SincWindow.h index 02de928679..bb35d90c20 100644 --- a/libs/qm-dsp/base/SincWindow.h +++ b/libs/qm-dsp/base/SincWindow.h @@ -37,7 +37,7 @@ public: } const double *getWindow() const { - return &m_window[0]; + return m_window.data(); } void cut(double *src) const { diff --git a/libs/qm-dsp/dsp/chromagram/Chromagram.cpp b/libs/qm-dsp/dsp/chromagram/Chromagram.cpp index a8597a5ddd..f42a4d2aad 100644 --- a/libs/qm-dsp/dsp/chromagram/Chromagram.cpp +++ b/libs/qm-dsp/dsp/chromagram/Chromagram.cpp @@ -34,7 +34,7 @@ int Chromagram::initialise( ChromaConfig Config ) m_normalise = Config.normalise; // if frame normalisation is required // No. of constant Q bins - m_uK = ( unsigned int ) ceil( m_BPO * log(m_FMax/m_FMin)/log(2.0)); + m_uK = (int) ceil( m_BPO * log(m_FMax/m_FMin)/log(2.0)); // Create array for chroma result m_chromadata = new double[ m_BPO ]; @@ -112,7 +112,7 @@ void Chromagram::unityNormalise(double *src) MathUtilities::getFrameMinMax( src, m_BPO, & min, &max ); - for( unsigned int i = 0; i < m_BPO; i++ ) + for (int i = 0; i < m_BPO; i++) { val = src[ i ] / max; @@ -153,19 +153,17 @@ double* Chromagram::process( const double *real, const double *imag ) } // initialise chromadata to 0 - for (unsigned i = 0; i < m_BPO; i++) m_chromadata[i] = 0; + for (int i = 0; i < m_BPO; i++) m_chromadata[i] = 0; - double cmax = 0.0; - double cval = 0; // Calculate ConstantQ frame m_ConstantQ->process( real, imag, m_CQRe, m_CQIm ); // add each octave of cq data into Chromagram - const unsigned octaves = (int)floor(double( m_uK/m_BPO))-1; - for (unsigned octave = 0; octave <= octaves; octave++) + const int octaves = (int)floor(double( m_uK/m_BPO))-1; + for (int octave = 0; octave <= octaves; octave++) { - unsigned firstBin = octave*m_BPO; - for (unsigned i = 0; i < m_BPO; i++) + int firstBin = octave*m_BPO; + for (int i = 0; i < m_BPO; i++) { m_chromadata[i] += kabs( m_CQRe[ firstBin + i ], m_CQIm[ firstBin + i ]); } diff --git a/libs/qm-dsp/dsp/chromagram/Chromagram.h b/libs/qm-dsp/dsp/chromagram/Chromagram.h index bd928f5d48..790abd8e23 100644 --- a/libs/qm-dsp/dsp/chromagram/Chromagram.h +++ b/libs/qm-dsp/dsp/chromagram/Chromagram.h @@ -21,10 +21,10 @@ #include "ConstantQ.h" struct ChromaConfig{ - unsigned int FS; + int FS; double min; double max; - unsigned int BPO; + int BPO; double CQThresh; MathUtilities::NormaliseType normalise; }; @@ -44,10 +44,10 @@ public: double kabs( double real, double imag ); // Results - unsigned int getK() { return m_uK;} - unsigned int getFrameSize() { return m_frameSize; } - unsigned int getHopSize() { return m_hopSize; } - + int getK() { return m_uK;} + int getFrameSize() { return m_frameSize; } + int getHopSize() { return m_hopSize; } + private: int initialise( ChromaConfig Config ); int deInitialise(); @@ -58,13 +58,13 @@ private: double* m_chromadata; double m_FMin; double m_FMax; - unsigned int m_BPO; - unsigned int m_uK; + int m_BPO; + int m_uK; MathUtilities::NormaliseType m_normalise; - unsigned int m_frameSize; - unsigned int m_hopSize; + int m_frameSize; + int m_hopSize; FFTReal* m_FFT; ConstantQ* m_ConstantQ; diff --git a/libs/qm-dsp/dsp/rateconversion/DecimatorB.cpp b/libs/qm-dsp/dsp/rateconversion/DecimatorB.cpp index 55df1e9fc0..5f27130c49 100644 --- a/libs/qm-dsp/dsp/rateconversion/DecimatorB.cpp +++ b/libs/qm-dsp/dsp/rateconversion/DecimatorB.cpp @@ -112,7 +112,6 @@ void DecimatorB::doProcess() { int filteridx = 0; int factorDone = 1; - int factorRemaining = m_decFactor; while (factorDone < m_decFactor) { diff --git a/libs/qm-dsp/dsp/signalconditioning/DFProcess.cpp b/libs/qm-dsp/dsp/signalconditioning/DFProcess.cpp index 52443fb384..7aa9649e6c 100644 --- a/libs/qm-dsp/dsp/signalconditioning/DFProcess.cpp +++ b/libs/qm-dsp/dsp/signalconditioning/DFProcess.cpp @@ -59,13 +59,11 @@ void DFProcess::initialise( DFProcConfig Config ) filtSrc = new double[ m_length ]; filtDst = new double[ m_length ]; - - //Low Pass Smoothing Filter Config - m_FilterConfigParams.ord = Config.LPOrd; - m_FilterConfigParams.ACoeffs = Config.LPACoeffs; - m_FilterConfigParams.BCoeffs = Config.LPBCoeffs; - - m_FiltFilt = new FiltFilt( m_FilterConfigParams ); + Filter::Parameters params; + params.a = std::vector(Config.LPACoeffs, Config.LPACoeffs + Config.LPOrd + 1); + params.b = std::vector(Config.LPBCoeffs, Config.LPBCoeffs + Config.LPOrd + 1); + + m_FiltFilt = new FiltFilt(params); //add delta threshold m_delta = Config.delta; @@ -193,7 +191,7 @@ void DFProcess::removeDCNormalize( double *src, double*dst ) MathUtilities::getAlphaNorm( src, m_length, m_alphaNormParam, &DFAlphaNorm ); - for( unsigned int i = 0; i< m_length; i++) + for (int i = 0; i < m_length; i++) { dst[ i ] = ( src[ i ] - DFMin ) / DFAlphaNorm; } diff --git a/libs/qm-dsp/dsp/signalconditioning/DFProcess.h b/libs/qm-dsp/dsp/signalconditioning/DFProcess.h index b93b535518..b256eca16e 100644 --- a/libs/qm-dsp/dsp/signalconditioning/DFProcess.h +++ b/libs/qm-dsp/dsp/signalconditioning/DFProcess.h @@ -81,8 +81,6 @@ private: double* m_filtScratchIn; double* m_filtScratchOut; - FilterConfig m_FilterConfigParams; - FiltFilt* m_FiltFilt; bool m_isMedianPositive; diff --git a/libs/qm-dsp/dsp/signalconditioning/FiltFilt.cpp b/libs/qm-dsp/dsp/signalconditioning/FiltFilt.cpp index 03b21138a4..c88641de7b 100644 --- a/libs/qm-dsp/dsp/signalconditioning/FiltFilt.cpp +++ b/libs/qm-dsp/dsp/signalconditioning/FiltFilt.cpp @@ -13,43 +13,22 @@ COPYING included with this distribution for more information. */ -#include #include "FiltFilt.h" ////////////////////////////////////////////////////////////////////// // Construction/Destruction ////////////////////////////////////////////////////////////////////// -FiltFilt::FiltFilt( FilterConfig Config ) +FiltFilt::FiltFilt(Filter::Parameters parameters) : + m_filter(parameters) { - m_filtScratchIn = NULL; - m_filtScratchOut = NULL; - m_ord = 0; - - initialise( Config ); + m_ord = m_filter.getOrder(); } FiltFilt::~FiltFilt() { - deInitialise(); } -void FiltFilt::initialise( FilterConfig Config ) -{ - m_ord = Config.ord; - m_filterConfig.ord = Config.ord; - m_filterConfig.ACoeffs = Config.ACoeffs; - m_filterConfig.BCoeffs = Config.BCoeffs; - - m_filter = new Filter( m_filterConfig ); -} - -void FiltFilt::deInitialise() -{ - delete m_filter; -} - - void FiltFilt::process(double *src, double *dst, unsigned int length) { unsigned int i; @@ -57,7 +36,6 @@ void FiltFilt::process(double *src, double *dst, unsigned int length) if (length == 0) return; if (length < 2) { - fprintf (stderr, "FiltFilt::process called for %d samples\n", length); for( i = 0; i < length; i++ ) { dst[i] = src [i]; } @@ -68,14 +46,13 @@ void FiltFilt::process(double *src, double *dst, unsigned int length) unsigned int nFact = 3 * ( nFilt - 1); unsigned int nExt = length + 2 * nFact; - m_filtScratchIn = new double[ nExt ]; - m_filtScratchOut = new double[ nExt ]; - + double *filtScratchIn = new double[ nExt ]; + double *filtScratchOut = new double[ nExt ]; for( i = 0; i< nExt; i++ ) { - m_filtScratchIn[ i ] = 0.0; - m_filtScratchOut[ i ] = 0.0; + filtScratchIn[ i ] = 0.0; + filtScratchOut[ i ] = 0.0; } // Edge transients reflection @@ -85,56 +62,56 @@ void FiltFilt::process(double *src, double *dst, unsigned int length) unsigned int index = 0; for( i = nFact; i > 0; i-- ) { - m_filtScratchIn[ index++ ] = sample0 - src[ i ]; + filtScratchIn[ index++ ] = sample0 - src[ i ]; } index = 0; for( i = 0; i < nFact && i + 2 < length; i++ ) { - m_filtScratchIn[ (nExt - nFact) + index++ ] = sampleN - src[ (length - 2) - i ]; + filtScratchIn[ (nExt - nFact) + index++ ] = sampleN - src[ (length - 2) - i ]; } for(; i < nFact; i++ ) { - m_filtScratchIn[ (nExt - nFact) + index++ ] = 0; + filtScratchIn[ (nExt - nFact) + index++ ] = 0; } index = 0; for( i = 0; i < length; i++ ) { - m_filtScratchIn[ i + nFact ] = src[ i ]; + filtScratchIn[ i + nFact ] = src[ i ]; } //////////////////////////////// // Do 0Ph filtering - m_filter->process( m_filtScratchIn, m_filtScratchOut, nExt); + m_filter.process( filtScratchIn, filtScratchOut, nExt); // reverse the series for FILTFILT for ( i = 0; i < nExt; i++) { - m_filtScratchIn[ i ] = m_filtScratchOut[ nExt - i - 1]; + filtScratchIn[ i ] = filtScratchOut[ nExt - i - 1]; } // do FILTER again - m_filter->process( m_filtScratchIn, m_filtScratchOut, nExt); + m_filter.process( filtScratchIn, filtScratchOut, nExt); // reverse the series back for ( i = 0; i < nExt; i++) { - m_filtScratchIn[ i ] = m_filtScratchOut[ nExt - i - 1 ]; + filtScratchIn[ i ] = filtScratchOut[ nExt - i - 1 ]; } for ( i = 0;i < nExt; i++) { - m_filtScratchOut[ i ] = m_filtScratchIn[ i ]; + filtScratchOut[ i ] = filtScratchIn[ i ]; } index = 0; for( i = 0; i < length; i++ ) { - dst[ index++ ] = m_filtScratchOut[ i + nFact ]; + dst[ index++ ] = filtScratchOut[ i + nFact ]; } - delete [] m_filtScratchIn; - delete [] m_filtScratchOut; + delete [] filtScratchIn; + delete [] filtScratchOut; } diff --git a/libs/qm-dsp/dsp/signalconditioning/FiltFilt.h b/libs/qm-dsp/dsp/signalconditioning/FiltFilt.h index e5a38124cf..9fc4395068 100644 --- a/libs/qm-dsp/dsp/signalconditioning/FiltFilt.h +++ b/libs/qm-dsp/dsp/signalconditioning/FiltFilt.h @@ -20,30 +20,21 @@ /** * Zero-phase digital filter, implemented by processing the data - * through a filter specified by the given FilterConfig structure (see + * through a filter specified by the given filter parameters (see * Filter) and then processing it again in reverse. */ class FiltFilt { public: - FiltFilt( FilterConfig Config ); + FiltFilt(Filter::Parameters); virtual ~FiltFilt(); void reset(); void process( double* src, double* dst, unsigned int length ); private: - void initialise( FilterConfig Config ); - void deInitialise(); - - unsigned int m_ord; - - Filter* m_filter; - - double* m_filtScratchIn; - double* m_filtScratchOut; - - FilterConfig m_filterConfig; + Filter m_filter; + int m_ord; }; #endif diff --git a/libs/qm-dsp/dsp/signalconditioning/Filter.cpp b/libs/qm-dsp/dsp/signalconditioning/Filter.cpp index fcc12e590a..e9523e229d 100644 --- a/libs/qm-dsp/dsp/signalconditioning/Filter.cpp +++ b/libs/qm-dsp/dsp/signalconditioning/Filter.cpp @@ -4,7 +4,6 @@ QM DSP Library Centre for Digital Music, Queen Mary, University of London. - This file 2005-2006 Christian Landone. This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as @@ -15,73 +14,112 @@ #include "Filter.h" -////////////////////////////////////////////////////////////////////// -// Construction/Destruction -////////////////////////////////////////////////////////////////////// +#include -Filter::Filter( FilterConfig Config ) +using namespace std; + +Filter::Filter(Parameters params) { - m_ord = 0; - m_outBuffer = NULL; - m_inBuffer = NULL; + if (params.a.empty()) { + m_fir = true; + if (params.b.empty()) { + throw logic_error("Filter must have at least one pair of coefficients"); + } + } else { + m_fir = false; + if (params.a.size() != params.b.size()) { + throw logic_error("Inconsistent numbers of filter coefficients"); + } + } + + m_sz = int(params.b.size()); + m_order = m_sz - 1; - initialise( Config ); + m_a = params.a; + m_b = params.b; + + // We keep some empty space at the start of the buffer, and + // encroach gradually into it as we add individual sample + // calculations at the start. Then when we run out of space, we + // move the buffer back to the end and begin again. This is + // significantly faster than moving the whole buffer along in + // 1-sample steps every time. + + m_offmax = 20; + m_offa = m_offmax; + m_offb = m_offmax; + + if (!m_fir) { + m_bufa.resize(m_order + m_offmax); + } + + m_bufb.resize(m_sz + m_offmax); } Filter::~Filter() { - deInitialise(); } -void Filter::initialise( FilterConfig Config ) +void +Filter::reset() { - m_ord = Config.ord; - m_ACoeffs = Config.ACoeffs; - m_BCoeffs = Config.BCoeffs; + m_offb = m_offmax; + m_offa = m_offmax; - m_inBuffer = new double[ m_ord + 1 ]; - m_outBuffer = new double[ m_ord + 1 ]; + if (!m_fir) { + m_bufa.assign(m_bufa.size(), 0.0); + } - reset(); + m_bufb.assign(m_bufb.size(), 0.0); } -void Filter::deInitialise() +void +Filter::process(const double *const __restrict__ in, + double *const __restrict__ out, + const int n) { - delete[] m_inBuffer; - delete[] m_outBuffer; + for (int s = 0; s < n; ++s) { + + if (m_offb > 0) --m_offb; + else { + for (int i = m_sz - 2; i >= 0; --i) { + m_bufb[i + m_offmax + 1] = m_bufb[i]; + } + m_offb = m_offmax; + } + m_bufb[m_offb] = in[s]; + + double b_sum = 0.0; + for (int i = 0; i < m_sz; ++i) { + b_sum += m_b[i] * m_bufb[i + m_offb]; + } + + double outval; + + if (m_fir) { + + outval = b_sum; + + } else { + + double a_sum = 0.0; + for (int i = 0; i < m_order; ++i) { + a_sum += m_a[i + 1] * m_bufa[i + m_offa]; + } + + outval = b_sum - a_sum; + + if (m_offa > 0) --m_offa; + else { + for (int i = m_order - 2; i >= 0; --i) { + m_bufa[i + m_offmax + 1] = m_bufa[i]; + } + m_offa = m_offmax; + } + m_bufa[m_offa] = outval; + } + + out[s] = outval; + } } -void Filter::reset() -{ - for( unsigned int i = 0; i < m_ord+1; i++ ){ m_inBuffer[ i ] = 0.0; } - for(unsigned int i = 0; i < m_ord+1; i++ ){ m_outBuffer[ i ] = 0.0; } -} - -void Filter::process( double *src, double *dst, unsigned int length ) -{ - unsigned int SP,i,j; - - double xin,xout; - - for (SP=0;SP -/** - * Filter specification. For a filter of order ord, the ACoeffs and - * BCoeffs arrays must point to ord+1 values each. ACoeffs provides - * the denominator and BCoeffs the numerator coefficients of the - * filter. - */ -struct FilterConfig{ - unsigned int ord; - double* ACoeffs; - double* BCoeffs; -}; - -/** - * Digital filter specified through FilterConfig structure. - */ -class Filter +class Filter { public: - Filter( FilterConfig Config ); - virtual ~Filter(); + struct Parameters { + std::vector a; + std::vector b; + }; + + /** + * Construct an IIR filter with numerators b and denominators + * a. The filter will have order b.size()-1. To make an FIR + * filter, leave the vector a in the param struct empty. + * Otherwise, a and b must have the same number of values. + */ + Filter(Parameters params); + + ~Filter(); void reset(); - void process( double *src, double *dst, unsigned int length ); + /** + * Filter the input sequence \arg in of length \arg n samples, and + * write the resulting \arg n samples into \arg out. There must be + * enough room in \arg out for \arg n samples to be written. + */ + void process(const double *const __restrict__ in, + double *const __restrict__ out, + const int n); + int getOrder() const { return m_order; } + private: - void initialise( FilterConfig Config ); - void deInitialise(); + int m_order; + int m_sz; + std::vector m_a; + std::vector m_b; + std::vector m_bufa; + std::vector m_bufb; + int m_offa; + int m_offb; + int m_offmax; + bool m_fir; - unsigned int m_ord; - - double* m_inBuffer; - double* m_outBuffer; - - double* m_ACoeffs; - double* m_BCoeffs; + Filter(const Filter &); // not supplied + Filter &operator=(const Filter &); // not supplied }; - + #endif diff --git a/libs/qm-dsp/dsp/tempotracking/TempoTrack.cpp b/libs/qm-dsp/dsp/tempotracking/TempoTrack.cpp index 389403edaa..7129b37966 100644 --- a/libs/qm-dsp/dsp/tempotracking/TempoTrack.cpp +++ b/libs/qm-dsp/dsp/tempotracking/TempoTrack.cpp @@ -70,10 +70,6 @@ void TempoTrack::initialise( TTParams Params ) m_tempoScratch = new double[ m_lagLength ]; m_smoothRCF = new double[ m_lagLength ]; - - unsigned int winPre = Params.WinT.pre; - unsigned int winPost = Params.WinT.post; - m_DFFramer.configure( m_winLength, m_lagLength ); m_DFPParams.length = m_winLength; @@ -120,9 +116,9 @@ void TempoTrack::deInitialise() } -void TempoTrack::createCombFilter(double* Filter, unsigned int winLength, unsigned int TSig, double beatLag) +void TempoTrack::createCombFilter(double* Filter, int winLength, int /* TSig */, double beatLag) { - unsigned int i; + int i; if( beatLag == 0 ) { @@ -147,15 +143,15 @@ double TempoTrack::tempoMM(double* ACF, double* weight, int tsig) double period = 0; double maxValRCF = 0.0; - unsigned int maxIndexRCF = 0; + int maxIndexRCF = 0; double* pdPeaks; - unsigned int maxIndexTemp; - double maxValTemp; - unsigned int count; + int maxIndexTemp; + double maxValTemp; + int count; - unsigned int numelem,i,j; + int numelem,i,j; int a, b; for( i = 0; i < m_lagLength; i++ ) @@ -476,7 +472,7 @@ void TempoTrack::constDetect( double* periodP, int currentIdx, int* flag ) } } -int TempoTrack::findMeter(double *ACF, unsigned int len, double period) +int TempoTrack::findMeter(double *ACF, int len, double period) { int i; int p = (int)MathUtilities::round( period ); @@ -491,7 +487,7 @@ int TempoTrack::findMeter(double *ACF, unsigned int len, double period) double temp4B = 0.0; double* dbf = new double[ len ]; int t = 0; - for( unsigned int u = 0; u < len; u++ ){ dbf[ u ] = 0.0; } + for( int u = 0; u < len; u++ ){ dbf[ u ] = 0.0; } if( (double)len < 6 * p + 2 ) { @@ -548,7 +544,7 @@ int TempoTrack::findMeter(double *ACF, unsigned int len, double period) return tsig; } -void TempoTrack::createPhaseExtractor(double *Filter, unsigned int winLength, double period, unsigned int fsp, unsigned int lastBeat) +void TempoTrack::createPhaseExtractor(double *Filter, int /* winLength */, double period, int fsp, int lastBeat) { int p = (int)MathUtilities::round( period ); int predictedOffset = 0; @@ -584,7 +580,7 @@ void TempoTrack::createPhaseExtractor(double *Filter, unsigned int winLength, do double sigma = (double)p/8; double PhaseMin = 0.0; double PhaseMax = 0.0; - unsigned int scratchLength = p*2; + int scratchLength = p*2; double temp = 0.0; for( int i = 0; i < scratchLength; i++ ) @@ -604,7 +600,7 @@ void TempoTrack::createPhaseExtractor(double *Filter, unsigned int winLength, do std::cerr << "predictedOffset = " << predictedOffset << std::endl; #endif - unsigned int index = 0; + int index = 0; for (int i = p - ( predictedOffset - 1); i < p + ( p - predictedOffset) + 1; i++) { #ifdef DEBUG_TEMPO_TRACK @@ -624,7 +620,7 @@ void TempoTrack::createPhaseExtractor(double *Filter, unsigned int winLength, do delete [] phaseScratch; } -int TempoTrack::phaseMM(double *DF, double *weighting, unsigned int winLength, double period) +int TempoTrack::phaseMM(double *DF, double *weighting, int winLength, double period) { int alignment = 0; int p = (int)MathUtilities::round( period ); @@ -667,7 +663,7 @@ int TempoTrack::phaseMM(double *DF, double *weighting, unsigned int winLength, d return alignment; } -int TempoTrack::beatPredict(unsigned int FSP0, double alignment, double period, unsigned int step ) +int TempoTrack::beatPredict(int FSP0, double alignment, double period, int step ) { int beat = 0; @@ -712,39 +708,39 @@ vector TempoTrack::process( vector DF, causalDF = DF; //Prepare Causal Extension DFData - unsigned int DFCLength = m_dataLength + m_winLength; +// int DFCLength = m_dataLength + m_winLength; - for( unsigned int j = 0; j < m_winLength; j++ ) + for( int j = 0; j < m_winLength; j++ ) { causalDF.push_back( 0 ); } double* RW = new double[ m_lagLength ]; - for( unsigned int clear = 0; clear < m_lagLength; clear++){ RW[ clear ] = 0.0;} + for (int clear = 0; clear < m_lagLength; clear++){ RW[ clear ] = 0.0;} double* GW = new double[ m_lagLength ]; - for(unsigned int clear = 0; clear < m_lagLength; clear++){ GW[ clear ] = 0.0;} + for (int clear = 0; clear < m_lagLength; clear++){ GW[ clear ] = 0.0;} double* PW = new double[ m_lagLength ]; - for(unsigned clear = 0; clear < m_lagLength; clear++){ PW[ clear ] = 0.0;} + for(int clear = 0; clear < m_lagLength; clear++){ PW[ clear ] = 0.0;} m_DFFramer.setSource( &causalDF[0], m_dataLength ); - unsigned int TTFrames = m_DFFramer.getMaxNoFrames(); + int TTFrames = m_DFFramer.getMaxNoFrames(); #ifdef DEBUG_TEMPO_TRACK std::cerr << "TTFrames = " << TTFrames << std::endl; #endif double* periodP = new double[ TTFrames ]; - for(unsigned clear = 0; clear < TTFrames; clear++){ periodP[ clear ] = 0.0;} + for(int clear = 0; clear < TTFrames; clear++){ periodP[ clear ] = 0.0;} double* periodG = new double[ TTFrames ]; - for(unsigned clear = 0; clear < TTFrames; clear++){ periodG[ clear ] = 0.0;} + for(int clear = 0; clear < TTFrames; clear++){ periodG[ clear ] = 0.0;} double* alignment = new double[ TTFrames ]; - for(unsigned clear = 0; clear < TTFrames; clear++){ alignment[ clear ] = 0.0;} + for(int clear = 0; clear < TTFrames; clear++){ alignment[ clear ] = 0.0;} m_beats.clear(); @@ -752,7 +748,7 @@ vector TempoTrack::process( vector DF, int TTLoopIndex = 0; - for( unsigned int i = 0; i < TTFrames; i++ ) + for( int i = 0; i < TTFrames; i++ ) { m_DFFramer.getFrame( m_rawDFFrame ); diff --git a/libs/qm-dsp/dsp/tempotracking/TempoTrack.h b/libs/qm-dsp/dsp/tempotracking/TempoTrack.h index 0c315717ba..72c2f756ef 100644 --- a/libs/qm-dsp/dsp/tempotracking/TempoTrack.h +++ b/libs/qm-dsp/dsp/tempotracking/TempoTrack.h @@ -30,16 +30,16 @@ using std::vector; struct WinThresh { - unsigned int pre; - unsigned int post; + int pre; + int post; }; struct TTParams { - unsigned int winLength; //Analysis window length - unsigned int lagLength; //Lag & Stride size - unsigned int alpha; //alpha-norm parameter - unsigned int LPOrd; // low-pass Filter order + int winLength; //Analysis window length + int lagLength; //Lag & Stride size + int alpha; //alpha-norm parameter + int LPOrd; // low-pass Filter order double* LPACoeffs; //low pass Filter den coefficients double* LPBCoeffs; //low pass Filter num coefficients WinThresh WinT;//window size in frames for adaptive thresholding [pre post]: @@ -59,22 +59,22 @@ private: void initialise( TTParams Params ); void deInitialise(); - int beatPredict( unsigned int FSP, double alignment, double period, unsigned int step); - int phaseMM( double* DF, double* weighting, unsigned int winLength, double period ); - void createPhaseExtractor( double* Filter, unsigned int winLength, double period, unsigned int fsp, unsigned int lastBeat ); - int findMeter( double* ACF, unsigned int len, double period ); + int beatPredict( int FSP, double alignment, double period, int step); + int phaseMM( double* DF, double* weighting, int winLength, double period ); + void createPhaseExtractor( double* Filter, int winLength, double period, int fsp, int lastBeat ); + int findMeter( double* ACF, int len, double period ); void constDetect( double* periodP, int currentIdx, int* flag ); void stepDetect( double* periodP, double* periodG, int currentIdx, int* flag ); - void createCombFilter( double* Filter, unsigned int winLength, unsigned int TSig, double beatLag ); + void createCombFilter( double* Filter, int winLength, int TSig, double beatLag ); double tempoMM( double* ACF, double* weight, int sig ); - unsigned int m_dataLength; - unsigned int m_winLength; - unsigned int m_lagLength; + int m_dataLength; + int m_winLength; + int m_lagLength; - double m_rayparam; - double m_sigma; - double m_DFWVNnorm; + double m_rayparam; + double m_sigma; + double m_DFWVNnorm; vector m_beats; // Vector of detected beats diff --git a/libs/qm-dsp/dsp/tonal/TCSgram.cpp b/libs/qm-dsp/dsp/tonal/TCSgram.cpp index 8840872576..8730dcf110 100644 --- a/libs/qm-dsp/dsp/tonal/TCSgram.cpp +++ b/libs/qm-dsp/dsp/tonal/TCSgram.cpp @@ -36,7 +36,7 @@ void TCSGram::getTCSVector(int iPosition, TCSVector& rTCSVector) const { if (iPosition < 0) rTCSVector = TCSVector(); - else if (iPosition >= m_VectorList.size()) + else if (iPosition >= int(m_VectorList.size())) rTCSVector = TCSVector(); else rTCSVector = m_VectorList[iPosition].second; diff --git a/libs/qm-dsp/dsp/tonal/TonalEstimator.h b/libs/qm-dsp/dsp/tonal/TonalEstimator.h index f555ef9e3b..b3c9c7b235 100644 --- a/libs/qm-dsp/dsp/tonal/TonalEstimator.h +++ b/libs/qm-dsp/dsp/tonal/TonalEstimator.h @@ -32,7 +32,7 @@ public: void printDebug() { - for (int i = 0; i < size(); i++) + for (int i = 0; i < int(size()); i++) { std::cout << (*this)[i] << ";"; } @@ -68,7 +68,7 @@ public: void printDebug() { - for (int i = 0; i < size(); i++) + for (int i = 0; i < int(size()); i++) { std::cout << (*this)[i] << ";"; } diff --git a/libs/qm-dsp/dsp/transforms/FFT.h b/libs/qm-dsp/dsp/transforms/FFT.h index c8956ac7f2..73af5616e7 100644 --- a/libs/qm-dsp/dsp/transforms/FFT.h +++ b/libs/qm-dsp/dsp/transforms/FFT.h @@ -4,6 +4,12 @@ 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 FFT_H diff --git a/libs/qm-dsp/dsp/wavelet/Wavelet.cpp b/libs/qm-dsp/dsp/wavelet/Wavelet.cpp index 764c84b24a..9c29b73cf5 100644 --- a/libs/qm-dsp/dsp/wavelet/Wavelet.cpp +++ b/libs/qm-dsp/dsp/wavelet/Wavelet.cpp @@ -1843,7 +1843,10 @@ Wavelet::createDecompositionFilters(Type wavelet, break; } - assert(flength == lpd.size()); - assert(flength == hpd.size()); + // avoid compiler warning for unused value if assert is not compiled in: + (void)flength; + + assert(flength == int(lpd.size())); + assert(flength == int(hpd.size())); } diff --git a/libs/qm-dsp/ext/kissfft/kiss_fft.c b/libs/qm-dsp/ext/kissfft/kiss_fft.c index 465d6c97a0..7824d34521 100644 --- a/libs/qm-dsp/ext/kissfft/kiss_fft.c +++ b/libs/qm-dsp/ext/kissfft/kiss_fft.c @@ -275,7 +275,8 @@ void kf_work( if (m==1) { do{ - *Fout = *f; + Fout->r = f->r; + Fout->i = f->i; f += fstride*in_stride; }while(++Fout != Fout_end ); }else{ diff --git a/libs/qm-dsp/gitrev.txt b/libs/qm-dsp/gitrev.txt new file mode 100644 index 0000000000..dde4075d05 --- /dev/null +++ b/libs/qm-dsp/gitrev.txt @@ -0,0 +1 @@ +qm-vamp-plugins-v1.7.1-20-g4d15479 diff --git a/libs/qm-dsp/maths/CosineDistance.cpp b/libs/qm-dsp/maths/CosineDistance.cpp index 13ab9ce0e8..4b06a7ad2a 100644 --- a/libs/qm-dsp/maths/CosineDistance.cpp +++ b/libs/qm-dsp/maths/CosineDistance.cpp @@ -34,7 +34,7 @@ double CosineDistance::distance(const vector &v1, } else { - for(int i=0; i #include +using namespace std; double MathUtilities::mod(double x, double y) { @@ -38,9 +39,9 @@ double MathUtilities::princarg(double ang) return ValOut; } -void MathUtilities::getAlphaNorm(const double *data, unsigned int len, unsigned int alpha, double* ANorm) +void MathUtilities::getAlphaNorm(const double *data, int len, int alpha, double* ANorm) { - unsigned int i; + int i; double temp = 0.0; double a=0.0; @@ -56,17 +57,16 @@ void MathUtilities::getAlphaNorm(const double *data, unsigned int len, unsigned *ANorm = a; } -double MathUtilities::getAlphaNorm( const std::vector &data, unsigned int alpha ) +double MathUtilities::getAlphaNorm( const vector &data, int alpha ) { - unsigned int i; - unsigned int len = data.size(); + int i; + int len = data.size(); double temp = 0.0; double a=0.0; for( i = 0; i < len; i++) { temp = data[ i ]; - a += ::pow( fabs(temp), double(alpha) ); } a /= ( double )len; @@ -84,13 +84,13 @@ double MathUtilities::round(double x) } } -double MathUtilities::median(const double *src, unsigned int len) +double MathUtilities::median(const double *src, int len) { if (len == 0) return 0; - std::vector scratch; + vector scratch; for (int i = 0; i < len; ++i) scratch.push_back(src[i]); - std::sort(scratch.begin(), scratch.end()); + sort(scratch.begin(), scratch.end()); int middle = len/2; if (len % 2 == 0) { @@ -100,9 +100,9 @@ double MathUtilities::median(const double *src, unsigned int len) } } -double MathUtilities::sum(const double *src, unsigned int len) +double MathUtilities::sum(const double *src, int len) { - unsigned int i ; + int i ; double retVal =0.0; for( i = 0; i < len; i++) @@ -113,7 +113,7 @@ double MathUtilities::sum(const double *src, unsigned int len) return retVal; } -double MathUtilities::mean(const double *src, unsigned int len) +double MathUtilities::mean(const double *src, int len) { double retVal =0.0; @@ -126,9 +126,9 @@ double MathUtilities::mean(const double *src, unsigned int len) return retVal; } -double MathUtilities::mean(const std::vector &src, - unsigned int start, - unsigned int count) +double MathUtilities::mean(const vector &src, + int start, + int count) { double sum = 0.; @@ -142,9 +142,9 @@ double MathUtilities::mean(const std::vector &src, return sum / count; } -void MathUtilities::getFrameMinMax(const double *data, unsigned int len, double *min, double *max) +void MathUtilities::getFrameMinMax(const double *data, int len, double *min, double *max) { - unsigned int i; + int i; double temp = 0.0; if (len == 0) { @@ -171,10 +171,10 @@ void MathUtilities::getFrameMinMax(const double *data, unsigned int len, double } } -int MathUtilities::getMax( double* pData, unsigned int Length, double* pMax ) +int MathUtilities::getMax( double* pData, int Length, double* pMax ) { - unsigned int index = 0; - unsigned int i; + int index = 0; + int i; double temp = 0.0; double max = pData[0]; @@ -197,15 +197,15 @@ int MathUtilities::getMax( double* pData, unsigned int Length, double* pMax ) return index; } -int MathUtilities::getMax( const std::vector & data, double* pMax ) +int MathUtilities::getMax( const vector & data, double* pMax ) { - unsigned int index = 0; - unsigned int i; + int index = 0; + int i; double temp = 0.0; double max = data[0]; - for( i = 0; i < data.size(); i++) + for( i = 0; i < int(data.size()); i++) { temp = data[ i ]; @@ -286,7 +286,7 @@ void MathUtilities::normalise(double *data, int length, NormaliseType type) } } -void MathUtilities::normalise(std::vector &data, NormaliseType type) +void MathUtilities::normalise(vector &data, NormaliseType type) { switch (type) { @@ -317,20 +317,46 @@ void MathUtilities::normalise(std::vector &data, NormaliseType type) } } -void MathUtilities::adaptiveThreshold(std::vector &data) +double MathUtilities::getLpNorm(const vector &data, int p) +{ + double tot = 0.0; + for (int i = 0; i < int(data.size()); ++i) { + tot += abs(pow(data[i], p)); + } + return pow(tot, 1.0 / p); +} + +vector MathUtilities::normaliseLp(const vector &data, + int p, + double threshold) +{ + int n = int(data.size()); + if (n == 0 || p == 0) return data; + double norm = getLpNorm(data, p); + if (norm < threshold) { + return vector(n, 1.0 / pow(n, 1.0 / p)); // unit vector + } + vector out(n); + for (int i = 0; i < n; ++i) { + out[i] = data[i] / norm; + } + return out; +} + +void MathUtilities::adaptiveThreshold(vector &data) { int sz = int(data.size()); if (sz == 0) return; - std::vector smoothed(sz); + vector smoothed(sz); int p_pre = 8; int p_post = 7; for (int i = 0; i < sz; ++i) { - int first = std::max(0, i - p_pre); - int last = std::min(sz - 1, i + p_post); + int first = max(0, i - p_pre); + int last = min(sz - 1, i + p_post); smoothed[i] = mean(data, first, last - first + 1); } diff --git a/libs/qm-dsp/maths/MathUtilities.h b/libs/qm-dsp/maths/MathUtilities.h index fac710af9a..ec86efcf2f 100644 --- a/libs/qm-dsp/maths/MathUtilities.h +++ b/libs/qm-dsp/maths/MathUtilities.h @@ -35,32 +35,32 @@ public: * Return through min and max pointers the highest and lowest * values in the given array of the given length. */ - static void getFrameMinMax( const double* data, unsigned int len, double* min, double* max ); + static void getFrameMinMax( const double* data, int len, double* min, double* max ); /** * Return the mean of the given array of the given length. */ - static double mean( const double* src, unsigned int len ); + static double mean( const double* src, int len ); /** * Return the mean of the subset of the given vector identified by * start and count. */ static double mean( const std::vector &data, - unsigned int start, unsigned int count ); + int start, int count ); /** * Return the sum of the values in the given array of the given * length. */ - static double sum( const double* src, unsigned int len ); + static double sum( const double* src, int len ); /** * Return the median of the values in the given array of the given * length. If the array is even in length, the returned value will * be half-way between the two values adjacent to median. */ - static double median( const double* src, unsigned int len ); + static double median( const double* src, int len ); /** * The principle argument function. Map the phase angle ang into @@ -73,14 +73,21 @@ public: */ static double mod( double x, double y); - static void getAlphaNorm(const double *data, unsigned int len, unsigned int alpha, double* ANorm); - static double getAlphaNorm(const std::vector &data, unsigned int alpha ); + /** + * The alpha norm is the alpha'th root of the mean alpha'th power + * magnitude. For example if alpha = 2 this corresponds to the RMS + * of the input data, and when alpha = 1 this is the mean + * magnitude. + */ + static void getAlphaNorm(const double *data, int len, int alpha, double* ANorm); - static void circShift( double* data, int length, int shift); - - static int getMax( double* data, unsigned int length, double* max = 0 ); - static int getMax( const std::vector &data, double* max = 0 ); - static int compareInt(const void * a, const void * b); + /** + * The alpha norm is the alpha'th root of the mean alpha'th power + * magnitude. For example if alpha = 2 this corresponds to the RMS + * of the input data, and when alpha = 1 this is the mean + * magnitude. + */ + static double getAlphaNorm(const std::vector &data, int alpha ); enum NormaliseType { NormaliseNone, @@ -94,12 +101,35 @@ public: static void normalise(std::vector &data, NormaliseType n = NormaliseUnitMax); + /** + * Calculate the L^p norm of a vector. Equivalent to MATLAB's + * norm(data, p). + */ + static double getLpNorm(const std::vector &data, + int p); + + /** + * Normalise a vector by dividing through by its L^p norm. If the + * norm is below the given threshold, the unit vector for that + * norm is returned. p may be 0, in which case no normalisation + * happens and the data is returned unchanged. + */ + static std::vector normaliseLp(const std::vector &data, + int p, + double threshold = 1e-6); + /** * Threshold the input/output vector data against a moving-mean * average filter. */ static void adaptiveThreshold(std::vector &data); + static void circShift( double* data, int length, int shift); + + static int getMax( double* data, int length, double* max = 0 ); + static int getMax( const std::vector &data, double* max = 0 ); + static int compareInt(const void * a, const void * b); + /** * Return true if x is 2^n for some integer n >= 0. */ diff --git a/libs/qm-dsp/maths/Polyfit.h b/libs/qm-dsp/maths/Polyfit.h index 5cf97d9df0..3405b5cb70 100644 --- a/libs/qm-dsp/maths/Polyfit.h +++ b/libs/qm-dsp/maths/Polyfit.h @@ -80,22 +80,36 @@ private: // some utility functions -namespace NSUtility +struct NSUtility { - inline void swap(double &a, double &b) {double t = a; a = b; b = t;} - void zeroise(vector &array, int n); - void zeroise(vector &array, int n); - void zeroise(vector > &matrix, int m, int n); - void zeroise(vector > &matrix, int m, int n); - inline double sqr(const double &x) {return x * x;} + static void swap(double &a, double &b) {double t = a; a = b; b = t;} + // fills a vector with zeros. + static void zeroise(vector &array, int n) { + array.clear(); + for(int j = 0; j < n; ++j) array.push_back(0); + } + // fills a vector with zeros. + static void zeroise(vector &array, int n) { + array.clear(); + for(int j = 0; j < n; ++j) array.push_back(0); + } + // fills a (m by n) matrix with zeros. + static void zeroise(vector > &matrix, int m, int n) { + vector zero; + zeroise(zero, n); + matrix.clear(); + for(int j = 0; j < m; ++j) matrix.push_back(zero); + } + // fills a (m by n) matrix with zeros. + static void zeroise(vector > &matrix, int m, int n) { + vector zero; + zeroise(zero, n); + matrix.clear(); + for(int j = 0; j < m; ++j) matrix.push_back(zero); + } + static double sqr(const double &x) {return x * x;} }; -//--------------------------------------------------------------------------- -// Implementation -//--------------------------------------------------------------------------- -using namespace NSUtility; -//------------------------------------------------------------------------------------------ - // main PolyFit routine @@ -113,9 +127,9 @@ double TPolyFit::PolyFit2 (const vector &x, const int npoints(x.size()); const int nterms(coefs.size()); double correl_coef; - zeroise(g, nterms); - zeroise(a, nterms, nterms); - zeroise(xmatr, npoints, nterms); + NSUtility::zeroise(g, nterms); + NSUtility::zeroise(a, nterms, nterms); + NSUtility::zeroise(xmatr, npoints, nterms); if (nterms < 1) { std::cerr << "ERROR: PolyFit called with less than one term" << std::endl; return 0; @@ -148,13 +162,13 @@ double TPolyFit::PolyFit2 (const vector &x, yc = 0.0; for(j = 0; j < nterms; ++j) yc += coefs [j] * xmatr [i][j]; - srs += sqr (yc - yi); + srs += NSUtility::sqr (yc - yi); sum_y += yi; sum_y2 += yi * yi; } // If all Y values are the same, avoid dividing by zero - correl_coef = sum_y2 - sqr (sum_y) / npoints; + correl_coef = sum_y2 - NSUtility::sqr (sum_y) / npoints; // Either return 0 or the correct value of correlation coefficient if (correl_coef != 0) correl_coef = srs / correl_coef; @@ -229,8 +243,8 @@ bool TPolyFit::GaussJordan (Matrix &b, vector >index; Matrix w; - zeroise(w, ncol, ncol); - zeroise(index, ncol, 3); + NSUtility::zeroise(w, ncol, ncol); + NSUtility::zeroise(index, ncol, 3); if(!GaussJordan2(b, y, w, index)) return false; @@ -278,7 +292,7 @@ bool TPolyFit::GaussJordan2(Matrix &b, double big, t; double pivot; double determ; - int irow, icol; + int irow = 0, icol = 0; int ncol(b.size()); int nv = 1; // single constant vector for(int i = 0; i < ncol; ++i) @@ -355,53 +369,6 @@ bool TPolyFit::GaussJordan2(Matrix &b, } // { i-loop } return true; } -//---------------------------------------------------------------------------------------------- - -//------------------------------------------------------------------------------------ - -// Utility functions -//-------------------------------------------------------------------------- - -// fills a vector with zeros. -void NSUtility::zeroise(vector &array, int n) -{ - array.clear(); - for(int j = 0; j < n; ++j) - array.push_back(0); -} -//-------------------------------------------------------------------------- - -// fills a vector with zeros. -void NSUtility::zeroise(vector &array, int n) -{ - array.clear(); - for(int j = 0; j < n; ++j) - array.push_back(0); -} -//-------------------------------------------------------------------------- - -// fills a (m by n) matrix with zeros. -void NSUtility::zeroise(vector > &matrix, int m, int n) -{ - vector zero; - zeroise(zero, n); - matrix.clear(); - for(int j = 0; j < m; ++j) - matrix.push_back(zero); -} -//-------------------------------------------------------------------------- - -// fills a (m by n) matrix with zeros. -void NSUtility::zeroise(vector > &matrix, int m, int n) -{ - vector zero; - zeroise(zero, n); - matrix.clear(); - for(int j = 0; j < m; ++j) - matrix.push_back(zero); -} -//-------------------------------------------------------------------------- - #endif diff --git a/libs/qm-dsp/maths/pca/pca.c b/libs/qm-dsp/maths/pca/pca.c index 9dadb8687d..1a7bfdd549 100644 --- a/libs/qm-dsp/maths/pca/pca.c +++ b/libs/qm-dsp/maths/pca/pca.c @@ -252,7 +252,7 @@ void tqli(double* d, double* e, int n, double** z) void pca_project(double** data, int n, int m, int ncomponents) { int i, j, k, k2; - double **symmat, **symmat2, *evals, *interm; + double **symmat, /* **symmat2, */ *evals, *interm; //TODO: assert ncomponents < m