Mercurial > hg > qm-dsp
changeset 195:4e57fb34d9e6
Merge
| author | Chris Cannam |
|---|---|
| date | Wed, 07 Oct 2015 11:15:13 +0100 |
| parents | 26daede606a8 (diff) 786d446fe22b (current diff) |
| children | be1fc3a6b901 |
| files | |
| diffstat | 17 files changed, 815 insertions(+), 168 deletions(-) [+] |
line wrap: on
line diff
--- a/.hgignore Wed Oct 07 09:33:28 2015 +0100 +++ b/.hgignore Wed Oct 07 11:15:13 2015 +0100 @@ -5,5 +5,7 @@ *.a doc/html/ tests/test-* +tests/test.log *.rej +ext/uncertain/
--- a/build/general/Makefile.inc Wed Oct 07 09:33:28 2015 +0100 +++ b/build/general/Makefile.inc Wed Oct 07 11:15:13 2015 +0100 @@ -38,6 +38,7 @@ dsp/tonal/ChangeDetectionFunction.h \ dsp/tonal/TCSgram.h \ dsp/tonal/TonalEstimator.h \ + dsp/transforms/DCT.h \ dsp/transforms/FFT.h \ dsp/wavelet/Wavelet.h \ hmm/hmm.h \ @@ -83,6 +84,7 @@ dsp/tonal/ChangeDetectionFunction.cpp \ dsp/tonal/TCSgram.cpp \ dsp/tonal/TonalEstimator.cpp \ + dsp/transforms/DCT.cpp \ dsp/transforms/FFT.cpp \ dsp/wavelet/Wavelet.cpp \ hmm/hmm.c \
--- a/dsp/signalconditioning/DFProcess.cpp Wed Oct 07 09:33:28 2015 +0100 +++ b/dsp/signalconditioning/DFProcess.cpp Wed Oct 07 11:15:13 2015 +0100 @@ -59,13 +59,11 @@ 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<double>(Config.LPACoeffs, Config.LPACoeffs + Config.LPOrd + 1); + params.b = std::vector<double>(Config.LPBCoeffs, Config.LPBCoeffs + Config.LPOrd + 1); + + m_FiltFilt = new FiltFilt(params); //add delta threshold m_delta = Config.delta;
--- a/dsp/signalconditioning/DFProcess.h Wed Oct 07 09:33:28 2015 +0100 +++ b/dsp/signalconditioning/DFProcess.h Wed Oct 07 11:15:13 2015 +0100 @@ -81,8 +81,6 @@ double* m_filtScratchIn; double* m_filtScratchOut; - FilterConfig m_FilterConfigParams; - FiltFilt* m_FiltFilt; bool m_isMedianPositive;
--- a/dsp/signalconditioning/FiltFilt.cpp Wed Oct 07 09:33:28 2015 +0100 +++ b/dsp/signalconditioning/FiltFilt.cpp Wed Oct 07 11:15:13 2015 +0100 @@ -19,36 +19,16 @@ // 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; @@ -59,14 +39,13 @@ 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 @@ -76,51 +55,51 @@ 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++ ) { - m_filtScratchIn[ (nExt - nFact) + index++ ] = sampleN - src[ (length - 2) - i ]; + filtScratchIn[ (nExt - nFact) + index++ ] = sampleN - src[ (length - 2) - i ]; } 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; }
--- a/dsp/signalconditioning/FiltFilt.h Wed Oct 07 09:33:28 2015 +0100 +++ b/dsp/signalconditioning/FiltFilt.h Wed Oct 07 11:15:13 2015 +0100 @@ -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
--- a/dsp/signalconditioning/Filter.cpp Wed Oct 07 09:33:28 2015 +0100 +++ b/dsp/signalconditioning/Filter.cpp Wed Oct 07 11:15:13 2015 +0100 @@ -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 <stdexcept> -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<length;SP++) - { - xin=src[SP]; - /* move buffer */ - for ( i = 0; i < m_ord; i++) {m_inBuffer[ m_ord - i ]=m_inBuffer[ m_ord - i - 1 ];} - m_inBuffer[0]=xin; - - xout=0.0; - for (j=0;j< m_ord + 1; j++) - xout = xout + m_BCoeffs[ j ] * m_inBuffer[ j ]; - for (j = 0; j < m_ord; j++) - xout= xout - m_ACoeffs[ j + 1 ] * m_outBuffer[ j ]; - - dst[ SP ] = xout; - for ( i = 0; i < m_ord - 1; i++ ) { m_outBuffer[ m_ord - i - 1 ] = m_outBuffer[ m_ord - i - 2 ];} - m_outBuffer[0]=xout; - - } /* end of SP loop */ -} - - -
--- a/dsp/signalconditioning/Filter.h Wed Oct 07 09:33:28 2015 +0100 +++ b/dsp/signalconditioning/Filter.h Wed Oct 07 11:15:13 2015 +0100 @@ -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 @@ -16,46 +15,53 @@ #ifndef FILTER_H #define FILTER_H -#ifndef NULL -#define NULL 0 -#endif +#include <vector> -/** - * 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<double> a; + std::vector<double> 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<double> m_a; + std::vector<double> m_b; + std::vector<double> m_bufa; + std::vector<double> 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
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/transforms/DCT.cpp Wed Oct 07 11:15:13 2015 +0100 @@ -0,0 +1,91 @@ +/* -*- 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 "DCT.h" + +#include <cmath> + +DCT::DCT(int n) : + m_n(n), + m_scaled(n, 0.0), + m_time2(n*4, 0.0), + m_freq2r(n*4, 0.0), + m_freq2i(n*4, 0.0), + m_fft(n*4) +{ + m_scale = m_n * sqrt(2.0 / m_n); +} + +DCT::~DCT() +{ +} + +void +DCT::forward(const double *in, double *out) +{ + for (int i = 0; i < m_n; ++i) { + m_time2[i*2 + 1] = in[i]; + m_time2[m_n*4 - i*2 - 1] = in[i]; + } + + m_fft.forward(m_time2.data(), m_freq2r.data(), m_freq2i.data()); + + for (int i = 0; i < m_n; ++i) { + out[i] = m_freq2r[i]; + } +} + +void +DCT::forwardUnitary(const double *in, double *out) +{ + forward(in, out); + for (int i = 0; i < m_n; ++i) { + out[i] /= m_scale; + } + out[0] /= sqrt(2.0); +} + +void +DCT::inverse(const double *in, double *out) +{ + for (int i = 0; i < m_n; ++i) { + m_freq2r[i] = in[i]; + } + for (int i = 0; i < m_n; ++i) { + m_freq2r[m_n*2 - i] = -in[i]; + } + m_freq2r[m_n] = 0.0; + + for (int i = 0; i <= m_n*2; ++i) { + m_freq2i[i] = 0.0; + } + + m_fft.inverse(m_freq2r.data(), m_freq2i.data(), m_time2.data()); + + for (int i = 0; i < m_n; ++i) { + out[i] = m_time2[i*2 + 1]; + } +} + +void +DCT::inverseUnitary(const double *in, double *out) +{ + for (int i = 0; i < m_n; ++i) { + m_scaled[i] = in[i] * m_scale; + } + m_scaled[0] *= sqrt(2.0); + inverse(m_scaled.data(), out); +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/transforms/DCT.h Wed Oct 07 11:15:13 2015 +0100 @@ -0,0 +1,85 @@ +/* -*- 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 DCT_H +#define DCT_H + +#include "FFT.h" + +#include <vector> + +class DCT +{ +public: + /** + * Construct a DCT object to calculate the Discrete Cosine + * Transform given input of size n samples. The transform is + * implemented using an FFT of size 4n, for simplicity. + */ + DCT(int n); + + ~DCT(); + + /** + * Carry out a type-II DCT of size n, where n is the value + * provided to the constructor above. + * + * The in and out pointers must point to (enough space for) n + * values. + */ + void forward(const double *in, double *out); + + /** + * Carry out a type-II unitary DCT of size n, where n is the value + * provided to the constructor above. This is a scaled version of + * the type-II DCT corresponding to a transform with an orthogonal + * matrix. This is the transform implemented by the dct() function + * in MATLAB. + * + * The in and out pointers must point to (enough space for) n + * values. + */ + void forwardUnitary(const double *in, double *out); + + /** + * Carry out a type-III (inverse) DCT of size n, where n is the + * value provided to the constructor above. + * + * The in and out pointers must point to (enough space for) n + * values. + */ + void inverse(const double *in, double *out); + + /** + * Carry out a type-III (inverse) unitary DCT of size n, where n + * is the value provided to the constructor above. This is the + * inverse of forwardUnitary(). + * + * The in and out pointers must point to (enough space for) n + * values. + */ + void inverseUnitary(const double *in, double *out); + +private: + int m_n; + double m_scale; + std::vector<double> m_scaled; + std::vector<double> m_time2; + std::vector<double> m_freq2r; + std::vector<double> m_freq2i; + FFTReal m_fft; +}; + +#endif
--- a/dsp/transforms/FFT.h Wed Oct 07 09:33:28 2015 +0100 +++ b/dsp/transforms/FFT.h Wed Oct 07 11:15:13 2015 +0100 @@ -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
--- a/maths/MathUtilities.cpp Wed Oct 07 09:33:28 2015 +0100 +++ b/maths/MathUtilities.cpp Wed Oct 07 11:15:13 2015 +0100 @@ -20,6 +20,7 @@ #include <vector> #include <cmath> +using namespace std; double MathUtilities::mod(double x, double y) { @@ -56,7 +57,7 @@ *ANorm = a; } -double MathUtilities::getAlphaNorm( const std::vector <double> &data, int alpha ) +double MathUtilities::getAlphaNorm( const vector <double> &data, int alpha ) { int i; int len = data.size(); @@ -66,7 +67,6 @@ for( i = 0; i < len; i++) { temp = data[ i ]; - a += ::pow( fabs(temp), double(alpha) ); } a /= ( double )len; @@ -88,9 +88,9 @@ { if (len == 0) return 0; - std::vector<double> scratch; + vector<double> 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) { @@ -126,7 +126,7 @@ return retVal; } -double MathUtilities::mean(const std::vector<double> &src, +double MathUtilities::mean(const vector<double> &src, int start, int count) { @@ -197,7 +197,7 @@ return index; } -int MathUtilities::getMax( const std::vector<double> & data, double* pMax ) +int MathUtilities::getMax( const vector<double> & data, double* pMax ) { int index = 0; int i; @@ -286,7 +286,7 @@ } } -void MathUtilities::normalise(std::vector<double> &data, NormaliseType type) +void MathUtilities::normalise(vector<double> &data, NormaliseType type) { switch (type) { @@ -317,20 +317,46 @@ } } -void MathUtilities::adaptiveThreshold(std::vector<double> &data) +double MathUtilities::getLpNorm(const vector<double> &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<double> MathUtilities::normaliseLp(const vector<double> &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<double>(n, 1.0 / pow(n, 1.0 / p)); // unit vector + } + vector<double> out(n); + for (int i = 0; i < n; ++i) { + out[i] = data[i] / norm; + } + return out; +} + +void MathUtilities::adaptiveThreshold(vector<double> &data) { int sz = int(data.size()); if (sz == 0) return; - std::vector<double> smoothed(sz); + vector<double> 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); }
--- a/maths/MathUtilities.h Wed Oct 07 09:33:28 2015 +0100 +++ b/maths/MathUtilities.h Wed Oct 07 11:15:13 2015 +0100 @@ -73,15 +73,22 @@ */ static double mod( double x, double y); + /** + * 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); + + /** + * 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 <double> &data, int alpha ); - 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<double> &data, double* max = 0 ); - static int compareInt(const void * a, const void * b); - enum NormaliseType { NormaliseNone, NormaliseUnitSum, @@ -95,11 +102,34 @@ NormaliseType n = NormaliseUnitMax); /** + * Calculate the L^p norm of a vector. Equivalent to MATLAB's + * norm(data, p). + */ + static double getLpNorm(const std::vector<double> &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<double> normaliseLp(const std::vector<double> &data, + int p, + double threshold = 1e-6); + + /** * Threshold the input/output vector data against a moving-mean * average filter. */ static void adaptiveThreshold(std::vector<double> &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<double> &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. */
--- a/tests/Makefile Wed Oct 07 09:33:28 2015 +0100 +++ b/tests/Makefile Wed Oct 07 11:15:13 2015 +0100 @@ -1,11 +1,11 @@ CFLAGS := -I.. $(CFLAGS) -CXXFLAGS := -I.. -Wall -Wextra -g $(CXXFLAGS) +CXXFLAGS := -I.. -Wall -Wextra -std=c++11 -g $(CXXFLAGS) LDFLAGS := $(LDFLAGS) -lboost_unit_test_framework -lpthread LIBS := ../libqm-dsp.a -TESTS := test-mathutilities test-window test-fft test-pvoc test-resampler test-medianfilter +TESTS := test-mathutilities test-window test-filter test-fft test-dct test-pvoc test-resampler test-medianfilter VG := valgrind -q @@ -21,9 +21,15 @@ test-window: TestWindow.o $(LIBS) $(CXX) -o $@ $^ $(LDFLAGS) +test-filter: TestFilter.o $(LIBS) + $(CXX) -o $@ $^ $(LDFLAGS) + test-fft: TestFFT.o $(LIBS) $(CXX) -o $@ $^ $(LDFLAGS) +test-dct: TestDCT.o $(LIBS) + $(CXX) -o $@ $^ $(LDFLAGS) + test-pvoc: TestPhaseVocoder.o $(LIBS) $(CXX) -o $@ $^ $(LDFLAGS) @@ -33,7 +39,9 @@ TestMathUtilities.o: $(LIBS) TestMedianFilter.o: $(LIBS) TestWindow.o: $(LIBS) +TestFilter.o: $(LIBS) TestFFT.o: $(LIBS) +TestDCT.o: $(LIBS) TestPhaseVocoder.o: $(LIBS) TestResampler.o: $(LIBS)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/TestDCT.cpp Wed Oct 07 11:15:13 2015 +0100 @@ -0,0 +1,196 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ + +#include "dsp/transforms/DCT.h" + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MAIN + +#include <boost/test/unit_test.hpp> + +#include <stdexcept> + +using namespace std; + +BOOST_AUTO_TEST_SUITE(TestDCT) + +BOOST_AUTO_TEST_CASE(forwardArrayBounds) +{ + // initialise bins to something recognisable, so we can tell if + // they haven't been written; and allocate the inputs on the heap + // so that, if running under valgrind, we get warnings about + // overruns + double *in = new double[4]; + in[0] = 1; + in[1] = 1; + in[2] = -1; + in[3] = -1; + double out[] = { 999, 999, 999, 999, 999, 999 }; + DCT(4).forward(in, out+1); + // And check we haven't overrun the arrays + BOOST_CHECK_EQUAL(out[0], 999.0); + BOOST_CHECK_EQUAL(out[5], 999.0); + delete[] in; +} + +BOOST_AUTO_TEST_CASE(u_forwardArrayBounds) +{ + // initialise bins to something recognisable, so we can tell if + // they haven't been written; and allocate the inputs on the heap + // so that, if running under valgrind, we get warnings about + // overruns + double *in = new double[4]; + in[0] = 1; + in[1] = 1; + in[2] = -1; + in[3] = -1; + double out[] = { 999, 999, 999, 999, 999, 999 }; + DCT(4).forwardUnitary(in, out+1); + // And check we haven't overrun the arrays + BOOST_CHECK_EQUAL(out[0], 999.0); + BOOST_CHECK_EQUAL(out[5], 999.0); + delete[] in; +} + +BOOST_AUTO_TEST_CASE(inverseArrayBounds) +{ + // initialise bins to something recognisable, so we can tell if + // they haven't been written; and allocate the inputs on the heap + // so that, if running under valgrind, we get warnings about + // overruns + double *in = new double[4]; + in[0] = 1; + in[1] = 1; + in[2] = -1; + in[3] = -1; + double out[] = { 999, 999, 999, 999, 999, 999 }; + DCT(4).inverse(in, out+1); + // And check we haven't overrun the arrays + BOOST_CHECK_EQUAL(out[0], 999.0); + BOOST_CHECK_EQUAL(out[5], 999.0); + delete[] in; + +} + +BOOST_AUTO_TEST_CASE(u_inverseArrayBounds) +{ + // initialise bins to something recognisable, so we can tell if + // they haven't been written; and allocate the inputs on the heap + // so that, if running under valgrind, we get warnings about + // overruns + double *in = new double[4]; + in[0] = 1; + in[1] = 1; + in[2] = -1; + in[3] = -1; + double out[] = { 999, 999, 999, 999, 999, 999 }; + DCT(4).inverseUnitary(in, out+1); + // And check we haven't overrun the arrays + BOOST_CHECK_EQUAL(out[0], 999.0); + BOOST_CHECK_EQUAL(out[5], 999.0); + delete[] in; +} + +BOOST_AUTO_TEST_CASE(nonU4) +{ + int n = 4; + vector<double> in { 1, 2, 3, 5 }; + vector<double> out(n); + vector<double> inv(n); + vector<double> expected { 22.0, -8.1564, 1.4142, -1.2137 }; + + DCT d(n); + + // our expected values are very approximate + double thresh = 1e-4; + + d.forward(in.data(), out.data()); + + for (int i = 0; i < n; ++i) { + BOOST_CHECK_SMALL(expected[i] - out[i], thresh); + } + + d.inverse(out.data(), inv.data()); + + for (int i = 0; i < n; ++i) { + BOOST_CHECK_SMALL(in[i] - inv[i], thresh); + } + + // do it again, just in case some internal state was modified in inverse + + d.forward(in.data(), out.data()); + + for (int i = 0; i < n; ++i) { + BOOST_CHECK_SMALL(expected[i] - out[i], thresh); + } +} + +BOOST_AUTO_TEST_CASE(u4) +{ + int n = 4; + vector<double> in { 1, 2, 3, 5 }; + vector<double> out(n); + vector<double> inv(n); + vector<double> expected { 5.5, -2.8837, 0.5, -0.4291 }; + + DCT d(n); + + // our expected values are very approximate + double thresh = 1e-4; + + d.forwardUnitary(in.data(), out.data()); + + for (int i = 0; i < n; ++i) { + BOOST_CHECK_SMALL(expected[i] - out[i], thresh); + } + + d.inverseUnitary(out.data(), inv.data()); + + for (int i = 0; i < n; ++i) { + BOOST_CHECK_SMALL(in[i] - inv[i], thresh); + } + + // do it again, just in case some internal state was modified in inverse + + d.forwardUnitary(in.data(), out.data()); + + for (int i = 0; i < n; ++i) { + BOOST_CHECK_SMALL(expected[i] - out[i], thresh); + } +} + +BOOST_AUTO_TEST_CASE(u5) +{ + int n = 5; + vector<double> in { 1, 2, 3, 5, 6 }; + vector<double> out(n); + vector<double> inv(n); + vector<double> expected { 7.6026, -4.1227, 0.3162, -0.0542, -0.3162 }; + + DCT d(n); + + // our expected values are very approximate + double thresh = 1e-4; + + d.forwardUnitary(in.data(), out.data()); + + for (int i = 0; i < n; ++i) { + BOOST_CHECK_SMALL(expected[i] - out[i], thresh); + } + + d.inverseUnitary(out.data(), inv.data()); + + for (int i = 0; i < n; ++i) { + BOOST_CHECK_SMALL(in[i] - inv[i], thresh); + } + + // do it again, just in case some internal state was modified in inverse + + d.forwardUnitary(in.data(), out.data()); + + for (int i = 0; i < n; ++i) { + BOOST_CHECK_SMALL(expected[i] - out[i], thresh); + } +} + +BOOST_AUTO_TEST_SUITE_END() +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/TestFilter.cpp Wed Oct 07 11:15:13 2015 +0100 @@ -0,0 +1,96 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ + +#include "dsp/signalconditioning/Filter.h" + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MAIN + +#include <boost/test/unit_test.hpp> + +#include <stdexcept> + +using namespace std; + +BOOST_AUTO_TEST_SUITE(TestFilter) + +static vector<double> in { 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16 }; + +static vector<double> iir_a { 1,5.75501989315662,16.326056867468,28.779190797823,34.2874379215653,28.137815126537,15.6064643257793,5.37874515231553,0.913800050254382,0.0,0.0 }; +static vector<double> iir_b { 0.0031954608137085,0.0180937089815597,0.0508407778575426,0.0895040074158415,0.107385387168148,0.0895040074158415,0.0508407778575426,0.0180937089815597,0.0031954608137085,0.0,0.0 }; +static vector<double> iir_expected { 0.003195460813709, 0.006094690058282, 0.009370240771381, 0.012857578361690, 0.015328760300750, 0.019107809614909, 0.022257958968869, 0.024598034053011, 0.029106103380941, 0.031152166476509, 0.034424013713795, 0.038775350541015, 0.039924063374886, 0.044846280036012, 0.047614917256999, 0.049338485830505 }; + +static vector<double> fir_b { -1.5511e-18,-0.022664,1.047e-17,0.27398,0.49737,0.27398,1.047e-17,-0.022664,-1.5511e-18 }; +static vector<double> fir_expected { -1.5511e-18,-0.022664,-0.045328,0.20599,0.95467,1.9773,3,4,5,6,7,8,9,10,11,12 }; + +BOOST_AUTO_TEST_CASE(iir) +{ + vector<double> a(iir_a); + vector<double> b(iir_b); + vector<double> expected(iir_expected); + + Filter f({ a, b }); + + int n = expected.size(); + vector<double> out(n, 0.0); + + f.process(in.data(), out.data(), n); + + double thresh = 1e-12; + + for (int i = 0; i < n; ++i) { + BOOST_CHECK_SMALL(out[i] - expected[i], thresh); + } +} + +BOOST_AUTO_TEST_CASE(iir_chunked) +{ + vector<double> a(iir_a); + vector<double> b(iir_b); + vector<double> expected(iir_expected); + + Filter f({ a, b }); + + int n = expected.size(); + vector<double> out(n, 0.0); + + int j = 0; + int i = 0; + while (j < n) { + if (++i == 4) { + i = 1; + } + if (j + i >= n) { + i = n - j; + } + f.process(in.data() + j, out.data() + j, i); + j += i; + } + + double thresh = 1e-12; + + for (int i = 0; i < n; ++i) { + BOOST_CHECK_SMALL(out[i] - expected[i], thresh); + } +} + +BOOST_AUTO_TEST_CASE(fir) +{ + vector<double> b(fir_b); + vector<double> expected(fir_expected); + + Filter f({ {}, b }); + + int n = expected.size(); + vector<double> out(n, 0.0); + + f.process(in.data(), out.data(), n); + + double thresh = 1e-4; + + for (int i = 0; i < n; ++i) { + BOOST_CHECK_SMALL(out[i] - expected[i], thresh); + } +} + +BOOST_AUTO_TEST_SUITE_END() +
--- a/tests/TestMathUtilities.cpp Wed Oct 07 09:33:28 2015 +0100 +++ b/tests/TestMathUtilities.cpp Wed Oct 07 11:15:13 2015 +0100 @@ -3,12 +3,15 @@ #include "maths/MathUtilities.h" #include <cmath> +#include <iostream> #define BOOST_TEST_DYN_LINK #define BOOST_TEST_MAIN #include <boost/test/unit_test.hpp> +using namespace std; + BOOST_AUTO_TEST_SUITE(TestMathUtilities) BOOST_AUTO_TEST_CASE(round) @@ -34,7 +37,7 @@ BOOST_CHECK_EQUAL(MathUtilities::mean(d0, 4), 1.5); double d1[] = { -2.6 }; BOOST_CHECK_EQUAL(MathUtilities::mean(d1, 1), -2.6); - std::vector<double> v; + vector<double> v; v.push_back(0); v.push_back(4); v.push_back(3); @@ -148,6 +151,98 @@ BOOST_CHECK_EQUAL(MathUtilities::gcd(37, 18), 1); } +BOOST_AUTO_TEST_CASE(getAlphaNorm1) +{ + vector<double> in { -1, 1.5, 3, 5 }; + double out = MathUtilities::getAlphaNorm(in, 1); + double expected = 2.6250; + double thresh = 1e-4; + BOOST_CHECK_SMALL(out - expected, thresh); +} + +BOOST_AUTO_TEST_CASE(getAlphaNorm2) +{ + vector<double> in { -1, 1.5, 3, 5 }; + double out = MathUtilities::getAlphaNorm(in, 2); + double expected = 3.0516; + double thresh = 1e-4; + BOOST_CHECK_SMALL(out - expected, thresh); +} + +BOOST_AUTO_TEST_CASE(getL1Norm) +{ + vector<double> in { -1, 1.5, 3, 5 }; + double out = MathUtilities::getLpNorm(in, 1); + // L1 norm is the sum of magnitudes + double expected = 10.5; + double thresh = 1e-5; + BOOST_CHECK_SMALL(out - expected, thresh); +} + +BOOST_AUTO_TEST_CASE(getL2Norm) +{ + vector<double> in { -1, 1.5, 3, 5 }; + double out = MathUtilities::getLpNorm(in, 2); + // L2 norm is the sqrt of sum of squared magnitudes + double expected = sqrt(37.25); + double thresh = 1e-5; + BOOST_CHECK_SMALL(out - expected, thresh); +} + +BOOST_AUTO_TEST_CASE(getL3Norm) +{ + vector<double> in { -1, 1.5, 3, 5 }; + double out = MathUtilities::getLpNorm(in, 3); + double expected = 5.3875; + double thresh = 1e-4; + BOOST_CHECK_SMALL(out - expected, thresh); +} + +BOOST_AUTO_TEST_CASE(normaliseL1) +{ + vector<double> in { -1, 1.5, 3, 5 }; + vector<double> expected { -0.095238, 0.142857, 0.285714, 0.476190 }; + vector<double> out = MathUtilities::normaliseLp(in, 1); + double thresh = 1e-5; + for (int i = 0; i < int(out.size()); ++i) { + BOOST_CHECK_SMALL(out[i] - expected[i], thresh); + } + out = MathUtilities::normaliseLp({ 0, 0, 0, 0 }, 1); + for (int i = 0; i < int(out.size()); ++i) { + BOOST_CHECK_EQUAL(out[i], 0.25); + } +} + +BOOST_AUTO_TEST_CASE(normaliseL2) +{ + vector<double> in { -1, 1.5, 3, 5 }; + vector<double> expected { -0.16385, 0.24577, 0.49154, 0.81923 }; + vector<double> out = MathUtilities::normaliseLp(in, 2); + double thresh = 1e-5; + for (int i = 0; i < int(out.size()); ++i) { + BOOST_CHECK_SMALL(out[i] - expected[i], thresh); + } + out = MathUtilities::normaliseLp({ 0, 0, 0, 0 }, 2); + for (int i = 0; i < int(out.size()); ++i) { + BOOST_CHECK_EQUAL(out[i], 0.5); + } +} + +BOOST_AUTO_TEST_CASE(normaliseL3) +{ + vector<double> in { -1, 1.5, 3, 5 }; + vector<double> expected { -0.18561, 0.27842, 0.55684, 0.92807 }; + vector<double> out = MathUtilities::normaliseLp(in, 3); + double thresh = 1e-5; + for (int i = 0; i < int(out.size()); ++i) { + BOOST_CHECK_SMALL(out[i] - expected[i], thresh); + } + out = MathUtilities::normaliseLp({ 0, 0, 0, 0 }, 3); + for (int i = 0; i < int(out.size()); ++i) { + BOOST_CHECK_SMALL(out[i] - 0.62996, 1e-5); + } +} + BOOST_AUTO_TEST_SUITE_END()