Mercurial > hg > qm-dsp
changeset 225:49844bc8a895
* Queen Mary C++ DSP library
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/base/Pitch.cpp Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,41 @@ +/* -*- 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 Copyright 2006 Chris Cannam. + All rights reserved. +*/ + +#include "Pitch.h" + +#include <cmath> + +float +Pitch::getFrequencyForPitch(int midiPitch, + float centsOffset, + float concertA) +{ + float p = float(midiPitch) + (centsOffset / 100); + return concertA * powf(2.0, (p - 69.0) / 12.0); +} + +int +Pitch::getPitchForFrequency(float frequency, + float *centsOffsetReturn, + float concertA) +{ + float p = 12.0 * (log(frequency / (concertA / 2.0)) / log(2.0)) + 57.0; + + int midiPitch = int(p + 0.00001); + float centsOffset = (p - midiPitch) * 100.0; + + if (centsOffset >= 50.0) { + midiPitch = midiPitch + 1; + centsOffset = -(100.0 - centsOffset); + } + + if (centsOffsetReturn) *centsOffsetReturn = centsOffset; + return midiPitch; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/base/Pitch.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,26 @@ +/* -*- 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 Copyright 2006 Chris Cannam. + All rights reserved. +*/ + +#ifndef _PITCH_H_ +#define _PITCH_H_ + +class Pitch +{ +public: + static float getFrequencyForPitch(int midiPitch, + float centsOffset = 0, + float concertA = 440.0); + + static int getPitchForFrequency(float frequency, + float *centsOffsetReturn = 0, + float concertA = 440.0); +}; + + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/base/Window.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,120 @@ +/* -*- 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 Copyright 2006 Chris Cannam. + All rights reserved. +*/ + +#ifndef _WINDOW_H_ +#define _WINDOW_H_ + +#include <cmath> +#include <iostream> +#include <map> + +enum WindowType { + RectangularWindow, + BartlettWindow, + HammingWindow, + HanningWindow, + BlackmanWindow, + GaussianWindow, + ParzenWindow +}; + +template <typename T> +class Window +{ +public: + /** + * Construct a windower of the given type. + */ + Window(WindowType type, size_t size) : m_type(type), m_size(size) { encache(); } + Window(const Window &w) : m_type(w.m_type), m_size(w.m_size) { encache(); } + Window &operator=(const Window &w) { + if (&w == this) return *this; + m_type = w.m_type; + m_size = w.m_size; + encache(); + return *this; + } + virtual ~Window() { delete m_cache; } + + void cut(T *src) const { cut(src, src); } + void cut(T *src, T *dst) const { + for (size_t i = 0; i < m_size; ++i) dst[i] = src[i] * m_cache[i]; + } + + WindowType getType() const { return m_type; } + size_t getSize() const { return m_size; } + +protected: + WindowType m_type; + size_t m_size; + T *m_cache; + + void encache(); +}; + +template <typename T> +void Window<T>::encache() +{ + size_t n = m_size; + T *mult = new T[n]; + size_t i; + for (i = 0; i < n; ++i) mult[i] = 1.0; + + switch (m_type) { + + case RectangularWindow: + for (i = 0; i < n; ++i) { + mult[i] = mult[i] * 0.5; + } + break; + + case BartlettWindow: + for (i = 0; i < n/2; ++i) { + mult[i] = mult[i] * (i / T(n/2)); + mult[i + n/2] = mult[i + n/2] * (1.0 - (i / T(n/2))); + } + break; + + case HammingWindow: + for (i = 0; i < n; ++i) { + mult[i] = mult[i] * (0.54 - 0.46 * cos(2 * M_PI * i / n)); + } + break; + + case HanningWindow: + for (i = 0; i < n; ++i) { + mult[i] = mult[i] * (0.50 - 0.50 * cos(2 * M_PI * i / n)); + } + break; + + case BlackmanWindow: + for (i = 0; i < n; ++i) { + mult[i] = mult[i] * (0.42 - 0.50 * cos(2 * M_PI * i / n) + + 0.08 * cos(4 * M_PI * i / n)); + } + break; + + case GaussianWindow: + for (i = 0; i < n; ++i) { + mult[i] = mult[i] * exp((-1.0 / (n*n)) * ((T(2*i) - n) * + (T(2*i) - n))); + } + break; + + case ParzenWindow: + for (i = 0; i < n; ++i) { + mult[i] = mult[i] * (1.0 - fabs((T(2*i) - n) / T(n + 1))); + } + break; + } + + m_cache = mult; +} + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/chromagram/ChromaProcess.cpp Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,213 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + + +#include "ChromaProcess.h" +#include "dsp/maths/Histogram.h" +#include <math.h> +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +ChromaProcess::ChromaProcess() +{ + +} + +ChromaProcess::~ChromaProcess() +{ + +} + +int ChromaProcess::findChromaBias( vector<double> chromaVector, unsigned int BPO, unsigned int frames ) +{ + vector<double> newChroma; + vector<int> peakIndex; + vector<int> modPeakIndex; + + unsigned int chromaLength = chromaVector.size(); + + unsigned int newLength = chromaLength + (2*BPO); + + newChroma.resize( newLength ); + newChroma.clear(); + + modPeakIndex.resize( newLength ); + modPeakIndex.clear(); + + //adds last row at the top and first row at the bottom to create + //circularity - effectively adds 2 to the bpo-length of the vectors: + + for( unsigned int i = 0; i < BPO; i++ ) + { + newChroma.push_back( chromaVector[ chromaLength - BPO + i ] ); + } + + for( unsigned i = 0; i < chromaLength; i++ ) + { + newChroma.push_back( chromaVector[ i ] ); + } + + for( unsigned i = 0; i < BPO; i++ ) + { + newChroma.push_back( chromaVector[ i ] ); + } + + // pick peaks in the chroma + peakIndex = getPeaks( newChroma, BPO ); + + // modularises to res = bpo/12 bins: + // corrects the mod value for bin 3 + modPeakIndex = mod( peakIndex, 3 ); + + // finds the highest concentration of peaks on the bpo/12 bin resolution + THistogram<int> m_hist(3); + + double ave, adev, sdev, var, skew, ccurt; + + m_hist.compute(modPeakIndex); + + m_hist.getMoments( modPeakIndex, ave, adev, sdev, var, skew, ccurt ); + + vector <double> histogram = m_hist.geTHistogramD(); + ////////////////////////////////////////////////////////////////////////////// + + /////////////////////////////////////////////////////////////////////////// + // Find actual bias from histogram + int minIdx, maxIdx; + double min, max; + + findHistMaxMin( histogram, &max, &maxIdx, &min, &minIdx ); + +/* + FILE* foutchroma = fopen("../testdata/newchroma.bin","wb"); + FILE* foutpeaks = fopen("../testdata/peaks.bin","wb"); + + + fwrite( &chromaVector[0], sizeof(double), chromaVector.size(), foutchroma ); + fwrite( &histogram[0], sizeof(double), histogram.size(), foutpeaks ); + + fclose( foutchroma ); + fclose( foutpeaks ); +*/ + return maxIdx - 1; +} + + +vector <int> ChromaProcess::getPeaks(vector <double> chroma, unsigned int BPO) +{ + vector <int> peaks; + + double pre = 0; + double post = 0; + double current = 0; + + unsigned int BPOCounter = 0; + unsigned int mult = 0; + unsigned int idx = 0; + + for( unsigned int i = 0; i < chroma.size() - 0; i++ ) + { + BPOCounter++; + + pre = chroma[ i ]; + current = chroma[ i + 1 ]; + post = chroma[ i + 2 ]; + + if( (current > 0) && (current > pre) && (current > post) ) + { + peaks.push_back( BPOCounter + 1); + } + + + if( BPOCounter == (BPO - 2 ) ) + { + BPOCounter = 0; + i+=2; + } + + } + + /* + for( unsigned int i = 1; i < chroma.size() - 1; i++ ) + { + BPOCounter++ ; + + pre = chroma[ i - 1 ]; + current = chroma[ i ]; + post = chroma[ i + 1 ]; + + if( (current > 0) && (current > pre) && (current > post) ) + { + peaks.push_back( BPOCounter + 1 ); + } + + if( BPOCounter == (PO - 1) ) + { + BPOCounter = 1; + i+=2; + } + } + */ + return peaks; +} + +vector <int> ChromaProcess::mod(vector <int> input, int res) +{ + vector <int> result; + + for( unsigned int i = 0; i < input.size(); i++ ) + { + int val = input[ i ]; + int res = val - res * floor( (double)val / (double)res ); + + if( val != 0 ) + { + if( res == 0 ) + res = 3; + + result.push_back( res ); + } + else + { + result.push_back( val ); + } + } + return result; +} + +void ChromaProcess::findHistMaxMin( vector<double> hist, double* max, int* maxIdx, double* min, int* minIdx ) +{ + double temp = 0.0; + unsigned int vecLength = hist.size(); + + *minIdx = 0; + *maxIdx = 0; + + *min = hist[0]; + *max = *min; + + for( unsigned int u = 0; u < vecLength; u++ ) + { + temp = hist[ u ]; + + if( temp < *min ) + { + *min = temp ; + *minIdx = u; + } + if( temp > *max ) + { + *max = temp ; + *maxIdx = u; + } + + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/chromagram/ChromaProcess.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,30 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#ifndef CHROMAPROCESS_H +#define CHROMAPROCESS_H + +#include <vector> + +using namespace std; + +class ChromaProcess +{ +public: + void findHistMaxMin( vector<double> hist, double* max, int*maxIdx, double* min, int* minIdx ); + vector <int> mod( vector <int> input, int res ); + vector <int> getPeaks( vector <double> chroma, unsigned int BPO ); + int findChromaBias( vector<double> chromaVector, unsigned int BPO, unsigned int frames ); + ChromaProcess(); + virtual ~ChromaProcess(); + +}; + +#endif // !defined(CHROMAPROCESS_H)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/chromagram/Chromagram.cpp Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,146 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#include <iostream> +#include <cmath> +#include "dsp/maths/MathUtilities.h" +#include "Chromagram.h" + +//---------------------------------------------------------------------------- + +Chromagram::Chromagram( ChromaConfig Config ) +{ + initialise( Config ); +} + +int Chromagram::initialise( ChromaConfig Config ) +{ + m_FMin = Config.min; // min freq + m_FMax = Config.max; // max freq + m_BPO = Config.BPO; // bins per octave + isNormalised = Config.isNormalised; // true 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)); + + // Create array for chroma result + m_chromadata = new double[ m_BPO ]; + + // Initialise FFT object + m_FFT = new FFT; + + // Create Config Structure for ConstantQ operator + CQConfig ConstantQConfig; + + // Populate CQ config structure with parameters + // inherited from the Chroma config + ConstantQConfig.FS = Config.FS; + ConstantQConfig.min = m_FMin; + ConstantQConfig.max = m_FMax; + ConstantQConfig.BPO = m_BPO; + ConstantQConfig.CQThresh = Config.CQThresh; + + // Initialise ConstantQ operator + m_ConstantQ = new ConstantQ( ConstantQConfig ); + + // Initialise working arrays + m_frameSize = m_ConstantQ->getfftlength(); + m_hopSize = m_ConstantQ->gethop(); + + m_FFTRe = new double[ m_frameSize ]; + m_FFTIm = new double[ m_frameSize ]; + m_CQRe = new double[ m_uK ]; + m_CQIm = new double[ m_uK ]; + + // Generate CQ Kernel + m_ConstantQ->sparsekernel(); + return 1; +} + +Chromagram::~Chromagram() +{ + deInitialise(); +} + +int Chromagram::deInitialise() +{ + delete [] m_chromadata; + + delete m_FFT; + + delete m_ConstantQ; + + delete [] m_FFTRe; + delete [] m_FFTIm; + delete [] m_CQRe; + delete [] m_CQIm; + return 1; +} + +//---------------------------------------------------------------------------------- +// returns the absolute value of complex number xx + i*yy +double Chromagram::kabs(double xx, double yy) +{ + double ab = sqrt(xx*xx + yy*yy); + return(ab); +} +//----------------------------------------------------------------------------------- + + +void Chromagram::unityNormalise(double *src) +{ + double min, max; + + double val = 0; + + MathUtilities::getFrameMinMax( src, m_BPO, & min, &max ); + + for( unsigned int i = 0; i < m_BPO; i++ ) + { + val = src[ i ] / max; + + src[ i ] = val; + } +} + + +double* Chromagram::process( double *data ) +{ + //initialise chromadata to 0 + for (unsigned i=0; i<m_BPO; i++) + m_chromadata[i]=0; + + double cmax = 0.0; + double cval = 0; + + // FFT of current frame + m_FFT->process( m_frameSize, 0, data, NULL, m_FFTRe, m_FFTIm ); + + // Calculate ConstantQ frame + m_ConstantQ->process( m_FFTRe, m_FFTIm, 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++) + { + unsigned firstBin = octave*m_BPO; + for (unsigned i=0; i<m_BPO; i++) + { + m_chromadata[i] += kabs( m_CQRe[ firstBin + i ], m_CQIm[ firstBin + i ]); + } + } + + if( isNormalised ) + unityNormalise( m_chromadata ); + + return m_chromadata; +} + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/chromagram/Chromagram.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,69 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#ifndef CHROMAGRAM_H +#define CHROMAGRAM_H + +#include "dsp/transforms/FFT.h" +#include "ConstantQ.h" + +struct ChromaConfig{ + unsigned int FS; + double min; + double max; + unsigned int BPO; + double CQThresh; + bool isNormalised; +}; + +class Chromagram +{ + +public: + Chromagram( ChromaConfig Config ); + ~Chromagram(); + + double* process( double *data ); + void unityNormalise( double* src ); + + // Complex arithmetic + 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; } + +private: + int initialise( ChromaConfig Config ); + int deInitialise(); + + double* m_chromadata; + double m_FMin; + double m_FMax; + unsigned int m_BPO; + unsigned int m_uK; + + bool isNormalised; + + unsigned int m_frameSize; + unsigned int m_hopSize; + + FFT* m_FFT; + ConstantQ* m_ConstantQ; + + double* m_FFTRe; + double* m_FFTIm; + double* m_CQRe; + double* m_CQIm; + +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/chromagram/ConstantQ.cpp Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,193 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#include "ConstantQ.h" +#include "dsp/transforms/FFT.h" + +//--------------------------------------------------------------------------- +// nextpow2 returns the smallest integer n such that 2^n >= x. +static double nextpow2(double x) { + double y = ceil(log(x)/log(2.0)); + return(y); +} + +static double squaredModule(const double & xx, const double & yy) { + return xx*xx + yy*yy; +} + +//---------------------------------------------------------------------------- + +ConstantQ::ConstantQ( CQConfig Config ) +{ + initialise( Config ); +} + +ConstantQ::~ConstantQ() +{ + deInitialise(); +} + +//---------------------------------------------------------------------------- +void ConstantQ::sparsekernel() +{ + //generates spectral kernel matrix (upside down?) + // initialise temporal kernel with zeros, twice length to deal w. complex numbers + + double* hammingWindowRe = new double [ m_FFTLength ]; + double* hammingWindowIm = new double [ m_FFTLength ]; + double* transfHammingWindowRe = new double [ m_FFTLength ]; + double* transfHammingWindowIm = new double [ m_FFTLength ]; + + for (unsigned u=0; u < m_FFTLength; u++) + { + hammingWindowRe[u] = 0; + hammingWindowIm[u] = 0; + } + + + // Here, fftleng*2 is a guess of the number of sparse cells in the matrix + // The matrix K x fftlength but the non-zero cells are an antialiased + // square root function. So mostly is a line, with some grey point. + m_sparseKernelIs.reserve( m_FFTLength*2 ); + m_sparseKernelJs.reserve( m_FFTLength*2 ); + m_sparseKernelRealValues.reserve( m_FFTLength*2 ); + m_sparseKernelImagValues.reserve( m_FFTLength*2 ); + + // for each bin value K, calculate temporal kernel, take its fft to + //calculate the spectral kernel then threshold it to make it sparse and + //add it to the sparse kernels matrix + double squareThreshold = m_CQThresh * m_CQThresh; + + FFT m_FFT; + + for (unsigned k = m_uK; k--; ) + { + // Computing a hamming window + const unsigned hammingLength = (int) ceil( m_dQ * m_FS / ( m_FMin * pow(2,((double)(k))/(double)m_BPO))); + for (unsigned i=0; i<hammingLength; i++) + { + const double angle = 2*PI*m_dQ*i/hammingLength; + const double real = cos(angle); + const double imag = sin(angle); + const double absol = hamming(hammingLength, i)/hammingLength; + hammingWindowRe[ i ] = absol*real; + hammingWindowIm[ i ] = absol*imag; + } + + //do fft of hammingWindow + m_FFT.process( m_FFTLength, 0, hammingWindowRe, hammingWindowIm, transfHammingWindowRe, transfHammingWindowIm ); + + + for (unsigned j=0; j<( m_FFTLength ); j++) + { + // perform thresholding + const double squaredBin = squaredModule( transfHammingWindowRe[ j ], transfHammingWindowIm[ j ]); + if (squaredBin <= squareThreshold) continue; + + // Insert non-zero position indexes, doubled because they are floats + m_sparseKernelIs.push_back(j); + m_sparseKernelJs.push_back(k); + + // take conjugate, normalise and add to array sparkernel + m_sparseKernelRealValues.push_back( transfHammingWindowRe[ j ]/m_FFTLength); + m_sparseKernelImagValues.push_back(-transfHammingWindowIm[ j ]/m_FFTLength); + } + + } + + delete [] hammingWindowRe; + delete [] hammingWindowIm; + delete [] transfHammingWindowRe; + delete [] transfHammingWindowIm; + +} + +//----------------------------------------------------------------------------- +double* ConstantQ::process( double* fftdata ) +{ + for (unsigned row=0; row<2*m_uK; row++) + { + m_CQdata[ row ] = 0; + m_CQdata[ row+1 ] = 0; + } + const unsigned *fftbin = &(m_sparseKernelIs[0]); + const unsigned *cqbin = &(m_sparseKernelJs[0]); + const double *real = &(m_sparseKernelRealValues[0]); + const double *imag = &(m_sparseKernelImagValues[0]); + const unsigned int sparseCells = m_sparseKernelRealValues.size(); + + for (unsigned i = 0; i<sparseCells; i++) + { + const unsigned row = cqbin[i]; + const unsigned col = fftbin[i]; + const double & r1 = real[i]; + const double & i1 = imag[i]; + const double & r2 = fftdata[ (2*m_FFTLength) - 2*col]; + const double & i2 = fftdata[ (2*m_FFTLength) - 2*col+1]; + // add the multiplication + m_CQdata[ 2*row ] += (r1*r2 - i1*i2); + m_CQdata[ 2*row+1] += (r1*i2 + i1*r2); + } + + return m_CQdata; +} + + +void ConstantQ::initialise( CQConfig Config ) +{ + m_FS = Config.FS; + m_FMin = Config.min; // min freq + m_FMax = Config.max; // max freq + m_BPO = Config.BPO; // bins per octave + m_CQThresh = Config.CQThresh;// ConstantQ threshold for kernel generation + + m_dQ = 1/(pow(2,(1/(double)m_BPO))-1); // Work out Q value for Filter bank + m_uK = (unsigned int) ceil(m_BPO * log(m_FMax/m_FMin)/log(2.0)); // No. of constant Q bins + + // work out length of fft required for this constant Q Filter bank + m_FFTLength = (int) pow(2, nextpow2(ceil( m_dQ*m_FS/m_FMin ))); + + m_hop = m_FFTLength/8; // <------ hop size is window length divided by 32 + + // allocate memory for cqdata + m_CQdata = new double [2*m_uK]; +} + +void ConstantQ::deInitialise() +{ + delete [] m_CQdata; +} + +void ConstantQ::process(double *FFTRe, double* FFTIm, double *CQRe, double *CQIm) +{ + for (unsigned row=0; row<m_uK; row++) + { + CQRe[ row ] = 0; + CQIm[ row ] = 0; + } + + const unsigned *fftbin = &(m_sparseKernelIs[0]); + const unsigned *cqbin = &(m_sparseKernelJs[0]); + const double *real = &(m_sparseKernelRealValues[0]); + const double *imag = &(m_sparseKernelImagValues[0]); + const unsigned int sparseCells = m_sparseKernelRealValues.size(); + + for (unsigned i = 0; i<sparseCells; i++) + { + const unsigned row = cqbin[i]; + const unsigned col = fftbin[i]; + const double & r1 = real[i]; + const double & i1 = imag[i]; + const double & r2 = FFTRe[ m_FFTLength- col]; + const double & i2 = FFTIm[ m_FFTLength - col]; + // add the multiplication + CQRe[ row ] += (r1*r2 - i1*i2); + CQIm[ row ] += (r1*i2 + i1*r2); + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/chromagram/ConstantQ.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,73 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#ifndef CONSTANTQ_H +#define CONSTANTQ_H + +#include <vector> +#include "dsp/maths/MathAliases.h" +#include "dsp/maths/MathUtilities.h" + +struct CQConfig{ + unsigned int FS; + double min; + double max; + unsigned int BPO; + double CQThresh; +}; + +class ConstantQ { + +//public functions incl. sparsekernel so can keep out of loop in main +public: + void process( double* FFTRe, double* FFTIm, double* CQRe, double* CQIm ); + + ConstantQ( CQConfig Config ); + ~ConstantQ(); + + double* process( double* FFTData ); + + void sparsekernel(); + + double hamming(int len, int n) { + double out = 0.54 - 0.46*cos(2*PI*n/len); + return(out); + } + + int getnumwin() { return m_numWin;} + double getQ() { return m_dQ;} + int getK() {return m_uK ;} + int getfftlength() { return m_FFTLength;} + int gethop() { return m_hop;} + +private: + void initialise( CQConfig Config ); + void deInitialise(); + + double* m_CQdata; + unsigned int m_FS; + double m_FMin; + double m_FMax; + double m_dQ; + double m_CQThresh; + unsigned int m_numWin; + unsigned int m_hop; + unsigned int m_BPO; + unsigned int m_FFTLength; + unsigned int m_uK; + std::vector<unsigned> m_sparseKernelIs; + std::vector<unsigned> m_sparseKernelJs; + std::vector<double> m_sparseKernelImagValues; + std::vector<double> m_sparseKernelRealValues; +}; + + +#endif//CONSTANTQ_H +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/maths/Correlation.cpp Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,51 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#include "Correlation.h" + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +Correlation::Correlation() +{ + +} + +Correlation::~Correlation() +{ + +} + +void Correlation::doAutoUnBiased(double *src, double *dst, unsigned int length) +{ + double tmp = 0.0; + double outVal = 0.0; + + unsigned int i,j; + + for( i = 0; i < length; i++) + { + for( j = i; j < length; j++) + { + tmp += src[ j-i ] * src[ j ]; + } + + + outVal = tmp / ( length - i ); + + if( outVal <= 0 ) + dst[ i ] = EPS; + else + dst[ i ] = outVal; + + tmp = 0.0; + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/maths/Correlation.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,25 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#ifndef CORRELATION_H +#define CORRELATION_H + +#define EPS 2.2204e-016 + +class Correlation +{ +public: + void doAutoUnBiased( double* src, double* dst, unsigned int length ); + Correlation(); + virtual ~Correlation(); + +}; + +#endif //
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/maths/Histogram.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,390 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ + +// Histogram.h: interface for the THistogram class. +// +////////////////////////////////////////////////////////////////////// + + +#ifndef HISTOGRAM_H +#define HISTOGRAM_H + + +#include <valarray> + +/*! \brief A histogram class + + This class computes the histogram of a vector. + +\par Template parameters + + - T type of input data (can be any: float, double, int, UINT, etc...) + - TOut type of output data: float or double. (default is double) + +\par Moments: + + The moments (average, standard deviation, skewness, etc.) are computed using +the algorithm of the Numerical recipies (see Numerical recipies in C, Chapter 14.1, pg 613). + +\par Example: + + This example shows the typical use of the class: +\code + // a vector containing the data + vector<float> data; + // Creating histogram using float data and with 101 containers, + THistogram<float> histo(101); + // computing the histogram + histo.compute(data); +\endcode + +Once this is done, you can get a vector with the histogram or the normalized histogram (such that it's area is 1): +\code + // getting normalized histogram + vector<float> v=histo.getNormalizedHistogram(); +\endcode + +\par Reference + + Equally spaced acsissa function integration (used in #GetArea): Numerical Recipies in C, Chapter 4.1, pg 130. + +\author Jonathan de Halleux, dehalleux@auto.ucl.ac.be, 2002 +*/ + +template<class T, class TOut = double> +class THistogram +{ +public: + //! \name Constructors + //@{ + /*! Default constructor + \param counters the number of histogram containers (default value is 10) + */ + THistogram(unsigned int counters = 10); + virtual ~THistogram() { clear();}; + //@} + + //! \name Histogram computation, update + //@{ + /*! Computes histogram of vector v + \param v vector to compute histogram + \param computeMinMax set to true if min/max of v have to be used to get the histogram limits + + This function computes the histogram of v and stores it internally. + \sa Update, GeTHistogram + */ + void compute( const std::vector<T>& v, bool computeMinMax = true); + //! Update histogram with the vector v + void update( const std::vector<T>& v); + //! Update histogram with t + void update( const T& t); + //@} + + //! \name Resetting functions + //@{ + //! Resize the histogram. Warning this function clear the histogram. + void resize( unsigned int counters ); + //! Clears the histogram + void clear() { m_counters.clear();}; + //@} + + //! \name Setters + //@{ + /*! This function sets the minimum of the histogram spectrum. + The spectrum is not recomputed, use it with care + */ + void setMinSpectrum( const T& min ) { m_min = min; computeStep();}; + /*! This function sets the minimum of the histogram spectrum. + The spectrum is not recomputed, use it with care + */ + void setMaxSpectrum( const T& max ) { m_max = max; computeStep();}; + //@} + //! \name Getters + //@{ + //! return minimum of histogram spectrum + const T& getMinSpectrum() const { return m_min;}; + //! return maximum of histogram spectrum + const T& getMaxSpectrum() const { return m_max;}; + //! return step size of histogram containers + TOut getStep() const { return m_step;}; + //! return number of points in histogram + unsigned int getSum() const; + /*! \brief returns area under the histogram + + The Simpson rule is used to integrate the histogram. + */ + TOut getArea() const; + + /*! \brief Computes the moments of the histogram + + \param data dataset + \param ave mean + \f[ \bar x = \frac{1}{N} \sum_{j=1}^N x_j\f] + \param adev mean absolute deviation + \f[ adev(X) = \frac{1}{N} \sum_{j=1}^N | x_j - \bar x |\f] + \param var average deviation: + \f[ \mbox{Var}(X) = \frac{1}{N-1} \sum_{j=1}^N (x_j - \bar x)^2\f] + \param sdev standard deviation: + \f[ \sigma(X) = \sqrt{var(\bar x) }\f] + \param skew skewness + \f[ \mbox{Skew}(X) = \frac{1}{N}\sum_{j=1}^N \left[ \frac{x_j - \bar x}{\sigma}\right]^3\f] + \param kurt kurtosis + \f[ \mbox{Kurt}(X) = \left\{ \frac{1}{N}\sum_{j=1}^N \left[ \frac{x_j - \bar x}{\sigma}\right]^4 \right\} - 3\f] + + */ + static void getMoments(const std::vector<T>& data, TOut& ave, TOut& adev, TOut& sdev, TOut& var, TOut& skew, TOut& kurt); + + //! return number of containers + unsigned int getSize() const { return m_counters.size();}; + //! returns i-th counter + unsigned int operator [] (unsigned int i) const { ASSERT( i < m_counters.size() ); return m_counters[i];}; + //! return the computed histogram + const std::vector<unsigned int>& geTHistogram() const { return m_counters;}; + //! return the computed histogram, in TOuts + std::vector<TOut> geTHistogramD() const; + /*! return the normalized computed histogram + + \return the histogram such that the area is equal to 1 + */ + std::vector<TOut> getNormalizedHistogram() const; + //! returns left containers position + std::vector<TOut> getLeftContainers() const; + //! returns center containers position + std::vector<TOut> getCenterContainers() const; + //@} +protected: + //! Computes the step + void computeStep() { m_step = (TOut)(((TOut)(m_max-m_min)) / (m_counters.size()-1));}; + //! Data accumulators + std::vector<unsigned int> m_counters; + //! minimum of dataset + T m_min; + //! maximum of dataset + T m_max; + //! width of container + TOut m_step; +}; + +template<class T, class TOut> +THistogram<T,TOut>::THistogram(unsigned int counters) + : m_counters(counters,0), m_min(0), m_max(0), m_step(0) +{ + +} + +template<class T, class TOut> +void THistogram<T,TOut>::resize( unsigned int counters ) +{ + clear(); + + m_counters.resize(counters,0); + + computeStep(); +} + +template<class T, class TOut> +void THistogram<T,TOut>::compute( const std::vector<T>& v, bool computeMinMax) +{ + using namespace std; + unsigned int i; + int index; + + if (m_counters.empty()) + return; + + if (computeMinMax) + { + m_max = m_min = v[0]; + for (i=1;i<v.size();i++) + { + m_max = std::max( m_max, v[i]); + m_min = std::min( m_min, v[i]); + } + } + + computeStep(); + + for (i = 0;i < v.size() ; i++) + { + index=(int) floor( ((TOut)(v[i]-m_min))/m_step ) ; + + if (index >= m_counters.size() || index < 0) + return; + + m_counters[index]++; + } +} + +template<class T, class TOut> +void THistogram<T,TOut>::update( const std::vector<T>& v) +{ + if (m_counters.empty()) + return; + + computeStep(); + + TOut size = m_counters.size(); + + int index; + for (unsigned int i = 0;i < size ; i++) + { + index = (int)floor(((TOut)(v[i]-m_min))/m_step); + + if (index >= m_counters.size() || index < 0) + return; + + m_counters[index]++; + } +} + +template<class T, class TOut> +void THistogram<T,TOut>::update( const T& t) +{ + int index=(int) floor( ((TOut)(t-m_min))/m_step ) ; + + if (index >= m_counters.size() || index < 0) + return; + + m_counters[index]++; +}; + +template<class T, class TOut> +std::vector<TOut> THistogram<T,TOut>::geTHistogramD() const +{ + std::vector<TOut> v(m_counters.size()); + for (unsigned int i = 0;i<m_counters.size(); i++) + v[i]=(TOut)m_counters[i]; + + return v; +} + +template <class T, class TOut> +std::vector<TOut> THistogram<T,TOut>::getLeftContainers() const +{ + std::vector<TOut> x( m_counters.size()); + + for (unsigned int i = 0;i<m_counters.size(); i++) + x[i]= m_min + i*m_step; + + return x; +} + +template <class T, class TOut> +std::vector<TOut> THistogram<T,TOut>::getCenterContainers() const +{ + std::vector<TOut> x( m_counters.size()); + + for (unsigned int i = 0;i<m_counters.size(); i++) + x[i]= m_min + (i+0.5)*m_step; + + return x; +} + +template <class T, class TOut> +unsigned int THistogram<T,TOut>::getSum() const +{ + unsigned int sum = 0; + for (unsigned int i = 0;i<m_counters.size(); i++) + sum+=m_counters[i]; + + return sum; +} + +template <class T, class TOut> +TOut THistogram<T,TOut>::getArea() const +{ + const size_t n=m_counters.size(); + TOut area=0; + + if (n>6) + { + area=3.0/8.0*(m_counters[0]+m_counters[n-1]) + +7.0/6.0*(m_counters[1]+m_counters[n-2]) + +23.0/24.0*(m_counters[2]+m_counters[n-3]); + for (unsigned int i=3;i<n-3;i++) + { + area+=m_counters[i]; + } + } + else if (n>4) + { + area=5.0/12.0*(m_counters[0]+m_counters[n-1]) + +13.0/12.0*(m_counters[1]+m_counters[n-2]); + for (unsigned int i=2;i<n-2;i++) + { + area+=m_counters[i]; + } + } + else if (n>1) + { + area=1/2.0*(m_counters[0]+m_counters[n-1]); + for (unsigned int i=1;i<n-1;i++) + { + area+=m_counters[i]; + } + } + else + area=0; + + return area*m_step; +} + +template <class T, class TOut> +std::vector<TOut> THistogram<T,TOut>::getNormalizedHistogram() const +{ + std::vector<TOut> normCounters( m_counters.size()); + TOut area = (TOut)getArea(); + + for (unsigned int i = 0;i<m_counters.size(); i++) + { + normCounters[i]= (TOut)m_counters[i]/area; + } + + return normCounters; +}; + +template <class T, class TOut> +void THistogram<T,TOut>::getMoments(const std::vector<T>& data, TOut& ave, TOut& adev, TOut& sdev, TOut& var, TOut& skew, TOut& kurt) +{ + int j; + double ep=0.0,s,p; + const size_t n = data.size(); + + if (n <= 1) + // nrerror("n must be at least 2 in moment"); + return; + + s=0.0; // First pass to get the mean. + for (j=0;j<n;j++) + s += data[j]; + + ave=s/(n); + adev=var=skew=kurt=0.0; + /* Second pass to get the first (absolute), second, + third, and fourth moments of the + deviation from the mean. */ + + for (j=0;j<n;j++) + { + adev += fabs(s=data[j]-(ave)); + ep += s; + var += (p=s*s); + skew += (p *= s); + kurt += (p *= s); + } + + + adev /= n; + var=(var-ep*ep/n)/(n-1); // Corrected two-pass formula. + sdev=sqrt(var); // Put the pieces together according to the conventional definitions. + if (var) + { + skew /= (n*(var)*(sdev)); + kurt=(kurt)/(n*(var)*(var))-3.0; + } + else + //nrerror("No skew/kurtosis when variance = 0 (in moment)"); + return; +} + +#endif +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/maths/MathAliases.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,55 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#ifndef MATHALIASES_H +#define MATHALIASES_H + +#include <cmath> +#include <complex> + +using namespace std; +typedef complex<double> ComplexData; + + +#ifndef PI +#define PI (3.14159265358979232846) +#endif + +#define TWO_PI (*2.PI) + +#define EPS 2.2204e-016 + +/* aliases to math.h functions */ +#define EXP exp +#define COS cos +#define SIN sin +#define ABS fabs +#define POW powf +#define SQRT sqrtf +#define LOG10 log10f +#define LOG logf +#define FLOOR floorf +#define TRUNC truncf + +/* aliases to complex.h functions */ +/** sample = EXPC(complex) */ +#define EXPC cexpf +/** complex = CEXPC(complex) */ +#define CEXPC cexp +/** sample = ARGC(complex) */ +#define ARGC cargf +/** sample = ABSC(complex) norm */ +#define ABSC cabsf +/** sample = REAL(complex) */ +#define REAL crealf +/** sample = IMAG(complex) */ +#define IMAG cimagf + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/maths/MathUtilities.cpp Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,170 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#include "MathUtilities.h" + +#include <iostream> +#include <cmath> + + +double MathUtilities::mod(double x, double y) +{ + double a = floor( x / y ); + + double b = x - ( y * a ); + return b; +} + +double MathUtilities::princarg(double ang) +{ + double ValOut; + + ValOut = mod( ang + M_PI, -2 * M_PI ) + M_PI; + + return ValOut; +} + +void MathUtilities::getAlphaNorm(const double *data, unsigned int len, unsigned int alpha, double* ANorm) +{ + unsigned int i; + double temp = 0.0; + double a=0.0; + + for( i = 0; i < len; i++) + { + temp = data[ i ]; + + a += ::pow( fabs(temp), alpha ); + } + a /= ( double )len; + a = ::pow( a, ( 1.0 / (double) alpha ) ); + + *ANorm = a; +} + +double MathUtilities::getAlphaNorm( const std::vector <double> &data, unsigned int alpha ) +{ + unsigned int i; + unsigned 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), alpha ); + } + a /= ( double )len; + a = ::pow( a, ( 1.0 / (double) alpha ) ); + + return a; +} + +double MathUtilities::round(double x) +{ + double val = (double)floor(x + 0.5); + + return val; +} + +double MathUtilities::median(const double *src, unsigned int len) +{ + unsigned int i, j; + double tmp = 0.0; + double tempMedian; + double medianVal; + + double* scratch = new double[ len ];//Vector < double > sortedX = Vector < double > ( size ); + + for ( i = 0; i < len; i++ ) + { + scratch[i] = src[i]; + } + + for ( i = 0; i < len - 1; i++ ) + { + for ( j = 0; j < len - 1 - i; j++ ) + { + if ( scratch[j + 1] < scratch[j] ) + { + // compare the two neighbors + tmp = scratch[j]; // swap a[j] and a[j+1] + scratch[j] = scratch[j + 1]; + scratch[j + 1] = tmp; + } + } + } + int middle; + if ( len % 2 == 0 ) + { + middle = len / 2; + tempMedian = ( scratch[middle] + scratch[middle - 1] ) / 2; + } + else + { + middle = ( int )floor( len / 2.0 ); + tempMedian = scratch[middle]; + } + + medianVal = tempMedian; + + delete [] scratch; + return medianVal; +} + +double MathUtilities::sum(const double *src, unsigned int len) +{ + unsigned int i ; + double retVal =0.0; + + for( i = 0; i < len; i++) + { + retVal += src[ i ]; + } + + return retVal; +} + +double MathUtilities::mean(const double *src, unsigned int len) +{ + double retVal =0.0; + + double s = sum( src, len ); + + retVal = s / (double)len; + + return retVal; +} + +void MathUtilities::getFrameMinMax(const double *data, unsigned int len, double *min, double *max) +{ + unsigned int i; + double temp = 0.0; + double a=0.0; + + *min = data[0]; + *max = data[0]; + + for( i = 0; i < len; i++) + { + temp = data[ i ]; + + if( temp < *min ) + { + *min = temp ; + } + if( temp > *max ) + { + *max = temp ; + } + + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/maths/MathUtilities.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,30 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#ifndef MATHUTILITIES_H +#define MATHUTILITIES_H + +#include <vector> + +class MathUtilities +{ +public: + static double round( double x ); + static void getFrameMinMax( const double* data, unsigned int len, double* min, double* max ); + static double mean( const double* src, unsigned int len ); + static double sum( const double* src, unsigned int len ); + static double princarg( double ang ); + static double median( const double* src, unsigned int len ); + 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 <double> &data, unsigned int alpha ); +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/maths/Polyfit.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,398 @@ +/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ +//--------------------------------------------------------------------------- + +#ifndef PolyfitHPP +#define PolyfitHPP +//--------------------------------------------------------------------------- +// Least-squares curve fitting class for arbitrary data types +/* + +{ ****************************************** + **** Scientific Subroutine Library **** + **** for C++ Builder **** + ****************************************** + + The following programs were written by Allen Miller and appear in the + book entitled "Pascal Programs For Scientists And Engineers" which is + published by Sybex, 1981. + They were originally typed and submitted to MTPUG in Oct. 1982 + Juergen Loewner + Hoher Heckenweg 3 + D-4400 Muenster + They have had minor corrections and adaptations for Turbo Pascal by + Jeff Weiss + 1572 Peacock Ave. + Sunnyvale, CA 94087. + + +2000 Oct 28 Updated for Delphi 4, open array parameters. + This allows the routine to be generalised so that it is no longer + hard-coded to make a specific order of best fit or work with a + specific number of points. +2001 Jan 07 Update Web site address + +Copyright © David J Taylor, Edinburgh and others listed above +Web site: www.satsignal.net +E-mail: davidtaylor@writeme.com +}*/ + + /////////////////////////////////////////////////////////////////////////////// + // Modified by CLandone for VC6 Aug 2004 + /////////////////////////////////////////////////////////////////////////////// + + +using std::vector; + + +class TPolyFit +{ + typedef vector<vector<double> > Matrix; +public: + + static double PolyFit2 (const vector<double> &x, // does the work + const vector<double> &y, + vector<double> &coef); + + +private: + TPolyFit &operator = (const TPolyFit &); // disable assignment + TPolyFit(); // and instantiation + TPolyFit(const TPolyFit&); // and copying + + + static void Square (const Matrix &x, // Matrix multiplication routine + const vector<double> &y, + Matrix &a, // A = transpose X times X + vector<double> &g, // G = Y times X + const int nrow, const int ncol); + // Forms square coefficient matrix + + static bool GaussJordan (Matrix &b, // square matrix of coefficients + const vector<double> &y, // constant vector + vector<double> &coef); // solution vector + // returns false if matrix singular + + static bool GaussJordan2(Matrix &b, + const vector<double> &y, + Matrix &w, + vector<vector<int> > &index); +}; + +// some utility functions + +namespace NSUtility +{ + inline void swap(double &a, double &b) {double t = a; a = b; b = t;} + void zeroise(vector<double> &array, int n); + void zeroise(vector<int> &array, int n); + void zeroise(vector<vector<double> > &matrix, int m, int n); + void zeroise(vector<vector<int> > &matrix, int m, int n); + inline double sqr(const double &x) {return x * x;} +}; + +//--------------------------------------------------------------------------- +// Implementation +//--------------------------------------------------------------------------- +using namespace NSUtility; +//------------------------------------------------------------------------------------------ + + +// main PolyFit routine + +double TPolyFit::PolyFit2 (const vector<double> &x, + const vector<double> &y, + vector<double> &coefs) +// nterms = coefs.size() +// npoints = x.size() +{ + int i, j; + double xi, yi, yc, srs, sum_y, sum_y2; + Matrix xmatr; // Data matrix + Matrix a; + vector<double> g; // Constant vector + const int npoints(x.size()); + const int nterms(coefs.size()); + double correl_coef; + zeroise(g, nterms); + zeroise(a, nterms, nterms); + zeroise(xmatr, npoints, nterms); + if(nterms < 1) + throw "PolyFit called with less than one term"; + if(npoints < 2) + throw "PolyFit called with less than two points"; + if(npoints != y.size()) + throw "PolyFit called with x and y of unequal size"; + for(i = 0; i < npoints; ++i) + { + // { setup x matrix } + xi = x[i]; + xmatr[i][0] = 1.0; // { first column } + for(j = 1; j < nterms; ++j) + xmatr[i][j] = xmatr [i][j - 1] * xi; + } + Square (xmatr, y, a, g, npoints, nterms); + if(!GaussJordan (a, g, coefs)) + return -1; + sum_y = 0.0; + sum_y2 = 0.0; + srs = 0.0; + for(i = 0; i < npoints; ++i) + { + yi = y[i]; + yc = 0.0; + for(j = 0; j < nterms; ++j) + yc += coefs [j] * xmatr [i][j]; + srs += 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; + // Either return 0 or the correct value of correlation coefficient + if (correl_coef != 0) + correl_coef = srs / correl_coef; + if (correl_coef >= 1) + correl_coef = 0.0; + else + correl_coef = sqrt (1.0 - correl_coef); + return correl_coef; +} + + +//------------------------------------------------------------------------ + +// Matrix multiplication routine +// A = transpose X times X +// G = Y times X + +// Form square coefficient matrix + +void TPolyFit::Square (const Matrix &x, + const vector<double> &y, + Matrix &a, + vector<double> &g, + const int nrow, + const int ncol) +{ + int i, k, l; + for(k = 0; k < ncol; ++k) + { + for(l = 0; l < k + 1; ++l) + { + a [k][l] = 0.0; + for(i = 0; i < nrow; ++i) + { + a[k][l] += x[i][l] * x [i][k]; + if(k != l) + a[l][k] = a[k][l]; + } + } + g[k] = 0.0; + for(i = 0; i < nrow; ++i) + g[k] += y[i] * x[i][k]; + } +} +//--------------------------------------------------------------------------------- + + +bool TPolyFit::GaussJordan (Matrix &b, + const vector<double> &y, + vector<double> &coef) +//b square matrix of coefficients +//y constant vector +//coef solution vector +//ncol order of matrix got from b.size() + + +{ +/* + { Gauss Jordan matrix inversion and solution } + + { B (n, n) coefficient matrix becomes inverse } + { Y (n) original constant vector } + { W (n, m) constant vector(s) become solution vector } + { DETERM is the determinant } + { ERROR = 1 if singular } + { INDEX (n, 3) } + { NV is number of constant vectors } +*/ + + int ncol(b.size()); + int irow, icol; + vector<vector<int> >index; + Matrix w; + + zeroise(w, ncol, ncol); + zeroise(index, ncol, 3); + + if(!GaussJordan2(b, y, w, index)) + return false; + + // Interchange columns + int m; + for (int i = 0; i < ncol; ++i) + { + m = ncol - i - 1; + if(index [m][0] != index [m][1]) + { + irow = index [m][0]; + icol = index [m][1]; + for(int k = 0; k < ncol; ++k) + swap (b[k][irow], b[k][icol]); + } + } + + for(int k = 0; k < ncol; ++k) + { + if(index [k][2] != 0) + { + throw "Error in GaussJordan: matrix is singular"; + } + } + + for( int i = 0; i < ncol; ++i) + coef[i] = w[i][0]; + + + return true; +} // end; { procedure GaussJordan } +//---------------------------------------------------------------------------------------------- + + +bool TPolyFit::GaussJordan2(Matrix &b, + const vector<double> &y, + Matrix &w, + vector<vector<int> > &index) +{ + //GaussJordan2; // first half of GaussJordan + // actual start of gaussj + + double big, t; + double pivot; + double determ; + int irow, icol; + int ncol(b.size()); + int nv = 1; // single constant vector + for(int i = 0; i < ncol; ++i) + { + w[i][0] = y[i]; // copy constant vector + index[i][2] = -1; + } + determ = 1.0; + for(int i = 0; i < ncol; ++i) + { + // Search for largest element + big = 0.0; + for(int j = 0; j < ncol; ++j) + { + if(index[j][2] != 0) + { + for(int k = 0; k < ncol; ++k) + { + if(index[k][2] > 0) + throw "Error in GaussJordan2: matrix is singular"; + + if(index[k][2] < 0 && fabs(b[j][k]) > big) + { + irow = j; + icol = k; + big = fabs(b[j][k]); + } + } // { k-loop } + } + } // { j-loop } + index [icol][2] = index [icol][2] + 1; + index [i][0] = irow; + index [i][1] = icol; + + // Interchange rows to put pivot on diagonal + // GJ3 + if(irow != icol) + { + determ = -determ; + for(int m = 0; m < ncol; ++m) + swap (b [irow][m], b[icol][m]); + if (nv > 0) + for (int m = 0; m < nv; ++m) + swap (w[irow][m], w[icol][m]); + } // end GJ3 + + // divide pivot row by pivot column + pivot = b[icol][icol]; + determ *= pivot; + b[icol][icol] = 1.0; + + for(int m = 0; m < ncol; ++m) + b[icol][m] /= pivot; + if(nv > 0) + for(int m = 0; m < nv; ++m) + w[icol][m] /= pivot; + + // Reduce nonpivot rows + for(int n = 0; n < ncol; ++n) + { + if(n != icol) + { + t = b[n][icol]; + b[n][icol] = 0.0; + for(int m = 0; m < ncol; ++m) + b[n][m] -= b[icol][m] * t; + if(nv > 0) + for(int m = 0; m < nv; ++m) + w[n][m] -= w[icol][m] * t; + } + } + } // { i-loop } + return true; +} +//---------------------------------------------------------------------------------------------- + +//------------------------------------------------------------------------------------ + +// Utility functions +//-------------------------------------------------------------------------- + +// fills a vector with zeros. +void NSUtility::zeroise(vector<double> &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<int> &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<vector<double> > &matrix, int m, int n) +{ + vector<double> 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<vector<int> > &matrix, int m, int n) +{ + vector<int> zero; + zeroise(zero, n); + matrix.clear(); + for(int j = 0; j < m; ++j) + matrix.push_back(zero); +} +//-------------------------------------------------------------------------- + + +#endif +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/onsets/DetectionFunction.cpp Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,204 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#include "DetectionFunction.h" + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +DetectionFunction::DetectionFunction( DFConfig Config ) : + m_window(0) +{ + magHistory = NULL; + phaseHistory = NULL; + phaseHistoryOld = NULL; + j = ComplexData( 0, 1 ); + + initialise( Config ); +} + +DetectionFunction::~DetectionFunction() +{ + deInitialise(); +} + + +void DetectionFunction::initialise( DFConfig Config ) +{ + m_dataLength = Config.frameLength; + m_halfLength = m_dataLength/2; + m_DFType = Config.DFType; + + magHistory = new double[ m_halfLength ]; + memset(magHistory,0, m_halfLength*sizeof(double)); + + phaseHistory = new double[ m_halfLength ]; + memset(phaseHistory,0, m_halfLength*sizeof(double)); + + phaseHistoryOld = new double[ m_halfLength ]; + memset(phaseHistoryOld,0, m_halfLength*sizeof(double)); + + m_phaseVoc = new PhaseVocoder; + + m_DFWindowedFrame = new double[ m_dataLength ]; + m_magnitude = new double[ m_halfLength ]; + m_thetaAngle = new double[ m_halfLength ]; + + m_window = new Window<double>(HanningWindow, m_dataLength); +} + +void DetectionFunction::deInitialise() +{ + delete [] magHistory ; + delete [] phaseHistory ; + delete [] phaseHistoryOld ; + + delete m_phaseVoc; + + delete [] m_DFWindowedFrame; + delete [] m_magnitude; + delete [] m_thetaAngle; + + delete m_window; +} + +double DetectionFunction::process( double *TDomain ) +{ + double retVal = 0; + + m_window->cut( TDomain, m_DFWindowedFrame ); + + m_phaseVoc->process( m_dataLength, m_DFWindowedFrame, m_magnitude, m_thetaAngle ); + + switch( m_DFType ) + { + case DF_HFC: + retVal = HFC( m_halfLength, m_magnitude); + break; + + case DF_SPECDIFF: + retVal = specDiff( m_halfLength, m_magnitude); + break; + + case DF_PHASEDEV: + retVal = phaseDev( m_halfLength, m_magnitude, m_thetaAngle); + break; + + case DF_COMPLEXSD: + retVal = complexSD( m_halfLength, m_magnitude, m_thetaAngle); + break; + } + + return retVal; +} + +double DetectionFunction::HFC(unsigned int length, double *src) +{ + unsigned int i; + double val = 0; + + for( i = 0; i < length; i++) + { + val += src[ i ] * ( i + 1); + } + return val; +} + +double DetectionFunction::specDiff(unsigned int length, double *src) +{ + unsigned int i; + double val = 0.0; + double temp = 0.0; + double diff = 0.0; + + for( i = 0; i < length; i++) + { + temp = fabs( (src[ i ] * src[ i ]) - (magHistory[ i ] * magHistory[ i ]) ); + + diff= sqrt(temp); + + if( src[ i ] > 0.1) + { + val += diff; + } + + magHistory[ i ] = src[ i ]; + } + + return val; +} + + +double DetectionFunction::phaseDev(unsigned int length, double *srcMagnitude, double *srcPhase) +{ + unsigned int i; + double tmpPhase = 0; + double tmpVal = 0; + double val = 0; + + double dev = 0; + + for( i = 0; i < length; i++) + { + tmpPhase = (srcPhase[ i ]- 2*phaseHistory[ i ]+phaseHistoryOld[ i ]); + dev = MathUtilities::princarg( tmpPhase ); + + if( srcMagnitude[ i ] > 0.1) + { + tmpVal = fabs( dev); + val += tmpVal ; + } + + phaseHistoryOld[ i ] = phaseHistory[ i ] ; + phaseHistory[ i ] = srcPhase[ i ]; + } + + + return val; +} + + +double DetectionFunction::complexSD(unsigned int length, double *srcMagnitude, double *srcPhase) +{ + unsigned int i; + double val = 0; + double tmpPhase = 0; + double tmpReal = 0; + double tmpImag = 0; + + double dev = 0; + ComplexData meas = ComplexData( 0, 0 ); + + for( i = 0; i < length; i++) + { + tmpPhase = (srcPhase[ i ]- 2*phaseHistory[ i ]+phaseHistoryOld[ i ]); + dev= MathUtilities::princarg( tmpPhase ); + + meas = magHistory[i] - ( srcMagnitude[ i ] * exp( j * dev) ); + + tmpReal = real( meas ); + tmpImag = imag( meas ); + + val += sqrt( (tmpReal * tmpReal) + (tmpImag * tmpImag) ); + + phaseHistoryOld[ i ] = phaseHistory[ i ] ; + phaseHistory[ i ] = srcPhase[ i ]; + magHistory[ i ] = srcMagnitude[ i ]; + } + + return val; +} + +double* DetectionFunction::getSpectrumMagnitude() +{ + return m_magnitude; +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/onsets/DetectionFunction.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,72 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#ifndef DETECTIONFUNCTION_H +#define DETECTIONFUNCTION_H + +#include "dsp/maths/MathUtilities.h" +#include "dsp/maths/MathAliases.h" +#include "dsp/phasevocoder/PhaseVocoder.h" +#include "base/Window.h" + +#define DF_HFC (1) +#define DF_SPECDIFF (2) +#define DF_PHASEDEV (3) +#define DF_COMPLEXSD (4) + +struct DFConfig{ + double stepSecs; // DF step in seconds + unsigned int stepSize; // DF step in samples + unsigned int frameLength; // DF analysis window - usually 2*step + int DFType; // type of detection function ( see defines ) +}; + +class DetectionFunction +{ +public: + double* getSpectrumMagnitude(); + DetectionFunction( DFConfig Config ); + virtual ~DetectionFunction(); + double process( double* TDomain ); + +private: + double HFC( unsigned int length, double* src); + double specDiff( unsigned int length, double* src); + double phaseDev(unsigned int length, double *srcMagnitude, double *srcPhase); + double complexSD(unsigned int length, double *srcMagnitude, double *srcPhase); + +private: + void initialise( DFConfig Config ); + void deInitialise(); + + int m_DFType; + unsigned int m_dataLength; + unsigned int m_halfLength; + + double* magHistory; + double* phaseHistory; + double* phaseHistoryOld; + + double* m_DFWindowedFrame; // Array for windowed analysis frame + double* m_magnitude; // Magnitude of analysis frame ( frequency domain ) + double* m_thetaAngle;// Phase of analysis frame ( frequency domain ) + + + vector < ComplexData > meas ; + + ComplexData j; + + Window<double> *m_window; + + PhaseVocoder* m_phaseVoc; // Phase Vocoder + +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/onsets/PeakPicking.cpp Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,137 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#include "PeakPicking.h" +#include "dsp/maths/Polyfit.h" + + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +PeakPicking::PeakPicking( PPickParams Config ) +{ + m_workBuffer = NULL; + initialise( Config ); +} + +PeakPicking::~PeakPicking() +{ + deInitialise(); +} + +void PeakPicking::initialise( PPickParams Config ) +{ + m_DFLength = Config.length ; + Qfilta = Config.QuadThresh.a ; + Qfiltb = Config.QuadThresh.b ; + Qfiltc = Config.QuadThresh.c ; + + m_DFProcessingParams.length = m_DFLength; + m_DFProcessingParams.LPOrd = Config.LPOrd; + m_DFProcessingParams.LPACoeffs = Config.LPACoeffs; + m_DFProcessingParams.LPBCoeffs = Config.LPBCoeffs; + m_DFProcessingParams.winPre = Config.WinT.pre; + m_DFProcessingParams.winPost = Config.WinT.post; + m_DFProcessingParams.AlphaNormParam = Config.alpha; + m_DFProcessingParams.isMedianPositive = false; + + m_DFSmoothing = new DFProcess( m_DFProcessingParams ); + + m_workBuffer = new double[ m_DFLength ]; + memset( m_workBuffer, 0, sizeof(double)*m_DFLength); +} + +void PeakPicking::deInitialise() +{ + delete [] m_workBuffer; + delete m_DFSmoothing; + m_workBuffer = NULL; +} + +void PeakPicking::process( double* src, unsigned int len, vector<int> &onsets ) +{ + vector <double> m_maxima; + + // Signal conditioning + m_DFSmoothing->process( src, m_workBuffer ); + + for( unsigned int u = 0; u < len; u++) + { + m_maxima.push_back( m_workBuffer[ u ] ); + } + + quadEval( m_maxima, onsets ); + + for( int b = 0; b < m_maxima.size(); b++) + { + src[ b ] = m_maxima[ b ]; + } +} + +int PeakPicking::quadEval( vector<double> &src, vector<int> &idx ) +{ + unsigned int maxLength; + + vector <int> m_maxIndex; + vector <int> m_onsetPosition; + + vector <double> m_maxFit; + vector <double> m_poly; + vector <double> m_err; + + double p; + + m_poly.push_back(0); + m_poly.push_back(0); + m_poly.push_back(0); + + for( int t = -2; t < 3; t++) + { + m_err.push_back( (double)t ); + } + for( unsigned int i = 2; i < src.size() - 2; i++) + { + if( (src[ i ] > src[ i - 1 ]) && (src[ i ] > src[ i + 1 ]) && (src[ i ] > 0) ) + { + m_maxIndex.push_back( i + 1 ); + } + } + + maxLength = m_maxIndex.size(); + + double selMax = 0; + + for( unsigned int j = 0; j < maxLength ; j++) + { + for( int k = -3; k < 2; k++) + { + selMax = src[ m_maxIndex[j] + k ] ; + m_maxFit.push_back(selMax); + } + + p = TPolyFit::PolyFit2( m_err, m_maxFit, m_poly); + + double f = m_poly[0]; + double g = m_poly[1]; + double h = m_poly[2]; + + int kk = m_poly.size(); + + if( h < -Qfilta || f > Qfiltc) + { + idx.push_back( m_maxIndex[j] ); + } + + m_maxFit.erase( m_maxFit.begin(), m_maxFit.end() ); + } + + return 1; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/onsets/PeakPicking.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,77 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +// PeakPicking.h: interface for the PeakPicking class. +// +////////////////////////////////////////////////////////////////////// + +#ifndef PEAKPICKING_H +#define PEAKPICKING_H + +#include "dsp/maths/MathUtilities.h" +#include "dsp/maths/MathAliases.h" +#include "dsp/signalconditioning/DFProcess.h" + + +struct WinThresh +{ + unsigned int pre; + unsigned int post; +}; + +struct QFitThresh +{ + double a; + double b; + double c; +}; + +struct PPickParams +{ + unsigned int length; //Detection FunctionLength + double tau; // time resolution of the detection function: + unsigned int alpha; //alpha-norm parameter + double cutoff;//low-pass Filter cutoff freq + unsigned 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]: + QFitThresh QuadThresh; +}; + +class PeakPicking +{ +public: + PeakPicking( PPickParams Config ); + virtual ~PeakPicking(); + + void process( double* src, unsigned int len, vector<int> &onsets ); + + +private: + void initialise( PPickParams Config ); + void deInitialise(); + int quadEval( vector<double> &src, vector<int> &idx ); + + DFProcConfig m_DFProcessingParams; + + unsigned int m_DFLength ; + double m_alphaNormParam ; + double Qfilta ; + double Qfiltb; + double Qfiltc; + + + double* m_workBuffer; + + DFProcess* m_DFSmoothing; +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/phasevocoder/PhaseVocoder.cpp Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,97 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#include "PhaseVocoder.h" +#include "dsp/transforms/FFT.h" +#include <math.h> + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +PhaseVocoder::PhaseVocoder() +{ + +} + +PhaseVocoder::~PhaseVocoder() +{ + +} + +void PhaseVocoder::FFTShift(unsigned int size, double *src) +{ + // IN-place Rotation of FFT arrays + unsigned int i; + + shiftBuffer = new double[size/2]; + + for( i = 0; i < size/2; i++) + { + shiftBuffer[ i ] = src[ i ]; + src[ i ] = src[ i + size/2]; + } + + for( i =size/2; i < size; i++) + { + src[ i ] = shiftBuffer[ i -(size/2)]; + } + + delete [] shiftBuffer; + +} + +void PhaseVocoder::process(unsigned int size, double *src, double *mag, double *theta) +{ + + // Primary Interface to Phase Vocoder + realOut = new double[ size ]; + imagOut = new double[ size ]; + + FFTShift( size, src); + + coreFFT( size, src, 0, realOut, imagOut); + + getMagnitude( size/2, mag, realOut, imagOut); + getPhase( size/2, theta, realOut, imagOut); + + delete [] realOut; + delete [] imagOut; +} + + +void PhaseVocoder::coreFFT( unsigned int NumSamples, double *RealIn, double* ImagIn, double *RealOut, double *ImagOut) +{ + // This function is taken from a standard freeware implementation defined in FFT.h + // TODO: Use FFTW + FFT::process( NumSamples,0, RealIn, ImagIn, RealOut, ImagOut ); +} + +void PhaseVocoder::getMagnitude(unsigned int size, double *mag, double *real, double *imag) +{ + unsigned int j; + + for( j = 0; j < size; j++) + { + mag[ j ] = sqrt( real[ j ] * real[ j ] + imag[ j ] * imag[ j ]); + } +} + +void PhaseVocoder::getPhase(unsigned int size, double *theta, double *real, double *imag) +{ + unsigned int k; + + // Phase Angle "matlab" style + //Watch out for quadrant mapping !!! + for( k = 0; k < size; k++) + { + theta[ k ] = atan2( -imag[ k ], real[ k ]); + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/phasevocoder/PhaseVocoder.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,35 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#ifndef PHASEVOCODER_H +#define PHASEVOCODER_H + + +class PhaseVocoder +{ +public: + PhaseVocoder(); + virtual ~PhaseVocoder(); + + void process( unsigned int size, double* src, double* mag, double* theta); + void FFTShift( unsigned int size, double* src); + +protected: + void getPhase(unsigned int size, double *theta, double *real, double *imag); + void coreFFT( unsigned int NumSamples, double *RealIn, double* ImagIn, double *RealOut, double *ImagOut); + void getMagnitude( unsigned int size, double* mag, double* real, double* imag); + + double* shiftBuffer; + double* imagOut; + double* realOut; + +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/rateconversion/Decimator.cpp Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,154 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#include "Decimator.h" + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +Decimator::Decimator( unsigned int inLength, unsigned int decFactor ) +{ + + m_inputLength = 0; + m_outputLength = 0; + m_decFactor = 1; + + initialise( inLength, decFactor ); +} + +Decimator::~Decimator() +{ + deInitialise(); +} + +void Decimator::initialise( unsigned int inLength, unsigned int decFactor) +{ + m_inputLength = inLength; + m_decFactor = decFactor; + m_outputLength = m_inputLength / m_decFactor; + + decBuffer = new double[ m_inputLength ]; + + if( m_decFactor == 4 ) + { + ////////////////////////////////////////////////// + b[ 0 ] = 0.10133306904918619; + b[ 1 ] = -0.2447523353702363; + b[ 2 ] = 0.33622528590120965; + b[ 3 ] = -0.13936581560633518; + b[ 4 ] = -0.13936581560633382; + b[ 5 ] = 0.3362252859012087; + b[ 6 ] = -0.2447523353702358; + b[ 7 ] = 0.10133306904918594; + + a[ 0 ] = 1; + a[ 1 ] = -3.9035590278139427; + a[ 2 ] = 7.5299379980621133; + a[ 3 ] = -8.6890803793177511; + a[ 4 ] = 6.4578667096099176; + a[ 5 ] = -3.0242979431223631; + a[ 6 ] = 0.83043385136748382; + a[ 7 ] = -0.094420800837809335; + ////////////////////////////////////////////////// + } + else if( m_decFactor == 2 ) + { + ////////////////////////////////////////////////// + b[ 0 ] = 0.20898944260075727; + b[ 1 ] = 0.40011234879814367; + b[ 2 ] = 0.819741973072733; + b[ 3 ] = 1.0087419911682323; + b[ 4 ] = 1.0087419911682325; + b[ 5 ] = 0.81974197307273156; + b[ 6 ] = 0.40011234879814295; + b[ 7 ] = 0.20898944260075661; + + a[ 0 ] = 1; + a[ 1 ] = 0.0077331184208358217; + a[ 2 ] = 1.9853971155964376; + a[ 3 ] = 0.19296739275341004; + a[ 4 ] = 1.2330748872852182; + a[ 5 ] = 0.18705341389316466; + a[ 6 ] = 0.23659265908013868; + a[ 7 ] = 0.032352924250533946; + } + else + { + ////////////////////////////////////////////////// + b[ 0 ] = 1; + b[ 1 ] = 0; + b[ 2 ] = 0; + b[ 3 ] = 0; + b[ 4 ] = 0; + b[ 5 ] = 0; + b[ 6 ] = 0; + b[ 7 ] = 0; + + a[ 0 ] = 1; + a[ 1 ] = 0; + a[ 2 ] = 0; + a[ 3 ] = 0; + a[ 4 ] = 0; + a[ 5 ] = 0; + a[ 6 ] = 0; + a[ 7 ] = 0; + } + + resetFilter(); +} + +void Decimator::deInitialise() +{ + delete [] decBuffer; +} + +void Decimator::resetFilter() +{ + Input = Output = 0; + + o1=o2=o3=o4=o5=o6=o7=0; +} + +void Decimator::doAntiAlias(double *src, double *dst, unsigned int length) +{ + + for( unsigned int i = 0; i < length; i++ ) + { + Input = (double)src[ i ]; + + Output = Input * b[ 0 ] + o1; + + o1 = Input * b[ 1 ] - Output * a[ 1 ] + o2; + o2 = Input * b[ 2 ] - Output * a[ 2 ] + o3; + o3 = Input * b[ 3 ] - Output * a[ 3 ] + o4; + o4 = Input * b[ 4 ] - Output * a[ 4 ] + o5; + o5 = Input * b[ 5 ] - Output * a[ 5 ] + o6; + o6 = Input * b[ 6 ] - Output * a[ 6 ] + o7; + o7 = Input * b[ 7 ] - Output * a[ 7 ] ; + + dst[ i ] = Output; + } + +} + +void Decimator::process(double *src, double *dst) +{ + if( m_decFactor != 1 ) + { + doAntiAlias( src, decBuffer, m_inputLength ); + } + unsigned idx = 0; + + for( unsigned int i = 0; i < m_outputLength; i++ ) + { + dst[ idx++ ] = decBuffer[ m_decFactor * i ]; + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/rateconversion/Decimator.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,42 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#ifndef DECIMATOR_H +#define DECIMATOR_H + +class Decimator +{ +public: + void process( double* src, double* dst ); + void doAntiAlias( double* src, double* dst, unsigned int length ); + + Decimator( unsigned int inLength, unsigned int decFactor ); + virtual ~Decimator(); + +private: + void resetFilter(); + void deInitialise(); + void initialise( unsigned int inLength, unsigned int decFactor ); + + unsigned int m_inputLength; + unsigned int m_outputLength; + unsigned int m_decFactor; + + double Input; + double Output ; + + double o1,o2,o3,o4,o5,o6,o7; + + double a[ 9 ]; + double b[ 9 ]; + + double* decBuffer; +}; + +#endif //
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/signalconditioning/DFProcess.cpp Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,172 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#include "DFProcess.h" +#include "dsp/maths/MathUtilities.h" + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +DFProcess::DFProcess( DFProcConfig Config ) +{ + filtSrc = NULL; + filtDst = NULL; + m_filtScratchIn = NULL; + m_filtScratchOut = NULL; + + m_FFOrd = 0; + + initialise( Config ); +} + +DFProcess::~DFProcess() +{ + deInitialise(); +} + +void DFProcess::initialise( DFProcConfig Config ) +{ + m_length = Config.length; + m_winPre = Config.winPre; + m_winPost = Config.winPost; + m_alphaNormParam = Config.AlphaNormParam; + + m_isMedianPositive = Config.isMedianPositive; + + 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 ); +} + +void DFProcess::deInitialise() +{ + delete [] filtSrc; + + delete [] filtDst; + + delete [] m_filtScratchIn; + + delete [] m_filtScratchOut; + + delete m_FiltFilt; +} + +void DFProcess::process(double *src, double* dst) +{ + removeDCNormalize( src, filtSrc ); + + m_FiltFilt->process( filtSrc, filtDst, m_length ); + + medianFilter( filtDst, dst ); +} + + +void DFProcess::medianFilter(double *src, double *dst) +{ + unsigned int i,k,j,l; + unsigned int index = 0; + + double val = 0; + + double* y = new double[ m_winPost + m_winPre + 1]; + memset( y, 0, sizeof( double ) * ( m_winPost + m_winPre + 1) ); + + double* scratch = new double[ m_length ]; + + for( i = 0; i < m_winPre; i++) + { + k = i + m_winPost + 1; + + for( j = 0; j < k; j++) + { + y[ j ] = src[ j ]; + } + scratch[ index ] = MathUtilities::median( y, k ); + index++; + } + + for( i = 0; i < ( m_length - ( m_winPost + m_winPre ) ); i ++) + { + + l = 0; + for( j = i; j < ( i + m_winPost + m_winPre + 1); j++) + { + y[ l ] = src[ j ]; + l++; + } + + scratch[ index++ ] = MathUtilities::median( y, (m_winPost + m_winPre + 1 )); + } + + for( i = std::max( m_length - m_winPost, (unsigned)1); i < m_length; i++) + { + k = std::max( i - m_winPre, (unsigned)1); + + l = 0; + for( j = k; j < m_length; j++) + { + y[ l ] = src[ j ]; + + l++; + } + + scratch[ index++ ] = MathUtilities::median( y, l); + } + + + for( i = 0; i < m_length; i++ ) + { + val = src[ i ] - scratch[ i ];// - 0.033; + + if( m_isMedianPositive ) + { + if( val > 0 ) + { + dst[ i ] = val; + } + else + { + dst[ i ] = 0; + } + } + else + { + dst[ i ] = val; + } + } + + delete [] y; + delete [] scratch; +} + + +void DFProcess::removeDCNormalize( double *src, double*dst ) +{ + double DFmax = 0; + double DFMin = 0; + double DFAlphaNorm = 0; + + MathUtilities::getFrameMinMax( src, m_length, &DFMin, &DFmax ); + + MathUtilities::getAlphaNorm( src, m_length, m_alphaNormParam, &DFAlphaNorm ); + + for( unsigned int i = 0; i< m_length; i++) + { + dst[ i ] = ( src[ i ] - DFMin ) / DFAlphaNorm; + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/signalconditioning/DFProcess.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,63 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#ifndef CDFPROCESS_H +#define CDFPROCESS_H + +#include <stdio.h> +#include "FiltFilt.h" + +struct DFProcConfig{ + unsigned int length; + unsigned int LPOrd; + double *LPACoeffs; + double *LPBCoeffs; + unsigned int winPre; + unsigned int winPost; + double AlphaNormParam; + bool isMedianPositive;}; + +class DFProcess +{ +public: + DFProcess( DFProcConfig Config ); + virtual ~DFProcess(); + + void process( double* src, double* dst ); + + +private: + void initialise( DFProcConfig Config ); + void deInitialise(); + void removeDCNormalize( double *src, double*dst ); + void medianFilter( double* src, double* dst ); + + unsigned int m_length; + unsigned int m_FFOrd; + + unsigned int m_winPre; + unsigned int m_winPost; + + double m_alphaNormParam; + + double* filtSrc; + double* filtDst; + + double* m_filtScratchIn; + double* m_filtScratchOut; + + FiltFiltConfig m_FilterConfigParams; + + FiltFilt* m_FiltFilt; + + bool m_isMedianPositive; +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/signalconditioning/FiltFilt.cpp Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,124 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#include "FiltFilt.h" + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +FiltFilt::FiltFilt( FiltFiltConfig Config ) +{ + m_filtScratchIn = NULL; + m_filtScratchOut = NULL; + m_ord = 0; + + initialise( Config ); +} + +FiltFilt::~FiltFilt() +{ + deInitialise(); +} + +void FiltFilt::initialise( FiltFiltConfig 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; + + unsigned int nFilt = m_ord + 1; + unsigned int nFact = 3 * ( nFilt - 1); + unsigned int nExt = length + 2 * nFact; + + + m_filtScratchIn = new double[ nExt ]; + m_filtScratchOut = new double[ nExt ]; + + + for( i = 0; i< nExt; i++ ) + { + m_filtScratchIn[ i ] = 0.0; + m_filtScratchOut[ i ] = 0.0; + } + + // Edge transients reflection + double sample0 = 2 * src[ 0 ]; + double sampleN = 2 * src[ length - 1 ]; + + unsigned int index = 0; + for( i = nFact; i > 0; i-- ) + { + m_filtScratchIn[ index++ ] = sample0 - src[ i ]; + } + index = 0; + for( i = 0; i < nFact; i++ ) + { + m_filtScratchIn[ (nExt - nFact) + index++ ] = sampleN - src[ (length - 2) - i ]; + } + + index = 0; + for( i = 0; i < length; i++ ) + { + m_filtScratchIn[ i + nFact ] = src[ i ]; + } + + //////////////////////////////// + // Do 0Ph filtering + m_filter->process( m_filtScratchIn, m_filtScratchOut, nExt); + + // reverse the series for FILTFILT + for ( i = 0; i < nExt; i++) + { + m_filtScratchIn[ i ] = m_filtScratchOut[ nExt - i - 1]; + } + + // do FILTER again + m_filter->process( m_filtScratchIn, m_filtScratchOut, nExt); + + // reverse the series back + for ( i = 0; i < nExt; i++) + { + m_filtScratchIn[ i ] = m_filtScratchOut[ nExt - i - 1 ]; + } + for ( i = 0;i < nExt; i++) + { + m_filtScratchOut[ i ] = m_filtScratchIn[ i ]; + } + + index = 0; + for( i = 0; i < length; i++ ) + { + dst[ index++ ] = m_filtScratchOut[ i + nFact ]; + } + + delete [] m_filtScratchIn; + delete [] m_filtScratchOut; + +} + +void FiltFilt::reset() +{ + +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/signalconditioning/FiltFilt.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,45 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#ifndef FILTFILT_H +#define FILTFILT_H + +#include "Filter.h" + +struct FiltFiltConfig{ + unsigned int ord; + double* ACoeffs; + double* BCoeffs; +}; + +class FiltFilt +{ +public: + FiltFilt( FiltFiltConfig Config ); + virtual ~FiltFilt(); + + void reset(); + void process( double* src, double* dst, unsigned int length ); + +private: + void initialise( FiltFiltConfig Config ); + void deInitialise(); + + unsigned int m_ord; + + Filter* m_filter; + + double* m_filtScratchIn; + double* m_filtScratchOut; + + FilterConfig m_filterConfig; +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/signalconditioning/Filter.cpp Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,82 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#include "Filter.h" + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +Filter::Filter( FilterConfig Config ) +{ + m_ord = 0; + m_outBuffer = NULL; + m_inBuffer = NULL; + + initialise( Config ); +} + +Filter::~Filter() +{ + deInitialise(); +} + +void Filter::initialise( FilterConfig Config ) +{ + m_ord = Config.ord; + m_ACoeffs = Config.ACoeffs; + m_BCoeffs = Config.BCoeffs; + + m_inBuffer = new double[ m_ord + 1 ]; + m_outBuffer = new double[ m_ord + 1 ]; + + reset(); +} + +void Filter::deInitialise() +{ + delete[] m_inBuffer; + delete[] m_outBuffer; +} + +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 */ +} + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/signalconditioning/Filter.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,48 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#ifndef FILTER_H +#define FILTER_H + +#ifndef NULL +#define NULL 0 +#endif + +struct FilterConfig{ + unsigned int ord; + double* ACoeffs; + double* BCoeffs; +}; + +class Filter +{ +public: + Filter( FilterConfig Config ); + virtual ~Filter(); + + void reset(); + + void process( double *src, double *dst, unsigned int length ); + + +private: + void initialise( FilterConfig Config ); + void deInitialise(); + + unsigned int m_ord; + + double* m_inBuffer; + double* m_outBuffer; + + double* m_ACoeffs; + double* m_BCoeffs; +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/signalconditioning/Framer.cpp Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,104 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#include "Framer.h" +#include <math.h> + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +Framer::Framer() +{ + m_dataFrame = NULL; + m_strideFrame = NULL; +} + +Framer::~Framer() +{ + if( m_dataFrame != NULL ) + delete [] m_dataFrame; + + if( m_strideFrame != NULL ) + delete [] m_strideFrame; +} + +void Framer::configure( unsigned int frameLength, unsigned int hop ) +{ + m_frameLength = frameLength; + m_stepSize = hop; + + resetCounters(); + + if( m_dataFrame != NULL ) + { + delete [] m_dataFrame; + m_dataFrame = NULL; + } + m_dataFrame = new double[ m_frameLength ]; + + if( m_strideFrame != NULL ) + { + delete [] m_strideFrame; + m_strideFrame = NULL; + } + m_strideFrame = new double[ m_stepSize ]; +} + +void Framer::getFrame(double *dst) +{ + + if( (m_ulSrcIndex + ( m_frameLength) ) < m_ulSampleLen ) + { + for( unsigned int u = 0; u < m_frameLength; u++) + { + dst[ u ] = m_srcBuffer[ m_ulSrcIndex++ ]; + } + m_ulSrcIndex -= ( m_frameLength - m_stepSize ); + } + else + { + unsigned int rem = (m_ulSampleLen - m_ulSrcIndex ); + unsigned int zero = m_frameLength - rem; + + for( unsigned int u = 0; u < rem; u++ ) + { + dst[ u ] = m_srcBuffer[ m_ulSrcIndex++ ]; + } + + for( unsigned int u = 0; u < zero; u++ ) + { + dst[ rem + u ] = 0; + } + + m_ulSrcIndex -= (( rem - m_stepSize ) ); + } + + m_framesRead++; +} + +void Framer::resetCounters() +{ + m_framesRead = 0; + m_ulSrcIndex = 0; +} + +unsigned int Framer::getMaxNoFrames() +{ + return m_maxFrames; +} + +void Framer::setSource(double *src, unsigned int length) +{ + m_srcBuffer = src; + m_ulSampleLen = length; + + m_maxFrames = (unsigned int)ceil( (double)m_ulSampleLen/(double)m_stepSize ) ; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/signalconditioning/Framer.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,47 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#ifndef FRAMER_H +#define FRAMER_H + +//#include <io.h> +#include <fcntl.h> +#include <stdio.h> + + +class Framer +{ +public: + void setSource( double* src, unsigned int length ); + unsigned int getMaxNoFrames(); + void getFrame( double* dst ); + void configure( unsigned int frameLength, unsigned int hop ); + Framer(); + virtual ~Framer(); + + void resetCounters(); + +private: + + unsigned long m_ulSampleLen; // DataLength (samples) + unsigned int m_framesRead; // Read Frames Index + + double* m_srcBuffer; + double* m_dataFrame; // Analysis Frame Buffer + double* m_strideFrame; // Stride Frame Buffer + unsigned int m_frameLength; // Analysis Frame Length + unsigned int m_stepSize; // Analysis Frame Stride + + unsigned int m_maxFrames; + + unsigned long m_ulSrcIndex; +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/tempotracking/TempoTrack.cpp Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,774 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#include "TempoTrack.h" + +#include "dsp/maths/MathAliases.h" +#include "dsp/maths/MathUtilities.h" + + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +TempoTrack::TempoTrack( TTParams Params ) +{ + m_tempoScratch = NULL; + m_rawDFFrame = NULL; + m_smoothDFFrame = NULL; + m_frameACF = NULL; + + m_dataLength = 0; + m_winLength = 0; + m_lagLength = 0; + + m_rayparam = 0; + m_sigma = 0; + m_DFWVNnorm = 0; + + initialise( Params ); +} + +TempoTrack::~TempoTrack() +{ + deInitialise(); +} + +void TempoTrack::initialise( TTParams Params ) +{ + m_winLength = Params.winLength; + m_lagLength = Params.lagLength; + + m_rayparam = 43.0; + m_sigma = sqrt(3.9017); + m_DFWVNnorm = exp( ( log( 2.0 ) / m_rayparam ) * ( m_winLength + 2 ) ); + + m_rawDFFrame = new double[ m_winLength ]; + m_smoothDFFrame = new double[ m_winLength ]; + m_frameACF = new double[ m_winLength ]; + m_tempoScratch = 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; + m_DFPParams.AlphaNormParam = Params.alpha; + m_DFPParams.LPOrd = Params.LPOrd; + m_DFPParams.LPACoeffs = Params.LPACoeffs; + m_DFPParams.LPBCoeffs = Params.LPBCoeffs; + m_DFPParams.winPre = Params.WinT.pre; + m_DFPParams.winPost = Params.WinT.post; + m_DFPParams.isMedianPositive = true; + + m_DFConditioning = new DFProcess( m_DFPParams ); + +} + +void TempoTrack::deInitialise() +{ + delete [] m_rawDFFrame; + + delete [] m_smoothDFFrame; + + delete [] m_frameACF; + + delete [] m_tempoScratch; + + delete m_DFConditioning; +} + +void TempoTrack::createCombFilter(double* Filter, unsigned int winLength, unsigned int TSig, double beatLag) +{ + unsigned int i; + + if( beatLag == 0 ) + { + for( i = 0; i < winLength; i++ ) + { + Filter[ i ] = ( ( i + 1 ) / pow( m_rayparam, 2.0) ) * exp( ( -pow(( i + 1 ),2.0 ) / ( 2.0 * pow( m_rayparam, 2.0)))); + } + } + else + { + m_sigma = beatLag/8; + for( i = 0; i < winLength; i++ ) + { + double dlag = (double)(i+1) - beatLag; + Filter[ i ] = exp(-0.5 * pow(( dlag / m_sigma), 2.0) ) / (sqrt( 2 * PI) * m_sigma); + } + } +} + +double TempoTrack::tempoMM(double* ACF, double* weight, int tsig) +{ + + double period = 0; + double maxValRCF = 0.0; + unsigned int maxIndexRCF = 0; + + double* pdPeaks; + + unsigned int maxIndexTemp; + double maxValTemp; + unsigned int count; + + unsigned int numelem; + int i, a, b; + + for( i = 0; i < m_lagLength; i++ ) + m_tempoScratch[ i ] = 0.0; + + if( tsig == 0 ) + { + //if time sig is unknown, use metrically unbiased version of Filterbank + numelem = 4; + } + else + { + numelem = tsig; + } + + for(i=1;i<m_lagLength-1;i++) + { + //first and last output values are left intentionally as zero + for (a=1;a<=numelem;a++) + { + for(b=(1-a);b<a;b++) + { + if( tsig == 0 ) + { + m_tempoScratch[i] += ACF[a*(i+1)+b-1] * (1.0 / (2.0 * (double)a-1)) * weight[i]; + } + else + { + m_tempoScratch[i] += ACF[a*(i+1)+b-1] * 1 * weight[i]; + } + } + } + } + + + //NOW FIND MAX INDEX OF ACFOUT + for( i = 0; i < m_lagLength; i++) + { + if( m_tempoScratch[ i ] > maxValRCF) + { + maxValRCF = m_tempoScratch[ i ]; + maxIndexRCF = i; + } + } + + if( tsig == 0 ) + tsig = 4; + + + if( tsig == 4 ) + { + pdPeaks = new double[ 4 ]; + for( i = 0; i < 4; i++ ){ pdPeaks[ i ] = 0.0;} + + pdPeaks[ 0 ] = ( double )maxIndexRCF + 1; + + maxIndexTemp = 0; + maxValTemp = 0.0; + count = 0; + + for( i = (2 * maxIndexRCF + 1) - 1; i < (2 * maxIndexRCF + 1) + 2; i++ ) + { + if( ACF[ i ] > maxValTemp ) + { + maxValTemp = ACF[ i ]; + maxIndexTemp = count; + } + count++; + } + pdPeaks[ 1 ] = (double)( maxIndexTemp + 1 + ( (2 * maxIndexRCF + 1 ) - 2 ) + 1 )/2; + + maxIndexTemp = 0; + maxValTemp = 0.0; + count = 0; + + for( i = (3 * maxIndexRCF + 2 ) - 2; i < (3 * maxIndexRCF + 2 ) + 3; i++ ) + { + if( ACF[ i ] > maxValTemp ) + { + maxValTemp = ACF[ i ]; + maxIndexTemp = count; + } + count++; + } + pdPeaks[ 2 ] = (double)( maxIndexTemp + 1 + ( (3 * maxIndexRCF + 2) - 4 ) + 1 )/3; + + maxIndexTemp = 0; + maxValTemp = 0.0; + count = 0; + + for( i = ( 4 * maxIndexRCF + 3) - 3; i < ( 4 * maxIndexRCF + 3) + 4; i++ ) + { + if( ACF[ i ] > maxValTemp ) + { + maxValTemp = ACF[ i ]; + maxIndexTemp = count; + } + count++; + } + pdPeaks[ 3 ] = (double)( maxIndexTemp + 1 + ( (4 * maxIndexRCF + 3) - 9 ) + 1 )/4 ; + + + period = MathUtilities::mean( pdPeaks, 4 ); + } + else + { + pdPeaks = new double[ 3 ]; + for( i = 0; i < 3; i++ ){ pdPeaks[ i ] = 0.0;} + + pdPeaks[ 0 ] = ( double )maxIndexRCF + 1; + + maxIndexTemp = 0; + maxValTemp = 0.0; + count = 0; + + for( i = (2 * maxIndexRCF + 1) - 1; i < (2 * maxIndexRCF + 1) + 2; i++ ) + { + if( ACF[ i ] > maxValTemp ) + { + maxValTemp = ACF[ i ]; + maxIndexTemp = count; + } + count++; + } + pdPeaks[ 1 ] = (double)( maxIndexTemp + 1 + ( (2 * maxIndexRCF + 1 ) - 2 ) + 1 )/2; + + maxIndexTemp = 0; + maxValTemp = 0.0; + count = 0; + + for( i = (3 * maxIndexRCF + 2 ) - 2; i < (3 * maxIndexRCF + 2 ) + 3; i++ ) + { + if( ACF[ i ] > maxValTemp ) + { + maxValTemp = ACF[ i ]; + maxIndexTemp = count; + } + count++; + } + pdPeaks[ 2 ] = (double)( maxIndexTemp + 1 + ( (3 * maxIndexRCF + 2) - 4 ) + 1 )/3; + + + period = MathUtilities::mean( pdPeaks, 3 ); + } + + delete [] pdPeaks; + + return period; +} + +void TempoTrack::stepDetect( double* periodP, double* periodG, int currentIdx, int* flag ) +{ + double stepthresh = 1 * 3.9017; + + if( *flag ) + { + if(abs(periodG[ currentIdx ] - periodP[ currentIdx ]) > stepthresh) + { + // do nuffin' + } + } + else + { + if(fabs(periodG[ currentIdx ]-periodP[ currentIdx ]) > stepthresh) + { + *flag = 3; + } + } +} + +void TempoTrack::constDetect( double* periodP, int currentIdx, int* flag ) +{ + double constthresh = 2 * 3.9017; + + if( fabs( 2 * periodP[ currentIdx ] - periodP[ currentIdx - 1] - periodP[ currentIdx - 2] ) < constthresh) + { + *flag = 1; + } + else + { + *flag = 0; + } +} + +int TempoTrack::findMeter(double *ACF, unsigned int len, double period) +{ + int i; + int p = (int)MathUtilities::round( period ); + int tsig; + + double Energy_3 = 0.0; + double Energy_4 = 0.0; + + double temp3A = 0.0; + double temp3B = 0.0; + double temp4A = 0.0; + double temp4B = 0.0; + + double* dbf = new double[ len ]; int t = 0; + for( unsigned int u = 0; u < len; u++ ){ dbf[ u ] = 0.0; } + + if( (double)len < 6 * p + 2 ) + { + for( i = ( 3 * p - 2 ); i < ( 3 * p + 2 ) + 1; i++ ) + { + temp3A += ACF[ i ]; + dbf[ t++ ] = ACF[ i ]; + } + + for( i = ( 4 * p - 2 ); i < ( 4 * p + 2 ) + 1; i++ ) + { + temp4A += ACF[ i ]; + } + + Energy_3 = temp3A; + Energy_4 = temp4A; + } + else + { + for( i = ( 3 * p - 2 ); i < ( 3 * p + 2 ) + 1; i++ ) + { + temp3A += ACF[ i ]; + } + + for( i = ( 4 * p - 2 ); i < ( 4 * p + 2 ) + 1; i++ ) + { + temp4A += ACF[ i ]; + } + + for( i = ( 6 * p - 2 ); i < ( 6 * p + 2 ) + 1; i++ ) + { + temp3B += ACF[ i ]; + } + + for( i = ( 2 * p - 2 ); i < ( 2 * p + 2 ) + 1; i++ ) + { + temp4B += ACF[ i ]; + } + + Energy_3 = temp3A + temp3B; + Energy_4 = temp4A + temp4B; + } + + if (Energy_3 > Energy_4) + { + tsig = 3; + } + else + { + tsig = 4; + } + + + return tsig; +} + +void TempoTrack::createPhaseExtractor(double *Filter, unsigned int winLength, double period, unsigned int fsp, unsigned int lastBeat) +{ + int p = (int)MathUtilities::round( period ); + int predictedOffset = 0; + + double* phaseScratch = new double[ p*2 ]; + + + if( lastBeat != 0 ) + { + lastBeat = (int)MathUtilities::round((double)lastBeat );///(double)winLength); + + predictedOffset = lastBeat + p - fsp; + + if (predictedOffset < 0) + { + lastBeat = 0; + } + } + + if( lastBeat != 0 ) + { + int mu = p; + double sigma = (double)p/4; + double PhaseMin = 0.0; + double PhaseMax = 0.0; + unsigned int scratchLength = p*2; + double temp = 0.0; + + for( int i = 0; i < scratchLength; i++ ) + { + phaseScratch[ i ] = exp( -0.5 * pow( ( i - mu ) / sigma, 2 ) ) / ( sqrt( 2*PI ) *sigma ); + } + + MathUtilities::getFrameMinMax( phaseScratch, scratchLength, &PhaseMin, &PhaseMax ); + + for(int i = 0; i < scratchLength; i ++) + { + temp = phaseScratch[ i ]; + phaseScratch[ i ] = (temp - PhaseMin)/PhaseMax; + } + + unsigned int index = 0; + for(int i = p - ( predictedOffset - 1); i < p + ( p - predictedOffset) + 1; i++) + { + Filter[ index++ ] = phaseScratch[ i ]; + } + } + else + { + for( int i = 0; i < p; i ++) + { + Filter[ i ] = 1; + } + } + + delete [] phaseScratch; +} + +int TempoTrack::phaseMM(double *DF, double *weighting, unsigned int winLength, double period) +{ + int alignment = 0; + int p = (int)MathUtilities::round( period ); + + double temp = 0.0; + + double* y = new double[ winLength ]; + double* align = new double[ p ]; + + for( int i = 0; i < winLength; i++ ) + { + y[ i ] = (double)( -i + winLength )/(double)winLength; + } + + for( int o = 0; o < p; o++ ) + { + temp = 0.0; + for(int i = 1 + (o - 1); i< winLength; i += (p + 1)) + { + temp = temp + DF[ i ] * y[ i ]; + } + align[ o ] = temp * weighting[ o ]; + } + + + double valTemp = 0.0; + for(int i = 0; i < p; i++) + { + if( align[ i ] > valTemp ) + { + valTemp = align[ i ]; + alignment = i; + } + } + + delete [] y; + delete [] align; + + return alignment; +} + +int TempoTrack::beatPredict(unsigned int FSP0, double alignment, double period, unsigned int step ) +{ + int beat = 0; + + int p = (int)MathUtilities::round( period ); + int align = (int)MathUtilities::round( alignment ); + int FSP = (int)MathUtilities::round( FSP0 ); + + int FEP = FSP + ( step ); + + beat = FSP + align; + + m_beats.push_back( beat ); + + while( beat + p < FEP ) + { + beat += p; + + m_beats.push_back( beat ); + } + + return beat; +} + +vector<int> TempoTrack::process(double *DF, unsigned int length) +{ + m_dataLength = length; + + double period = 0.0; + int stepFlag = 0; + int constFlag = 0; + int FSP = 0; + int tsig = 0; + int lastBeat = 0; + + + double* RW = new double[ m_lagLength ]; + for( unsigned 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;} + + double* PW = new double[ m_lagLength ]; + for(unsigned int clear = 0; clear < m_lagLength; clear++){ PW[ clear ] = 0.0;} + + m_DFFramer.setSource( DF, m_dataLength ); + + unsigned int TTFrames = m_DFFramer.getMaxNoFrames(); + + double* periodP = new double[ TTFrames ]; + for(unsigned int clear = 0; clear < TTFrames; clear++){ periodP[ clear ] = 0.0;} + + double* periodG = new double[ TTFrames ]; + for(unsigned int clear = 0; clear < TTFrames; clear++){ periodG[ clear ] = 0.0;} + + double* alignment = new double[ TTFrames ]; + for(unsigned int clear = 0; clear < TTFrames; clear++){ alignment[ clear ] = 0.0;} + + m_beats.clear(); + + createCombFilter( RW, m_lagLength, 0, 0 ); + + int TTLoopIndex = 0; + + for( unsigned int i = 0; i < TTFrames; i++ ) + { + m_DFFramer.getFrame( m_rawDFFrame ); + + m_DFConditioning->process( m_rawDFFrame, m_smoothDFFrame ); + + m_correlator.doAutoUnBiased( m_smoothDFFrame, m_frameACF, m_winLength ); + + periodP[ TTLoopIndex ] = tempoMM( m_frameACF, RW, 0 ); + + if( GW[ 0 ] != 0 ) + { + periodG[ TTLoopIndex ] = tempoMM( m_frameACF, GW, tsig ); + } + else + { + periodG[ TTLoopIndex ] = 0.0; + } + + stepDetect( periodP, periodG, TTLoopIndex, &stepFlag ); + + if( stepFlag == 1) + { + constDetect( periodP, TTLoopIndex, &constFlag ); + stepFlag = 0; + } + else + { + stepFlag -= 1; + } + + if( stepFlag < 0 ) + { + stepFlag = 0; + } + + if( constFlag != 0) + { + tsig = findMeter( m_frameACF, m_winLength, periodP[ TTLoopIndex ] ); + + createCombFilter( GW, m_lagLength, tsig, periodP[ TTLoopIndex ] ); + + periodG[ TTLoopIndex ] = tempoMM( m_frameACF, GW, tsig ); + + period = periodG[ TTLoopIndex ]; + + createPhaseExtractor( PW, m_winLength, period, FSP, 0 ); + + constFlag = 0; + + } + else + { + if( GW[ 0 ] != 0 ) + { + period = periodG[ TTLoopIndex ]; + createPhaseExtractor( PW, m_winLength, period, FSP, lastBeat ); + + } + else + { + period = periodP[ TTLoopIndex ]; + createPhaseExtractor( PW, m_winLength, period, FSP, 0 ); + } + } + + alignment[ TTLoopIndex ] = phaseMM( m_rawDFFrame, PW, m_winLength, period ); + + lastBeat = beatPredict(FSP, alignment[ TTLoopIndex ], period, m_lagLength ); + + FSP += (m_lagLength); + + TTLoopIndex++; + } + + + delete [] periodP; + delete [] periodG; + delete [] alignment; + + delete [] RW; + delete [] GW; + delete [] PW; + + return m_beats; +} + + + + + +vector<int> TempoTrack::process( vector <double> DF ) +{ + m_dataLength = DF.size(); + + double period = 0.0; + int stepFlag = 0; + int constFlag = 0; + int FSP = 0; + int tsig = 0; + int lastBeat = 0; + + vector <double> causalDF; + + causalDF = DF; + + //Prepare Causal Extension DFData + unsigned int DFCLength = m_dataLength + m_winLength; + + for( unsigned 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;} + + double* GW = new double[ m_lagLength ]; + for(unsigned 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;} + + m_DFFramer.setSource( &causalDF[0], m_dataLength ); + + unsigned int TTFrames = m_DFFramer.getMaxNoFrames(); + + double* periodP = new double[ TTFrames ]; + for(unsigned 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;} + + double* alignment = new double[ TTFrames ]; + for(unsigned clear = 0; clear < TTFrames; clear++){ alignment[ clear ] = 0.0;} + + m_beats.clear(); + + createCombFilter( RW, m_lagLength, 0, 0 ); + + int TTLoopIndex = 0; + + for( unsigned int i = 0; i < TTFrames; i++ ) + { + m_DFFramer.getFrame( m_rawDFFrame ); + + m_DFConditioning->process( m_rawDFFrame, m_smoothDFFrame ); + + m_correlator.doAutoUnBiased( m_smoothDFFrame, m_frameACF, m_winLength ); + + periodP[ TTLoopIndex ] = tempoMM( m_frameACF, RW, 0 ); + + if( GW[ 0 ] != 0 ) + { + periodG[ TTLoopIndex ] = tempoMM( m_frameACF, GW, tsig ); + } + else + { + periodG[ TTLoopIndex ] = 0.0; + } + + stepDetect( periodP, periodG, TTLoopIndex, &stepFlag ); + + if( stepFlag == 1) + { + constDetect( periodP, TTLoopIndex, &constFlag ); + stepFlag = 0; + } + else + { + stepFlag -= 1; + } + + if( stepFlag < 0 ) + { + stepFlag = 0; + } + + if( constFlag != 0) + { + tsig = findMeter( m_frameACF, m_winLength, periodP[ TTLoopIndex ] ); + + createCombFilter( GW, m_lagLength, tsig, periodP[ TTLoopIndex ] ); + + periodG[ TTLoopIndex ] = tempoMM( m_frameACF, GW, tsig ); + + period = periodG[ TTLoopIndex ]; + + createPhaseExtractor( PW, m_winLength, period, FSP, 0 ); + + constFlag = 0; + + } + else + { + if( GW[ 0 ] != 0 ) + { + period = periodG[ TTLoopIndex ]; + createPhaseExtractor( PW, m_winLength, period, FSP, lastBeat ); + + } + else + { + period = periodP[ TTLoopIndex ]; + createPhaseExtractor( PW, m_winLength, period, FSP, 0 ); + } + } + + alignment[ TTLoopIndex ] = phaseMM( m_rawDFFrame, PW, m_winLength, period ); + + lastBeat = beatPredict(FSP, alignment[ TTLoopIndex ], period, m_lagLength ); + + FSP += (m_lagLength); + + TTLoopIndex++; + } + + + delete [] periodP; + delete [] periodG; + delete [] alignment; + + delete [] RW; + delete [] GW; + delete [] PW; + + return m_beats; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/tempotracking/TempoTrack.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,97 @@ +/* -*- 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 copyright 2005-2006 Christian Landone. + All rights reserved. +*/ + +#ifndef TEMPOTRACK_H +#define TEMPOTRACK_H + + +#include <stdio.h> +#include <vector> + +#include "dsp/signalconditioning/DFProcess.h" +#include "dsp/maths/Correlation.h" +#include "dsp/signalconditioning/Framer.h" + + + +using std::vector; + +struct WinThresh +{ + unsigned int pre; + unsigned 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 + 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]: +}; + + +class TempoTrack +{ +public: + TempoTrack( TTParams Params ); + virtual ~TempoTrack(); + + vector<int> process( vector <double> DF ); + vector<int> process( double* DF, unsigned int length ); + + +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 ); + 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 ); + double tempoMM( double* ACF, double* weight, int sig ); + + unsigned int m_dataLength; + unsigned int m_winLength; + unsigned int m_lagLength; + + double m_rayparam; + double m_sigma; + double m_DFWVNnorm; + + vector<int> m_beats; // Vector of detected beats + + double* m_tempoScratch; + + // Processing Buffers + double* m_rawDFFrame; // Original Detection Function Analysis Frame + double* m_smoothDFFrame; // Smoothed Detection Function Analysis Frame + double* m_frameACF; // AutoCorrelation of Smoothed Detection Function + + //Low Pass Coefficients for DF Smoothing + double* m_ACoeffs; + double* m_BCoeffs; + + // Objetcs/operators declaration + Framer m_DFFramer; + DFProcess* m_DFConditioning; + Correlation m_correlator; + + // Config structure for DFProcess + DFProcConfig m_DFPParams; +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/tonal/ChangeDetectionFunction.cpp Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,139 @@ +/* -*- 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 copyright 2006 Martin Gasser. + All rights reserved. +*/ + +#include "ChangeDetectionFunction.h" + +#ifndef PI +#define PI (3.14159265358979232846) +#endif + + + +ChangeDetectionFunction::ChangeDetectionFunction(ChangeDFConfig config) : + m_dFilterSigma(0.0), m_iFilterWidth(0) +{ + setFilterWidth(config.smoothingWidth); +} + +ChangeDetectionFunction::~ChangeDetectionFunction() +{ +} + +void ChangeDetectionFunction::setFilterWidth(const int iWidth) +{ + m_iFilterWidth = iWidth*2+1; + + // it is assumed that the gaussian is 0 outside of +/- FWHM + // => filter width = 2*FWHM = 2*2.3548*sigma + m_dFilterSigma = double(m_iFilterWidth) / double(2*2.3548); + m_vaGaussian.resize(m_iFilterWidth); + + double dScale = 1.0 / (m_dFilterSigma*sqrt(2*PI)); + + for (int x = -(m_iFilterWidth-1)/2; x <= (m_iFilterWidth-1)/2; x++) + { + double w = dScale * std::exp ( -(x*x)/(2*m_dFilterSigma*m_dFilterSigma) ); + m_vaGaussian[x + (m_iFilterWidth-1)/2] = w; + } + +#ifdef DEBUG_CHANGE_DETECTION_FUNCTION + std::cout << "Filter sigma: " << m_dFilterSigma << std::endl; + std::cout << "Filter width: " << m_iFilterWidth << std::endl; +#endif +} + + +ChangeDistance ChangeDetectionFunction::process(const TCSGram& rTCSGram) +{ + ChangeDistance retVal; + retVal.resize(rTCSGram.getSize(), 0.0); + + TCSGram smoothedTCSGram; + + for (int iPosition = 0; iPosition < rTCSGram.getSize(); iPosition++) + { + int iSkipLower = 0; + + int iLowerPos = iPosition - (m_iFilterWidth-1)/2; + int iUpperPos = iPosition + (m_iFilterWidth-1)/2; + + if (iLowerPos < 0) + { + iSkipLower = -iLowerPos; + iLowerPos = 0; + } + + if (iUpperPos >= rTCSGram.getSize()) + { + int iMaxIndex = rTCSGram.getSize() - 1; + iUpperPos = iMaxIndex; + } + + TCSVector smoothedVector; + + // for every bin of the vector, calculate the smoothed value + for (int iPC = 0; iPC < 6; iPC++) + { + size_t j = 0; + double dSmoothedValue = 0.0; + TCSVector rCV; + + for (int i = iLowerPos; i <= iUpperPos; i++) + { + rTCSGram.getTCSVector(i, rCV); + dSmoothedValue += m_vaGaussian[iSkipLower + j++] * rCV[iPC]; + } + + smoothedVector[iPC] = dSmoothedValue; + } + + smoothedTCSGram.addTCSVector(smoothedVector); + } + + for (int iPosition = 0; iPosition < rTCSGram.getSize(); iPosition++) + { + /* + TODO: calculate a confidence measure for the current estimation + if the current estimate is not confident enough, look further into the future/the past + e.g., High frequency content, zero crossing rate, spectral flatness + */ + + TCSVector nextTCS; + TCSVector previousTCS; + + int iWindow = 1; + + // while (previousTCS.magnitude() < 0.1 && (iPosition-iWindow) > 0) + { + smoothedTCSGram.getTCSVector(iPosition-iWindow, previousTCS); + // std::cout << previousTCS.magnitude() << std::endl; + iWindow++; + } + + iWindow = 1; + + // while (nextTCS.magnitude() < 0.1 && (iPosition+iWindow) < (rTCSGram.getSize()-1) ) + { + smoothedTCSGram.getTCSVector(iPosition+iWindow, nextTCS); + iWindow++; + } + + double distance = 0.0; + // Euclidean distance + for (size_t j = 0; j < 6; j++) + { + distance += std::pow(nextTCS[j] - previousTCS[j], 2.0); + } + + retVal[iPosition] = std::pow(distance, 0.5); + } + + return retVal; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/tonal/ChangeDetectionFunction.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,43 @@ +/* -*- 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 copyright 2006 Martin Gasser. + All rights reserved. +*/ + +#ifndef _CHANGEDETECTIONFUNCTION_ +#define _CHANGEDETECTIONFUNCTION_ + +#define DEBUG_CHANGE_DETECTION_FUNCTION 1 + +#include "TCSgram.h" + +#include <valarray> +using std::valarray; + +typedef valarray<double> ChangeDistance; + +struct ChangeDFConfig +{ + int smoothingWidth; +}; + +class ChangeDetectionFunction +{ +public: + ChangeDetectionFunction(ChangeDFConfig); + ~ChangeDetectionFunction(); + ChangeDistance process(const TCSGram& rTCSGram); +private: + void setFilterWidth(const int iWidth); + +private: + valarray<double> m_vaGaussian; + double m_dFilterSigma; + int m_iFilterWidth; +}; + +#endif // _CHANGDETECTIONFUNCTION_
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/tonal/TCSgram.cpp Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,72 @@ +/* -*- 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 copyright 2006 Martin Gasser. + All rights reserved. +*/ + +#include "TCSgram.h" + +#include <valarray> +#include <cmath> +#include <iostream> +#include <limits> + +#include "dsp/maths/MathUtilities.h" + +TCSGram::TCSGram() : + m_uNumBins(6) +{ +} + +TCSGram::~TCSGram() +{ +} + + +void TCSGram::getTCSVector(int iPosition, TCSVector& rTCSVector) const +{ + if (iPosition < 0) + rTCSVector = TCSVector(); + else if (iPosition >= m_VectorList.size()) + rTCSVector = TCSVector(); + else + rTCSVector = m_VectorList[iPosition].second; +} + +long TCSGram::getTime(size_t uPosition) const +{ + return m_VectorList[uPosition].first; +} + + +void TCSGram::addTCSVector(const TCSVector& rTCSVector) +{ + size_t uSize = m_VectorList.size(); + long lMilliSeconds = static_cast<long>(uSize*m_dFrameDurationMS); + std::pair<long, TCSVector> p; + p.first = lMilliSeconds; + p.second = rTCSVector; + + m_VectorList.push_back(p); +} + +long TCSGram::getDuration() const +{ + size_t uSize = m_VectorList.size(); + return static_cast<long>(uSize*m_dFrameDurationMS); +} + +void TCSGram::printDebug() +{ + vectorlist_t::iterator vectorIterator = m_VectorList.begin(); + + while (vectorIterator != m_VectorList.end()) + { + vectorIterator->second.printDebug(); + vectorIterator++; + } +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/tonal/TCSgram.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,44 @@ +/* -*- 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 copyright 2006 Martin Gasser. + All rights reserved. +*/ + +#ifndef _TCSGram_ +#define _TCSGram_ + +#include <vector> +#include <valarray> +#include <utility> + +#include "TonalEstimator.h" + +typedef std::vector<std::pair<long, TCSVector> > vectorlist_t; + +class TCSGram +{ +public: + TCSGram(); + ~TCSGram(); + void getTCSVector(int, TCSVector&) const; + void addTCSVector(const TCSVector&); + long getTime(size_t) const; + long getDuration() const; + void printDebug(); + int getSize() const { return m_VectorList.size(); } + void reserve(size_t uSize) { m_VectorList.reserve(uSize); } + void clear() { m_VectorList.clear(); } + void setFrameDuration(const double dFrameDurationMS) { m_dFrameDurationMS = dFrameDurationMS; } + void setNumBins(const unsigned int uNumBins) { m_uNumBins = uNumBins; } + void normalize(); +protected: + vectorlist_t m_VectorList; + unsigned int m_uNumBins; + double m_dFrameDurationMS; +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/tonal/TonalEstimator.cpp Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,98 @@ +/* -*- 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 copyright 2006 Martin Gasser. + All rights reserved. +*/ + +#include "TonalEstimator.h" + +#include <cmath> +#include <iostream> + +#ifndef PI +#define PI (3.14159265358979232846) +#endif + +TonalEstimator::TonalEstimator() +{ + m_Basis.resize(6); + + int i = 0; + + + // circle of fifths + m_Basis[i].resize(12); + for (int iP = 0; iP < 12; iP++) + { + m_Basis[i][iP] = std::sin( (7.0 / 6.0) * iP * PI); + } + + i++; + + m_Basis[i].resize(12); + for (int iP = 0; iP < 12; iP++) + { + m_Basis[i][iP] = std::cos( (7.0 / 6.0) * iP * PI); + } + + i++; + + + // circle of major thirds + m_Basis[i].resize(12); + for (int iP = 0; iP < 12; iP++) + { + m_Basis[i][iP] = 0.6 * std::sin( (2.0 / 3.0) * iP * PI); + } + + i++; + + m_Basis[i].resize(12); + for (int iP = 0; iP < 12; iP++) + { + m_Basis[i][iP] = 0.6 * std::cos( (2.0 / 3.0) * iP * PI); + } + + i++; + + + // circle of minor thirds + m_Basis[i].resize(12); + for (int iP = 0; iP < 12; iP++) + { + m_Basis[i][iP] = 1.1 * std::sin( (3.0 / 2.0) * iP * PI); + } + + i++; + + m_Basis[i].resize(12); + for (int iP = 0; iP < 12; iP++) + { + m_Basis[i][iP] = 1.1 * std::cos( (3.0 / 2.0) * iP * PI); + } + +} + +TonalEstimator::~TonalEstimator() +{ +} + +TCSVector TonalEstimator::transform2TCS(const ChromaVector& rVector) +{ + TCSVector vaRetVal; + vaRetVal.resize(6, 0.0); + + for (int i = 0; i < 6; i++) + { + for (int iP = 0; iP < 12; iP++) + { + vaRetVal[i] += m_Basis[i][iP] * rVector[iP]; + } + } + + return vaRetVal; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/tonal/TonalEstimator.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,94 @@ +/* -*- 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 copyright 2006 Martin Gasser. + All rights reserved. +*/ + +#ifndef _TONALESTIMATOR_ +#define _TONALESTIMATOR_ + + +#include <valarray> +#include <numeric> +#include <algorithm> +#include <iostream> + +class ChromaVector : public std::valarray<double> +{ +public: + ChromaVector(size_t uSize = 12) : std::valarray<double>() + { resize(uSize, 0.0f); } + + virtual ~ChromaVector() {}; + + void printDebug() + { + for (int i = 0; i < size(); i++) + { + std::cout << (*this)[i] << ";"; + } + + std::cout << std::endl; + } + + void normalizeL1() + { + // normalize the chroma vector (L1 norm) + double dSum = 0.0; + + for (size_t i = 0; i < 12; (dSum += std::abs((*this)[i++]))); + for (size_t i = 0; i < 12; dSum > 0.0000001?((*this)[i] /= dSum):(*this)[i]=0.0, i++); + + } + +}; + +class TCSVector : public std::valarray<double> +{ +public: + TCSVector() : std::valarray<double>() + { resize(6, 0.0f); } + + virtual ~TCSVector() {}; + + void printDebug() + { + for (int i = 0; i < size(); i++) + { + std::cout << (*this)[i] << ";"; + } + + std::cout << std::endl; + } + + double magnitude() const + { + double dMag = 0.0; + + for (size_t i = 0; i < 6; i++) + { + dMag += std::pow((*this)[i], 2.0); + } + + return std::sqrt(dMag); + } + +}; + + + +class TonalEstimator +{ +public: + TonalEstimator(); + virtual ~TonalEstimator(); + TCSVector transform2TCS(const ChromaVector& rVector); +protected: + std::valarray< std::valarray<double> > m_Basis; +}; + +#endif // _TONALESTIMATOR_
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/transforms/FFT.cpp Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,154 @@ +/* -*- 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 is based on Don Cross's public domain FFT implementation. +*/ + +#include "FFT.h" +#include <cmath> + +////////////////////////////////////////////////////////////////////// +// Construction/Destruction +////////////////////////////////////////////////////////////////////// + +FFT::FFT() +{ + +} + +FFT::~FFT() +{ + +} + +void FFT::process(unsigned int p_nSamples, bool p_bInverseTransform, double *p_lpRealIn, double *p_lpImagIn, double *p_lpRealOut, double *p_lpImagOut) +{ + + if(!p_lpRealIn || !p_lpRealOut || !p_lpImagOut) return; + + + unsigned int NumBits; + unsigned int i, j, k, n; + unsigned int BlockSize, BlockEnd; + + double angle_numerator = 2.0 * M_PI; + double tr, ti; + + if( !isPowerOfTwo(p_nSamples) ) + { + return; + } + + if( p_bInverseTransform ) angle_numerator = -angle_numerator; + + NumBits = numberOfBitsNeeded ( p_nSamples ); + + + for( i=0; i < p_nSamples; i++ ) + { + j = reverseBits ( i, NumBits ); + p_lpRealOut[j] = p_lpRealIn[i]; + p_lpImagOut[j] = (p_lpImagIn == 0) ? 0.0 : p_lpImagIn[i]; + } + + + BlockEnd = 1; + for( BlockSize = 2; BlockSize <= p_nSamples; BlockSize <<= 1 ) + { + double delta_angle = angle_numerator / (double)BlockSize; + double sm2 = -sin ( -2 * delta_angle ); + double sm1 = -sin ( -delta_angle ); + double cm2 = cos ( -2 * delta_angle ); + double cm1 = cos ( -delta_angle ); + double w = 2 * cm1; + double ar[3], ai[3]; + + for( i=0; i < p_nSamples; i += BlockSize ) + { + + ar[2] = cm2; + ar[1] = cm1; + + ai[2] = sm2; + ai[1] = sm1; + + for ( j=i, n=0; n < BlockEnd; j++, n++ ) + { + + ar[0] = w*ar[1] - ar[2]; + ar[2] = ar[1]; + ar[1] = ar[0]; + + ai[0] = w*ai[1] - ai[2]; + ai[2] = ai[1]; + ai[1] = ai[0]; + + k = j + BlockEnd; + tr = ar[0]*p_lpRealOut[k] - ai[0]*p_lpImagOut[k]; + ti = ar[0]*p_lpImagOut[k] + ai[0]*p_lpRealOut[k]; + + p_lpRealOut[k] = p_lpRealOut[j] - tr; + p_lpImagOut[k] = p_lpImagOut[j] - ti; + + p_lpRealOut[j] += tr; + p_lpImagOut[j] += ti; + + } + } + + BlockEnd = BlockSize; + + } + + + if( p_bInverseTransform ) + { + double denom = (double)p_nSamples; + + for ( i=0; i < p_nSamples; i++ ) + { + p_lpRealOut[i] /= denom; + p_lpImagOut[i] /= denom; + } + } +} + +bool FFT::isPowerOfTwo(unsigned int p_nX) +{ + if( p_nX < 2 ) return false; + + if( p_nX & (p_nX-1) ) return false; + + return true; +} + +unsigned int FFT::numberOfBitsNeeded(unsigned int p_nSamples) +{ + int i; + + if( p_nSamples < 2 ) + { + return 0; + } + + for ( i=0; ; i++ ) + { + if( p_nSamples & (1 << i) ) return i; + } +} + +unsigned int FFT::reverseBits(unsigned int p_nIndex, unsigned int p_nBits) +{ + unsigned int i, rev; + + for(i=rev=0; i < p_nBits; i++) + { + rev = (rev << 1) | (p_nIndex & 1); + p_nIndex >>= 1; + } + + return rev; +}
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dsp/transforms/FFT.h Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,28 @@ +/* -*- 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 is based on Don Cross's public domain FFT implementation. +*/ + +#ifndef FFT_H +#define FFT_H + +class FFT +{ +public: + static void process(unsigned int p_nSamples, bool p_bInverseTransform, + double *p_lpRealIn, double *p_lpImagIn, + double *p_lpRealOut, double *p_lpImagOut); + FFT(); + virtual ~FFT(); + +protected: + static unsigned int reverseBits(unsigned int p_nIndex, unsigned int p_nBits); + static unsigned int numberOfBitsNeeded( unsigned int p_nSamples ); + static bool isPowerOfTwo( unsigned int p_nX ); +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/qm-dsp.pro Wed Apr 05 17:35:59 2006 +0000 @@ -0,0 +1,76 @@ +TEMPLATE = lib +CONFIG += release warn_on staticlib +OBJECTS_DIR = tmp_obj +MOC_DIR = tmp_moc +DEPENDPATH += base \ + dsp/chromagram \ + dsp/maths \ + dsp/onsets \ + dsp/phasevocoder \ + dsp/rateconversion \ + dsp/signalconditioning \ + dsp/tempotracking \ + dsp/tonal \ + dsp/transforms +INCLUDEPATH += . \ + base \ + dsp/maths \ + dsp/chromagram \ + dsp/transforms \ + dsp/onsets \ + dsp/phasevocoder \ + dsp/signalconditioning \ + dsp/rateconversion \ + dsp/tempotracking \ + dsp/tonal + +# Input +HEADERS += base/Pitch.h \ + base/Window.h \ + dsp/chromagram/Chromagram.h \ + dsp/chromagram/ChromaProcess.h \ + dsp/chromagram/ConstantQ.h \ + dsp/maths/Correlation.h \ + dsp/maths/Histogram.h \ + dsp/maths/MathAliases.h \ + dsp/maths/MathUtilities.h \ + dsp/maths/Polyfit.h \ + dsp/maths/SimpleMatrix.h \ + dsp/onsets/DetectionFunction.h \ + dsp/onsets/PeakPicking.h \ + dsp/phasevocoder/PhaseVocoder.h \ + dsp/rateconversion/Decimator.h \ + dsp/signalconditioning/DFProcess.h \ + dsp/signalconditioning/Filter.h \ + dsp/signalconditioning/FiltFilt.h \ + dsp/signalconditioning/Framer.h \ + dsp/signalconditioning/Windowing.h \ + dsp/tempotracking/TempoTrack.h \ + dsp/tonal/ChangeDetectionFunction.h \ + dsp/tonal/ChromaMatrix.h \ + dsp/tonal/TCSgram.h \ + dsp/tonal/TonalEstimator.h \ + dsp/transforms/FFT.h \ + dsp/transforms/Fourier.h +SOURCES += base/Pitch.cpp \ + dsp/chromagram/Chromagram.cpp \ + dsp/chromagram/ChromaProcess.cpp \ + dsp/chromagram/ConstantQ.cpp \ + dsp/maths/Correlation.cpp \ + dsp/maths/MathUtilities.cpp \ + dsp/onsets/DetectionFunction.cpp \ + dsp/onsets/PeakPicking.cpp \ + dsp/phasevocoder/PhaseVocoder.cpp \ + dsp/rateconversion/Decimator.cpp \ + dsp/signalconditioning/DFProcess.cpp \ + dsp/signalconditioning/Filter.cpp \ + dsp/signalconditioning/FiltFilt.cpp \ + dsp/signalconditioning/Framer.cpp \ + dsp/signalconditioning/Windowing.cpp \ + dsp/tempotracking/TempoTrack.cpp \ + dsp/tonal/ChangeDetectionFunction.cpp \ + dsp/tonal/ChromaMatrix.cpp \ + dsp/tonal/TCSgram.cpp \ + dsp/tonal/TonalEstimator.cpp \ + dsp/transforms/FFT.cpp \ + dsp/transforms/Fourier.cpp