Mercurial > hg > qm-dsp
changeset 241:a98dd8ec96f8
* Move dsp/maths to maths ; bring PCA and HMM across from Soundbite
| author | Chris Cannam <c.cannam@qmul.ac.uk> |
|---|---|
| date | Wed, 09 Jan 2008 10:31:29 +0000 |
| parents | 1a406914b3a9 |
| children | 6060110dc3c6 |
| files | dsp/chromagram/ChromaProcess.cpp dsp/chromagram/Chromagram.cpp dsp/chromagram/ConstantQ.h dsp/keydetection/GetKeyMode.cpp dsp/maths/Correlation.cpp dsp/maths/Correlation.h dsp/maths/Histogram.h dsp/maths/MathAliases.h dsp/maths/MathUtilities.cpp dsp/maths/MathUtilities.h dsp/maths/Polyfit.h dsp/onsets/DetectionFunction.cpp dsp/onsets/DetectionFunction.h dsp/onsets/PeakPicking.cpp dsp/onsets/PeakPicking.h dsp/signalconditioning/DFProcess.cpp dsp/tempotracking/TempoTrack.cpp dsp/tempotracking/TempoTrack.h dsp/tonal/TCSgram.cpp hmm/hmm.c hmm/hmm.h maths/Correlation.cpp maths/Correlation.h maths/Histogram.h maths/MathAliases.h maths/MathUtilities.cpp maths/MathUtilities.h maths/Polyfit.h maths/pca/data.txt maths/pca/pca.c maths/pca/pca.h qm-dsp.pro |
| diffstat | 32 files changed, 2447 insertions(+), 1215 deletions(-) [+] |
line wrap: on
line diff
--- a/dsp/chromagram/ChromaProcess.cpp Wed Nov 21 16:53:51 2007 +0000 +++ b/dsp/chromagram/ChromaProcess.cpp Wed Jan 09 10:31:29 2008 +0000 @@ -10,7 +10,7 @@ #include "ChromaProcess.h" -#include "dsp/maths/Histogram.h" +#include "maths/Histogram.h" #include <math.h> ////////////////////////////////////////////////////////////////////// // Construction/Destruction
--- a/dsp/chromagram/Chromagram.cpp Wed Nov 21 16:53:51 2007 +0000 +++ b/dsp/chromagram/Chromagram.cpp Wed Jan 09 10:31:29 2008 +0000 @@ -10,7 +10,7 @@ #include <iostream> #include <cmath> -#include "dsp/maths/MathUtilities.h" +#include "maths/MathUtilities.h" #include "Chromagram.h" //----------------------------------------------------------------------------
--- a/dsp/chromagram/ConstantQ.h Wed Nov 21 16:53:51 2007 +0000 +++ b/dsp/chromagram/ConstantQ.h Wed Jan 09 10:31:29 2008 +0000 @@ -12,8 +12,8 @@ #define CONSTANTQ_H #include <vector> -#include "dsp/maths/MathAliases.h" -#include "dsp/maths/MathUtilities.h" +#include "maths/MathAliases.h" +#include "maths/MathUtilities.h" struct CQConfig{ unsigned int FS;
--- a/dsp/keydetection/GetKeyMode.cpp Wed Nov 21 16:53:51 2007 +0000 +++ b/dsp/keydetection/GetKeyMode.cpp Wed Jan 09 10:31:29 2008 +0000 @@ -3,7 +3,7 @@ ////////////////////////////////////////////////////////////////////// #include "GetKeyMode.h" -#include "dsp/maths/MathUtilities.h" +#include "maths/MathUtilities.h" #include "base/Pitch.h" #include <iostream>
--- a/dsp/maths/Correlation.cpp Wed Nov 21 16:53:51 2007 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,51 +0,0 @@ -/* -*- 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; - } -}
--- a/dsp/maths/Correlation.h Wed Nov 21 16:53:51 2007 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,25 +0,0 @@ -/* -*- 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 //
--- a/dsp/maths/Histogram.h Wed Nov 21 16:53:51 2007 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,393 +0,0 @@ -/* -*- 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 -
--- a/dsp/maths/MathAliases.h Wed Nov 21 16:53:51 2007 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,55 +0,0 @@ -/* -*- 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
--- a/dsp/maths/MathUtilities.cpp Wed Nov 21 16:53:51 2007 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,221 +0,0 @@ -/* -*- 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 ; - } - - } -} - -int MathUtilities::getMax( double* pData, unsigned int Length, double* pMax ) -{ - unsigned int index = 0; - unsigned int i; - double temp = 0.0; - - double max = pData[0]; - - for( i = 0; i < Length; i++) - { - temp = pData[ i ]; - - if( temp > max ) - { - max = temp ; - index = i; - } - - } - - *pMax = max; - - - return index; -} - -void MathUtilities::circShift( double* pData, int length, int shift) -{ - shift = shift % length; - double temp; - int i,n; - - for( i = 0; i < shift; i++) - { - temp=*(pData + length - 1); - - for( n = length-2; n >= 0; n--) - { - *(pData+n+1)=*(pData+n); - } - - *pData = temp; - } -} - -int MathUtilities::compareInt (const void * a, const void * b) -{ - return ( *(int*)a - *(int*)b ); -} -
--- a/dsp/maths/MathUtilities.h Wed Nov 21 16:53:51 2007 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,33 +0,0 @@ -/* -*- 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 ); - static void circShift( double* data, int length, int shift); - static int getMax( double* data, unsigned int length, double* max ); - static int compareInt(const void * a, const void * b); -}; - -#endif
--- a/dsp/maths/Polyfit.h Wed Nov 21 16:53:51 2007 +0000 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,398 +0,0 @@ -/* -*- 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 -
--- a/dsp/onsets/DetectionFunction.cpp Wed Nov 21 16:53:51 2007 +0000 +++ b/dsp/onsets/DetectionFunction.cpp Wed Jan 09 10:31:29 2008 +0000 @@ -144,10 +144,6 @@ case DF_BROADBAND: retVal = broadband( m_halfLength, m_magnitude); break; - - case DF_POWER: - retVal = power( m_halfLength, m_magnitude ); - break; } return retVal; @@ -269,15 +265,6 @@ return val; } -double DetectionFunction::power(unsigned int length, double *src) -{ - double val = 0; - for (unsigned int i = 0; i < length; ++i) { - val += src[i]; - } - return val; -} - double* DetectionFunction::getSpectrumMagnitude() { return m_magnitude;
--- a/dsp/onsets/DetectionFunction.h Wed Nov 21 16:53:51 2007 +0000 +++ b/dsp/onsets/DetectionFunction.h Wed Jan 09 10:31:29 2008 +0000 @@ -11,8 +11,8 @@ #ifndef DETECTIONFUNCTION_H #define DETECTIONFUNCTION_H -#include "dsp/maths/MathUtilities.h" -#include "dsp/maths/MathAliases.h" +#include "maths/MathUtilities.h" +#include "maths/MathAliases.h" #include "dsp/phasevocoder/PhaseVocoder.h" #include "base/Window.h" @@ -21,7 +21,6 @@ #define DF_PHASEDEV (3) #define DF_COMPLEXSD (4) #define DF_BROADBAND (5) -#define DF_POWER (6) struct DFConfig{ double stepSecs; // DF step in seconds @@ -52,7 +51,6 @@ double phaseDev(unsigned int length, double *srcPhase); double complexSD(unsigned int length, double *srcMagnitude, double *srcPhase); double broadband(unsigned int length, double *srcMagnitude); - double power(unsigned int length, double *src); private: void initialise( DFConfig Config );
--- a/dsp/onsets/PeakPicking.cpp Wed Nov 21 16:53:51 2007 +0000 +++ b/dsp/onsets/PeakPicking.cpp Wed Jan 09 10:31:29 2008 +0000 @@ -9,7 +9,7 @@ */ #include "PeakPicking.h" -#include "dsp/maths/Polyfit.h" +#include "maths/Polyfit.h" #include <iostream>
--- a/dsp/onsets/PeakPicking.h Wed Nov 21 16:53:51 2007 +0000 +++ b/dsp/onsets/PeakPicking.h Wed Jan 09 10:31:29 2008 +0000 @@ -15,8 +15,8 @@ #ifndef PEAKPICKING_H #define PEAKPICKING_H -#include "dsp/maths/MathUtilities.h" -#include "dsp/maths/MathAliases.h" +#include "maths/MathUtilities.h" +#include "maths/MathAliases.h" #include "dsp/signalconditioning/DFProcess.h"
--- a/dsp/signalconditioning/DFProcess.cpp Wed Nov 21 16:53:51 2007 +0000 +++ b/dsp/signalconditioning/DFProcess.cpp Wed Jan 09 10:31:29 2008 +0000 @@ -9,7 +9,7 @@ */ #include "DFProcess.h" -#include "dsp/maths/MathUtilities.h" +#include "maths/MathUtilities.h" ////////////////////////////////////////////////////////////////////// // Construction/Destruction
--- a/dsp/tempotracking/TempoTrack.cpp Wed Nov 21 16:53:51 2007 +0000 +++ b/dsp/tempotracking/TempoTrack.cpp Wed Jan 09 10:31:29 2008 +0000 @@ -10,8 +10,8 @@ #include "TempoTrack.h" -#include "dsp/maths/MathAliases.h" -#include "dsp/maths/MathUtilities.h" +#include "maths/MathAliases.h" +#include "maths/MathUtilities.h" #include <iostream>
--- a/dsp/tempotracking/TempoTrack.h Wed Nov 21 16:53:51 2007 +0000 +++ b/dsp/tempotracking/TempoTrack.h Wed Jan 09 10:31:29 2008 +0000 @@ -16,7 +16,7 @@ #include <vector> #include "dsp/signalconditioning/DFProcess.h" -#include "dsp/maths/Correlation.h" +#include "maths/Correlation.h" #include "dsp/signalconditioning/Framer.h"
--- a/dsp/tonal/TCSgram.cpp Wed Nov 21 16:53:51 2007 +0000 +++ b/dsp/tonal/TCSgram.cpp Wed Jan 09 10:31:29 2008 +0000 @@ -15,7 +15,7 @@ #include <iostream> #include <limits> -#include "dsp/maths/MathUtilities.h" +#include "maths/MathUtilities.h" TCSGram::TCSGram() : m_uNumBins(6)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/hmm/hmm.c Wed Jan 09 10:31:29 2008 +0000 @@ -0,0 +1,800 @@ +/* + * hmm.c + * soundbite + * + * Created by Mark Levy on 12/02/2006. + * Copyright 2006 Centre for Digital Music, Queen Mary, University of London. All rights reserved. + * + */ + +#include <stdio.h> +#include <math.h> +#include <stdlib.h> +#include <float.h> +#include <time.h> /* to seed random number generator */ +#include "clapack.h" /* LAPACK for matrix inversion */ +#ifdef _MAC_OS_X +#include <vecLib/cblas.h> +#else +#include <cblas.h> /* BLAS for matrix multiplication */ +#endif + +#include "hmm.h" + +model_t* hmm_init(double** x, int T, int L, int N) +{ + int i, j, d, e, t; + double s, ss; + + model_t* model; + model = (model_t*) malloc(sizeof(model_t)); + model->N = N; + model->L = L; + model->p0 = (double*) malloc(N*sizeof(double)); + model->a = (double**) malloc(N*sizeof(double*)); + model->mu = (double**) malloc(N*sizeof(double*)); + for (i = 0; i < N; i++) + { + model->a[i] = (double*) malloc(N*sizeof(double)); + model->mu[i] = (double*) malloc(L*sizeof(double)); + } + model->cov = (double**) malloc(L*sizeof(double*)); + for (i = 0; i < L; i++) + model->cov[i] = (double*) malloc(L*sizeof(double)); + + srand(time(0)); + double* global_mean = (double*) malloc(L*sizeof(double)); + + /* find global mean */ + for (d = 0; d < L; d++) + { + global_mean[d] = 0; + for (t = 0; t < T; t++) + global_mean[d] += x[t][d]; + global_mean[d] /= T; + } + + /* calculate global diagonal covariance */ + for (d = 0; d < L; d++) + { + for (e = 0; e < L; e++) + model->cov[d][e] = 0; + for (t = 0; t < T; t++) + model->cov[d][d] += (x[t][d] - global_mean[d]) * (x[t][d] - global_mean[d]); + model->cov[d][d] /= T-1; + } + + /* set all means close to global mean */ + for (i = 0; i < N; i++) + { + for (d = 0; d < L; d++) + { + /* add some random noise related to covariance */ + /* ideally the random number would be Gaussian(0,1), as a hack we make it uniform on [-0.25,0.25] */ + model->mu[i][d] = global_mean[d] + (0.5 * rand() / (double) RAND_MAX - 0.25) * sqrt(model->cov[d][d]); + } + } + + /* random intial and transition probs */ + s = 0; + for (i = 0; i < N; i++) + { + model->p0[i] = 1 + rand() / (double) RAND_MAX; + s += model->p0[i]; + ss = 0; + for (j = 0; j < N; j++) + { + model->a[i][j] = 1 + rand() / (double) RAND_MAX; + ss += model->a[i][j]; + } + for (j = 0; j < N; j++) + { + model->a[i][j] /= ss; + } + } + for (i = 0; i < N; i++) + model->p0[i] /= s; + + free(global_mean); + + return model; +} + +void hmm_close(model_t* model) +{ + int i; + + for (i = 0; i < model->N; i++) + { + free(model->a[i]); + free(model->mu[i]); + } + free(model->a); + free(model->mu); + for (i = 0; i < model->L; i++) + free(model->cov[i]); + free(model->cov); + free(model); +} + +void hmm_train(double** x, int T, model_t* model) +{ + int i, t; + double loglik; /* overall log-likelihood at each iteration */ + + int N = model->N; + int L = model->L; + double* p0 = model->p0; + double** a = model->a; + double** mu = model->mu; + double** cov = model->cov; + + /* allocate memory */ + double** gamma = (double**) malloc(T*sizeof(double*)); + double*** xi = (double***) malloc(T*sizeof(double**)); + for (t = 0; t < T; t++) + { + gamma[t] = (double*) malloc(N*sizeof(double)); + xi[t] = (double**) malloc(N*sizeof(double*)); + for (i = 0; i < N; i++) + xi[t][i] = (double*) malloc(N*sizeof(double)); + } + + /* temporary memory */ + double* gauss_y = (double*) malloc(L*sizeof(double)); + double* gauss_z = (double*) malloc(L*sizeof(double)); + + /* obs probs P(j|{x}) */ + double** b = (double**) malloc(T*sizeof(double*)); + for (t = 0; t < T; t++) + b[t] = (double*) malloc(N*sizeof(double)); + + /* inverse covariance and its determinant */ + double** icov = (double**) malloc(L*sizeof(double*)); + for (i = 0; i < L; i++) + icov[i] = (double*) malloc(L*sizeof(double)); + double detcov; + + double thresh = 0.0001; + int niter = 50; + int iter = 0; + double loglik1, loglik2; + while(iter < niter && !(iter > 1 && (loglik - loglik1) < thresh * (loglik1 - loglik2))) + { + ++iter; + + fprintf(stderr, "calculating obsprobs...\n"); + fflush(stderr); + + /* precalculate obs probs */ + invert(cov, L, icov, &detcov); + + for (t = 0; t < T; t++) + { + //int allzero = 1; + for (i = 0; i < N; i++) + { + b[t][i] = exp(loggauss(x[t], L, mu[i], icov, detcov, gauss_y, gauss_z)); + + //if (b[t][i] != 0) + // allzero = 0; + } + /* + if (allzero) + { + printf("all the b[t][i] were zero for t = %d, correcting...\n", t); + for (i = 0; i < N; i++) + { + b[t][i] = 0.00001; + } + } + */ + } + + fprintf(stderr, "forwards-backwards...\n"); + fflush(stderr); + + forward_backwards(xi, gamma, &loglik, &loglik1, &loglik2, iter, N, T, p0, a, b); + + fprintf(stderr, "iteration %d: loglik = %f\n", iter, loglik); + fprintf(stderr, "re-estimation...\n"); + fflush(stderr); + + baum_welch(p0, a, mu, cov, N, T, L, x, xi, gamma); + + /* + printf("a:\n"); + for (i = 0; i < model->N; i++) + { + for (j = 0; j < model->N; j++) + printf("%f ", model->a[i][j]); + printf("\n"); + } + printf("\n\n"); + */ + //hmm_print(model); + } + + /* deallocate memory */ + for (t = 0; t < T; t++) + { + free(gamma[t]); + free(b[t]); + for (i = 0; i < N; i++) + free(xi[t][i]); + free(xi[t]); + } + free(gamma); + free(xi); + free(b); + + for (i = 0; i < L; i++) + free(icov[i]); + free(icov); + + free(gauss_y); + free(gauss_z); +} + +void mlss_reestimate(double* p0, double** a, double** mu, double** cov, int N, int T, int L, int* q, double** x) +{ + /* fit a single Gaussian to observations in each state */ + + /* calculate the mean observation in each state */ + + /* calculate the overall covariance */ + + /* count transitions */ + + /* estimate initial probs from transitions (???) */ +} + +void baum_welch(double* p0, double** a, double** mu, double** cov, int N, int T, int L, double** x, double*** xi, double** gamma) +{ + int i, j, t; + + double* sum_gamma = (double*) malloc(N*sizeof(double)); + + /* temporary memory */ + double* u = (double*) malloc(L*L*sizeof(double)); + double* yy = (double*) malloc(T*L*sizeof(double)); + double* yy2 = (double*) malloc(T*L*sizeof(double)); + + /* re-estimate transition probs */ + for (i = 0; i < N; i++) + { + sum_gamma[i] = 0; + for (t = 0; t < T-1; t++) + sum_gamma[i] += gamma[t][i]; + } + + for (i = 0; i < N; i++) + { + if (sum_gamma[i] == 0) + { + fprintf(stderr, "sum_gamma[%d] was zero...\n", i); + } + //double s = 0; + for (j = 0; j < N; j++) + { + a[i][j] = 0; + for (t = 0; t < T-1; t++) + a[i][j] += xi[t][i][j]; + //s += a[i][j]; + a[i][j] /= sum_gamma[i]; + } + /* + for (j = 0; j < N; j++) + { + a[i][j] /= s; + } + */ + } + + /* NB: now we need to sum gamma over all t */ + for (i = 0; i < N; i++) + sum_gamma[i] += gamma[T-1][i]; + + /* re-estimate initial probs */ + for (i = 0; i < N; i++) + p0[i] = gamma[0][i]; + + /* re-estimate covariance */ + int d, e; + double sum_sum_gamma = 0; + for (i = 0; i < N; i++) + sum_sum_gamma += sum_gamma[i]; + + /* + for (d = 0; d < L; d++) + { + for (e = d; e < L; e++) + { + cov[d][e] = 0; + for (t = 0; t < T; t++) + for (j = 0; j < N; j++) + cov[d][e] += gamma[t][j] * (x[t][d] - mu[j][d]) * (x[t][e] - mu[j][e]); + + cov[d][e] /= sum_sum_gamma; + + if (isnan(cov[d][e])) + { + printf("cov[%d][%d] was nan\n", d, e); + for (j = 0; j < N; j++) + for (i = 0; i < L; i++) + if (isnan(mu[j][i])) + printf("mu[%d][%d] was nan\n", j, i); + for (t = 0; t < T; t++) + for (j = 0; j < N; j++) + if (isnan(gamma[t][j])) + printf("gamma[%d][%d] was nan\n", t, j); + exit(-1); + } + } + } + for (d = 0; d < L; d++) + for (e = 0; e < d; e++) + cov[d][e] = cov[e][d]; + */ + + /* using BLAS */ + for (d = 0; d < L; d++) + for (e = 0; e < L; e++) + cov[d][e] = 0; + + for (j = 0; j < N; j++) + { + for (d = 0; d < L; d++) + for (t = 0; t < T; t++) + { + yy[d*T+t] = x[t][d] - mu[j][d]; + yy2[d*T+t] = gamma[t][j] * (x[t][d] - mu[j][d]); + } + + cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, L, L, T, 1.0, yy, T, yy2, T, 0, u, L); + + for (e = 0; e < L; e++) + for (d = 0; d < L; d++) + cov[d][e] += u[e*L+d]; + } + + for (d = 0; d < L; d++) + for (e = 0; e < L; e++) + cov[d][e] /= T; /* sum_sum_gamma; */ + + //printf("sum_sum_gamma = %f\n", sum_sum_gamma); /* fine, = T IS THIS ALWAYS TRUE with pooled cov?? */ + + /* re-estimate means */ + for (j = 0; j < N; j++) + { + for (d = 0; d < L; d++) + { + mu[j][d] = 0; + for (t = 0; t < T; t++) + mu[j][d] += gamma[t][j] * x[t][d]; + mu[j][d] /= sum_gamma[j]; + } + } + + /* deallocate memory */ + free(sum_gamma); + free(yy); + free(yy2); + free(u); +} + +void forward_backwards(double*** xi, double** gamma, double* loglik, double* loglik1, double* loglik2, int iter, int N, int T, double* p0, double** a, double** b) +{ + /* forwards-backwards with scaling */ + int i, j, t; + + double** alpha = (double**) malloc(T*sizeof(double*)); + double** beta = (double**) malloc(T*sizeof(double*)); + for (t = 0; t < T; t++) + { + alpha[t] = (double*) malloc(N*sizeof(double)); + beta[t] = (double*) malloc(N*sizeof(double)); + } + + /* scaling coefficients */ + double* c = (double*) malloc(T*sizeof(double)); + + /* calculate forward probs and scale coefficients */ + c[0] = 0; + for (i = 0; i < N; i++) + { + alpha[0][i] = p0[i] * b[0][i]; + c[0] += alpha[0][i]; + + //printf("p0[%d] = %f, b[0][%d] = %f\n", i, p0[i], i, b[0][i]); + } + c[0] = 1 / c[0]; + for (i = 0; i < N; i++) + { + alpha[0][i] *= c[0]; + + //printf("alpha[0][%d] = %f\n", i, alpha[0][i]); /* OK agrees with Matlab */ + } + + *loglik1 = *loglik; + *loglik = -log(c[0]); + if (iter == 2) + *loglik2 = *loglik; + + for (t = 1; t < T; t++) + { + c[t] = 0; + for (j = 0; j < N; j++) + { + alpha[t][j] = 0; + for (i = 0; i < N; i++) + alpha[t][j] += alpha[t-1][i] * a[i][j]; + alpha[t][j] *= b[t][j]; + + c[t] += alpha[t][j]; + } + + /* + if (c[t] == 0) + { + printf("c[%d] = 0, going to blow up so exiting\n", t); + for (i = 0; i < N; i++) + if (b[t][i] == 0) + fprintf(stderr, "b[%d][%d] was zero\n", t, i); + fprintf(stderr, "x[t] was \n"); + for (i = 0; i < L; i++) + fprintf(stderr, "%f ", x[t][i]); + fprintf(stderr, "\n\n"); + exit(-1); + } + */ + + c[t] = 1 / c[t]; + for (j = 0; j < N; j++) + alpha[t][j] *= c[t]; + + //printf("c[%d] = %e\n", t, c[t]); + + *loglik -= log(c[t]); + } + + /* calculate backwards probs using same coefficients */ + for (i = 0; i < N; i++) + beta[T-1][i] = 1; + t = T - 1; + while (1) + { + for (i = 0; i < N; i++) + beta[t][i] *= c[t]; + + if (t == 0) + break; + + for (i = 0; i < N; i++) + { + beta[t-1][i] = 0; + for (j = 0; j < N; j++) + beta[t-1][i] += a[i][j] * b[t][j] * beta[t][j]; + } + + t--; + } + + /* + printf("alpha:\n"); + for (t = 0; t < T; t++) + { + for (i = 0; i < N; i++) + printf("%4.4e\t\t", alpha[t][i]); + printf("\n"); + } + printf("\n\n");printf("beta:\n"); + for (t = 0; t < T; t++) + { + for (i = 0; i < N; i++) + printf("%4.4e\t\t", beta[t][i]); + printf("\n"); + } + printf("\n\n"); + */ + + /* calculate posterior probs */ + double tot; + for (t = 0; t < T; t++) + { + tot = 0; + for (i = 0; i < N; i++) + { + gamma[t][i] = alpha[t][i] * beta[t][i]; + tot += gamma[t][i]; + } + for (i = 0; i < N; i++) + { + gamma[t][i] /= tot; + + //printf("gamma[%d][%d] = %f\n", t, i, gamma[t][i]); + } + } + + for (t = 0; t < T-1; t++) + { + tot = 0; + for (i = 0; i < N; i++) + { + for (j = 0; j < N; j++) + { + xi[t][i][j] = alpha[t][i] * a[i][j] * b[t+1][j] * beta[t+1][j]; + tot += xi[t][i][j]; + } + } + for (i = 0; i < N; i++) + for (j = 0; j < N; j++) + xi[t][i][j] /= tot; + } + + /* + // CHECK - fine + // gamma[t][i] = \sum_j{xi[t][i][j]} + tot = 0; + for (j = 0; j < N; j++) + tot += xi[3][1][j]; + printf("gamma[3][1] = %f, sum_j(xi[3][1][j]) = %f\n", gamma[3][1], tot); + */ + + for (t = 0; t < T; t++) + { + free(alpha[t]); + free(beta[t]); + } + free(alpha); + free(beta); + free(c); +} + +void viterbi_decode(double** x, int T, model_t* model, int* q) +{ + int i, j, t; + double max; + + int N = model->N; + int L = model->L; + double* p0 = model->p0; + double** a = model->a; + double** mu = model->mu; + double** cov = model->cov; + + /* inverse covariance and its determinant */ + double** icov = (double**) malloc(L*sizeof(double*)); + for (i = 0; i < L; i++) + icov[i] = (double*) malloc(L*sizeof(double)); + double detcov; + + double** logb = (double**) malloc(T*sizeof(double*)); + double** phi = (double**) malloc(T*sizeof(double*)); + int** psi = (int**) malloc(T*sizeof(int*)); + for (t = 0; t < T; t++) + { + logb[t] = (double*) malloc(N*sizeof(double)); + phi[t] = (double*) malloc(N*sizeof(double)); + psi[t] = (int*) malloc(N*sizeof(int)); + } + + /* temporary memory */ + double* gauss_y = (double*) malloc(L*sizeof(double)); + double* gauss_z = (double*) malloc(L*sizeof(double)); + + /* calculate observation logprobs */ + invert(cov, L, icov, &detcov); + for (t = 0; t < T; t++) + for (i = 0; i < N; i++) + logb[t][i] = loggauss(x[t], L, mu[i], icov, detcov, gauss_y, gauss_z); + + /* initialise */ + for (i = 0; i < N; i++) + { + phi[0][i] = log(p0[i]) + logb[0][i]; + psi[0][i] = 0; + } + + for (t = 1; t < T; t++) + { + for (j = 0; j < N; j++) + { + max = -1000000; // TODO: what should this be?? = smallest possible sumlogprob + psi[t][j] = 0; + for (i = 0; i < N; i++) + { + if (phi[t-1][i] + log(a[i][j]) > max) + { + max = phi[t-1][i] + log(a[i][j]); + phi[t][j] = max + logb[t][j]; + psi[t][j] = i; + } + } + } + } + + /* find maximising state at time T-1 */ + max = phi[T-1][0]; + q[T-1] = 0; + for (i = 1; i < N; i++) + { + if (phi[T-1][i] > max) + { + max = phi[T-1][i]; + q[T-1] = i; + } + } + + + /* track back */ + t = T - 2; + while (t >= 0) + { + q[t] = psi[t+1][q[t+1]]; + t--; + } + + /* de-allocate memory */ + for (i = 0; i < L; i++) + free(icov[i]); + free(icov); + for (t = 0; t < T; t++) + { + free(logb[t]); + free(phi[t]); + free(psi[t]); + } + free(logb); + free(phi); + free(psi); + + free(gauss_y); + free(gauss_z); +} + +/* invert matrix and calculate determinant using LAPACK */ +void invert(double** cov, int L, double** icov, double* detcov) +{ + /* copy square matrix into a vector in column-major order */ + double* a = (double*) malloc(L*L*sizeof(double)); + int i, j; + for(j=0; j < L; j++) + for (i=0; i < L; i++) + a[j*L+i] = cov[i][j]; + + long M = (long) L; + long* ipiv = (long*) malloc(L*L*sizeof(int)); + long ret; + + /* LU decomposition */ + ret = dgetrf_(&M, &M, a, &M, ipiv, &ret); /* ret should be zero, negative if cov is singular */ + if (ret < 0) + { + fprintf(stderr, "Covariance matrix was singular, couldn't invert\n"); + exit(-1); + } + + /* find determinant */ + double det = 1; + for(i = 0; i < L; i++) + det *= a[i*L+i]; + // TODO: get this to work!!! If detcov < 0 then cov is bad anyway... + /* + int sign = 1; + for (i = 0; i < L; i++) + if (ipiv[i] != i) + sign = -sign; + det *= sign; + */ + if (det < 0) + det = -det; + *detcov = det; + + /* allocate required working storage */ + long lwork = -1; + double lwbest; + dgetri_(&M, a, &M, ipiv, &lwbest, &lwork, &ret); + lwork = (long) lwbest; + double* work = (double*) malloc(lwork*sizeof(double)); + + /* find inverse */ + dgetri_(&M, a, &M, ipiv, work, &lwork, &ret); + + for(j=0; j < L; j++) + for (i=0; i < L; i++) + icov[i][j] = a[j*L+i]; + + free(work); + free(a); +} + +/* probability of multivariate Gaussian given mean, inverse and determinant of covariance */ +double gauss(double* x, int L, double* mu, double** icov, double detcov, double* y, double* z) +{ + int i, j; + double s = 0; + for (i = 0; i < L; i++) + y[i] = x[i] - mu[i]; + for (i = 0; i < L; i++) + { + //z[i] = 0; + //for (j = 0; j < L; j++) + // z[i] += icov[i][j] * y[j]; + z[i] = cblas_ddot(L, &icov[i][0], 1, y, 1); + } + s = cblas_ddot(L, z, 1, y, 1); + //for (i = 0; i < L; i++) + // s += z[i] * y[i]; + + return exp(-s/2.0) / (pow(2*PI, L/2.0) * sqrt(detcov)); +} + +/* log probability of multivariate Gaussian given mean, inverse and determinant of covariance */ +double loggauss(double* x, int L, double* mu, double** icov, double detcov, double* y, double* z) +{ + int i, j; + double s = 0; + double ret; + for (i = 0; i < L; i++) + y[i] = x[i] - mu[i]; + for (i = 0; i < L; i++) + { + //z[i] = 0; + //for (j = 0; j < L; j++) + // z[i] += icov[i][j] * y[j]; + z[i] = cblas_ddot(L, &icov[i][0], 1, y, 1); + } + s = cblas_ddot(L, z, 1, y, 1); + //for (i = 0; i < L; i++) + // s += z[i] * y[i]; + + ret = -0.5 * (s + L * log(2*PI) + log(detcov)); + + /* + // TEST + if (isinf(ret) > 0) + printf("loggauss returning infinity\n"); + if (isinf(ret) < 0) + printf("loggauss returning -infinity\n"); + if (isnan(ret)) + printf("loggauss returning nan\n"); + */ + + return ret; +} + +void hmm_print(model_t* model) +{ + int i, j; + printf("p0:\n"); + for (i = 0; i < model->N; i++) + printf("%f ", model->p0[i]); + printf("\n\n"); + printf("a:\n"); + for (i = 0; i < model->N; i++) + { + for (j = 0; j < model->N; j++) + printf("%f ", model->a[i][j]); + printf("\n"); + } + printf("\n\n"); + printf("mu:\n"); + for (i = 0; i < model->N; i++) + { + for (j = 0; j < model->L; j++) + printf("%f ", model->mu[i][j]); + printf("\n"); + } + printf("\n\n"); + printf("cov:\n"); + for (i = 0; i < model->L; i++) + { + for (j = 0; j < model->L; j++) + printf("%f ", model->cov[i][j]); + printf("\n"); + } + printf("\n\n"); +} + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/hmm/hmm.h Wed Jan 09 10:31:29 2008 +0000 @@ -0,0 +1,39 @@ +#ifndef _HMM_H +#define _HMM_H + +/* + * hmm.h + * soundbite + * + * Created by Mark Levy on 12/02/2006. + * Copyright 2006 Centre for Digital Music, Queen Mary, University of London. All rights reserved. + * + */ + +#ifndef PI +#define PI 3.14159265358979323846264338327950288 +#endif + +typedef struct _model_t { + int N; /* number of states */ + double* p0; /* initial probs */ + double** a; /* transition probs */ + int L; /* dimensionality of data */ + double** mu; /* state means */ + double** cov; /* covariance, tied between all states */ +} model_t; + +void hmm_train(double** x, int T, model_t* model); /* with scaling */ +void forward_backwards(double*** xi, double** gamma, double* loglik, double* loglik1, double* loglik2, int iter, + int N, int T, double* p0, double** a, double** b); +void baum_welch(double* p0, double** a, double** mu, double** cov, int N, int T, int L, double** x, double*** xi, double** gamma); +void viterbi_decode(double** x, int T, model_t* model, int* q); /* using logs */ +model_t* hmm_init(double** x, int T, int L, int N); +void hmm_close(model_t* model); +void invert(double** cov, int L, double** icov, double* detcov); /* uses LAPACK (included with Mac OSX) */ +double gauss(double* x, int L, double* mu, double** icov, double detcov, double* y, double* z); +double loggauss(double* x, int L, double* mu, double** icov, double detcov, double* y, double* z); +void hmm_print(model_t* model); + +#endif +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/maths/Correlation.cpp Wed Jan 09 10:31:29 2008 +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/maths/Correlation.h Wed Jan 09 10:31:29 2008 +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/maths/Histogram.h Wed Jan 09 10:31:29 2008 +0000 @@ -0,0 +1,393 @@ +/* -*- 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/maths/MathAliases.h Wed Jan 09 10:31:29 2008 +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/maths/MathUtilities.cpp Wed Jan 09 10:31:29 2008 +0000 @@ -0,0 +1,221 @@ +/* -*- 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 ; + } + + } +} + +int MathUtilities::getMax( double* pData, unsigned int Length, double* pMax ) +{ + unsigned int index = 0; + unsigned int i; + double temp = 0.0; + + double max = pData[0]; + + for( i = 0; i < Length; i++) + { + temp = pData[ i ]; + + if( temp > max ) + { + max = temp ; + index = i; + } + + } + + *pMax = max; + + + return index; +} + +void MathUtilities::circShift( double* pData, int length, int shift) +{ + shift = shift % length; + double temp; + int i,n; + + for( i = 0; i < shift; i++) + { + temp=*(pData + length - 1); + + for( n = length-2; n >= 0; n--) + { + *(pData+n+1)=*(pData+n); + } + + *pData = temp; + } +} + +int MathUtilities::compareInt (const void * a, const void * b) +{ + return ( *(int*)a - *(int*)b ); +} +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/maths/MathUtilities.h Wed Jan 09 10:31:29 2008 +0000 @@ -0,0 +1,33 @@ +/* -*- 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 ); + static void circShift( double* data, int length, int shift); + static int getMax( double* data, unsigned int length, double* max ); + static int compareInt(const void * a, const void * b); +}; + +#endif
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/maths/Polyfit.h Wed Jan 09 10:31:29 2008 +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/maths/pca/data.txt Wed Jan 09 10:31:29 2008 +0000 @@ -0,0 +1,36 @@ + 3.0 3.0 3.0 3.0 3.0 3.0 35.0 45.0 + 53.0 55.0 58.0 113.0 113.0 86.0 67.0 90.0 + 3.5 3.5 4.0 4.0 4.5 4.5 46.0 59.0 + 63.0 58.0 58.0 125.0 126.0 110.0 78.0 97.0 + 4.0 4.0 4.5 4.5 5.0 5.0 48.0 60.0 + 68.0 65.0 65.0 123.0 123.0 117.0 87.0 108.0 + 5.0 5.0 5.0 5.5 5.5 5.5 46.0 63.0 + 70.0 64.0 63.0 116.0 119.0 115.0 97.0 112.0 + 6.0 6.0 6.0 6.0 6.5 6.5 51.0 69.0 + 77.0 70.0 71.0 120.0 122.0 122.0 96.0 123.0 + 11.0 11.0 11.0 11.0 11.0 11.0 64.0 75.0 + 81.0 79.0 79.0 112.0 114.0 113.0 98.0 115.0 + 20.0 20.0 20.0 20.0 20.0 20.0 76.0 86.0 + 93.0 92.0 91.0 104.0 104.5 107.0 97.5 104.0 + 30.0 30.0 30.0 30.0 30.1 30.2 84.0 96.0 + 98.0 99.0 96.0 101.0 102.0 99.0 94.0 99.0 + 30.0 33.4 36.8 40.0 43.0 45.6 100.0 106.0 + 106.0 108.0 101.0 99.0 98.0 99.0 95.0 95.0 + 42.0 44.0 46.0 48.0 50.0 51.0 109.0 111.0 + 110.0 110.0 103.0 95.5 95.5 95.0 92.5 92.0 + 60.0 61.7 63.5 65.5 67.3 69.2 122.0 124.0 + 124.0 121.0 103.0 93.2 92.5 92.2 90.0 90.8 + 70.0 70.1 70.2 70.3 70.4 70.5 137.0 132.0 + 134.0 128.0 101.0 91.7 90.2 88.8 87.3 85.8 + 78.0 77.6 77.2 76.8 76.4 76.0 167.0 159.0 + 152.0 144.0 103.0 89.8 87.7 85.7 83.7 81.8 + 98.9 97.8 96.7 95.5 94.3 93.2 183.0 172.0 + 162.0 152.0 102.0 87.5 85.3 83.3 81.3 79.3 + 160.0 157.0 155.0 152.0 149.0 147.0 186.0 175.0 + 165.0 156.0 120.0 87.0 84.9 82.8 80.8 79.0 + 272.0 266.0 260.0 254.0 248.0 242.0 192.0 182.0 + 170.0 159.0 131.0 88.0 85.8 83.7 81.6 79.6 + 382.0 372.0 362.0 352.0 343.0 333.0 205.0 192.0 + 178.0 166.0 138.0 86.2 84.0 82.0 79.8 77.5 + 770.0 740.0 710.0 680.0 650.0 618.0 226.0 207.0 + 195.0 180.0 160.0 82.9 80.2 77.7 75.2 72.7 \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/maths/pca/pca.c Wed Jan 09 10:31:29 2008 +0000 @@ -0,0 +1,356 @@ +/*********************************/ +/* Principal Components Analysis */ +/*********************************/ + +/*********************************************************************/ +/* Principal Components Analysis or the Karhunen-Loeve expansion is a + classical method for dimensionality reduction or exploratory data + analysis. One reference among many is: F. Murtagh and A. Heck, + Multivariate Data Analysis, Kluwer Academic, Dordrecht, 1987. + + Author: + F. Murtagh + Phone: + 49 89 32006298 (work) + + 49 89 965307 (home) + Earn/Bitnet: fionn@dgaeso51, fim@dgaipp1s, murtagh@stsci + Span: esomc1::fionn + Internet: murtagh@scivax.stsci.edu + + F. Murtagh, Munich, 6 June 1989 */ +/*********************************************************************/ + +#include <stdio.h> +#include <stdlib.h> +#include <math.h> + +#include "pca.h" + +#define SIGN(a, b) ( (b) < 0 ? -fabs(a) : fabs(a) ) + +/** Variance-covariance matrix: creation *****************************/ + +/* Create m * m covariance matrix from given n * m data matrix. */ +void covcol(double** data, int n, int m, double** symmat) +{ +double *mean; +int i, j, j1, j2; + +/* Allocate storage for mean vector */ + +mean = (double*) malloc(m*sizeof(double)); + +/* Determine mean of column vectors of input data matrix */ + +for (j = 0; j < m; j++) + { + mean[j] = 0.0; + for (i = 0; i < n; i++) + { + mean[j] += data[i][j]; + } + mean[j] /= (double)n; + } + +/* +printf("\nMeans of column vectors:\n"); +for (j = 0; j < m; j++) { + printf("%12.1f",mean[j]); } printf("\n"); + */ + +/* Center the column vectors. */ + +for (i = 0; i < n; i++) + { + for (j = 0; j < m; j++) + { + data[i][j] -= mean[j]; + } + } + +/* Calculate the m * m covariance matrix. */ +for (j1 = 0; j1 < m; j1++) + { + for (j2 = j1; j2 < m; j2++) + { + symmat[j1][j2] = 0.0; + for (i = 0; i < n; i++) + { + symmat[j1][j2] += data[i][j1] * data[i][j2]; + } + symmat[j2][j1] = symmat[j1][j2]; + } + } + +free(mean); + +return; + +} + +/** Error handler **************************************************/ + +void erhand(char* err_msg) +{ + fprintf(stderr,"Run-time error:\n"); + fprintf(stderr,"%s\n", err_msg); + fprintf(stderr,"Exiting to system.\n"); + exit(1); +} + + +/** Reduce a real, symmetric matrix to a symmetric, tridiag. matrix. */ + +/* Householder reduction of matrix a to tridiagonal form. +Algorithm: Martin et al., Num. Math. 11, 181-195, 1968. +Ref: Smith et al., Matrix Eigensystem Routines -- EISPACK Guide +Springer-Verlag, 1976, pp. 489-494. +W H Press et al., Numerical Recipes in C, Cambridge U P, +1988, pp. 373-374. */ +void tred2(double** a, int n, double* d, double* e) +{ + int l, k, j, i; + double scale, hh, h, g, f; + + for (i = n-1; i >= 1; i--) + { + l = i - 1; + h = scale = 0.0; + if (l > 0) + { + for (k = 0; k <= l; k++) + scale += fabs(a[i][k]); + if (scale == 0.0) + e[i] = a[i][l]; + else + { + for (k = 0; k <= l; k++) + { + a[i][k] /= scale; + h += a[i][k] * a[i][k]; + } + f = a[i][l]; + g = f>0 ? -sqrt(h) : sqrt(h); + e[i] = scale * g; + h -= f * g; + a[i][l] = f - g; + f = 0.0; + for (j = 0; j <= l; j++) + { + a[j][i] = a[i][j]/h; + g = 0.0; + for (k = 0; k <= j; k++) + g += a[j][k] * a[i][k]; + for (k = j+1; k <= l; k++) + g += a[k][j] * a[i][k]; + e[j] = g / h; + f += e[j] * a[i][j]; + } + hh = f / (h + h); + for (j = 0; j <= l; j++) + { + f = a[i][j]; + e[j] = g = e[j] - hh * f; + for (k = 0; k <= j; k++) + a[j][k] -= (f * e[k] + g * a[i][k]); + } + } + } + else + e[i] = a[i][l]; + d[i] = h; + } + d[0] = 0.0; + e[0] = 0.0; + for (i = 0; i < n; i++) + { + l = i - 1; + if (d[i]) + { + for (j = 0; j <= l; j++) + { + g = 0.0; + for (k = 0; k <= l; k++) + g += a[i][k] * a[k][j]; + for (k = 0; k <= l; k++) + a[k][j] -= g * a[k][i]; + } + } + d[i] = a[i][i]; + a[i][i] = 1.0; + for (j = 0; j <= l; j++) + a[j][i] = a[i][j] = 0.0; + } +} + +/** Tridiagonal QL algorithm -- Implicit **********************/ + +void tqli(double* d, double* e, int n, double** z) +{ + int m, l, iter, i, k; + double s, r, p, g, f, dd, c, b; + + for (i = 1; i < n; i++) + e[i-1] = e[i]; + e[n-1] = 0.0; + for (l = 0; l < n; l++) + { + iter = 0; + do + { + for (m = l; m < n-1; m++) + { + dd = fabs(d[m]) + fabs(d[m+1]); + if (fabs(e[m]) + dd == dd) break; + } + if (m != l) + { + if (iter++ == 30) erhand("No convergence in TLQI."); + g = (d[l+1] - d[l]) / (2.0 * e[l]); + r = sqrt((g * g) + 1.0); + g = d[m] - d[l] + e[l] / (g + SIGN(r, g)); + s = c = 1.0; + p = 0.0; + for (i = m-1; i >= l; i--) + { + f = s * e[i]; + b = c * e[i]; + if (fabs(f) >= fabs(g)) + { + c = g / f; + r = sqrt((c * c) + 1.0); + e[i+1] = f * r; + c *= (s = 1.0/r); + } + else + { + s = f / g; + r = sqrt((s * s) + 1.0); + e[i+1] = g * r; + s *= (c = 1.0/r); + } + g = d[i+1] - p; + r = (d[i] - g) * s + 2.0 * c * b; + p = s * r; + d[i+1] = g + p; + g = c * r - b; + for (k = 0; k < n; k++) + { + f = z[k][i+1]; + z[k][i+1] = s * z[k][i] + c * f; + z[k][i] = c * z[k][i] - s * f; + } + } + d[l] = d[l] - p; + e[l] = g; + e[m] = 0.0; + } + } while (m != l); + } +} + +/* In place projection onto basis vectors */ +void pca_project(double** data, int n, int m, int ncomponents) +{ + int i, j, k, k2; + double **symmat, **symmat2, *evals, *interm; + + //TODO: assert ncomponents < m + + symmat = (double**) malloc(m*sizeof(double*)); + for (i = 0; i < m; i++) + symmat[i] = (double*) malloc(m*sizeof(double)); + + covcol(data, n, m, symmat); + + /********************************************************************* + Eigen-reduction + **********************************************************************/ + + /* Allocate storage for dummy and new vectors. */ + evals = (double*) malloc(m*sizeof(double)); /* Storage alloc. for vector of eigenvalues */ + interm = (double*) malloc(m*sizeof(double)); /* Storage alloc. for 'intermediate' vector */ + //MALLOC_ARRAY(symmat2,m,m,double); + //for (i = 0; i < m; i++) { + // for (j = 0; j < m; j++) { + // symmat2[i][j] = symmat[i][j]; /* Needed below for col. projections */ + // } + //} + tred2(symmat, m, evals, interm); /* Triangular decomposition */ +tqli(evals, interm, m, symmat); /* Reduction of sym. trid. matrix */ +/* evals now contains the eigenvalues, +columns of symmat now contain the associated eigenvectors. */ + +/* + printf("\nEigenvalues:\n"); + for (j = m-1; j >= 0; j--) { + printf("%18.5f\n", evals[j]); } + printf("\n(Eigenvalues should be strictly positive; limited\n"); + printf("precision machine arithmetic may affect this.\n"); + printf("Eigenvalues are often expressed as cumulative\n"); + printf("percentages, representing the 'percentage variance\n"); + printf("explained' by the associated axis or principal component.)\n"); + + printf("\nEigenvectors:\n"); + printf("(First three; their definition in terms of original vbes.)\n"); + for (j = 0; j < m; j++) { + for (i = 1; i <= 3; i++) { + printf("%12.4f", symmat[j][m-i]); } + printf("\n"); } + */ + +/* Form projections of row-points on prin. components. */ +/* Store in 'data', overwriting original data. */ +for (i = 0; i < n; i++) { + for (j = 0; j < m; j++) { + interm[j] = data[i][j]; } /* data[i][j] will be overwritten */ + for (k = 0; k < ncomponents; k++) { + data[i][k] = 0.0; + for (k2 = 0; k2 < m; k2++) { + data[i][k] += interm[k2] * symmat[k2][m-k-1]; } + } +} + +/* +printf("\nProjections of row-points on first 3 prin. comps.:\n"); + for (i = 0; i < n; i++) { + for (j = 0; j < 3; j++) { + printf("%12.4f", data[i][j]); } + printf("\n"); } + */ + +/* Form projections of col.-points on first three prin. components. */ +/* Store in 'symmat2', overwriting what was stored in this. */ +//for (j = 0; j < m; j++) { +// for (k = 0; k < m; k++) { +// interm[k] = symmat2[j][k]; } /*symmat2[j][k] will be overwritten*/ +// for (i = 0; i < 3; i++) { +// symmat2[j][i] = 0.0; +// for (k2 = 0; k2 < m; k2++) { +// symmat2[j][i] += interm[k2] * symmat[k2][m-i-1]; } +// if (evals[m-i-1] > 0.0005) /* Guard against zero eigenvalue */ +// symmat2[j][i] /= sqrt(evals[m-i-1]); /* Rescale */ +// else +// symmat2[j][i] = 0.0; /* Standard kludge */ +// } +// } + +/* + printf("\nProjections of column-points on first 3 prin. comps.:\n"); + for (j = 0; j < m; j++) { + for (k = 0; k < 3; k++) { + printf("%12.4f", symmat2[j][k]); } + printf("\n"); } + */ + + +for (i = 0; i < m; i++) + free(symmat[i]); +free(symmat); +//FREE_ARRAY(symmat2,m); +free(evals); +free(interm); + +} + + +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/maths/pca/pca.h Wed Jan 09 10:31:29 2008 +0000 @@ -0,0 +1,16 @@ +#ifndef _PCA_H +#define _PCA_H + +/* + * pca.h + * soundbite + * + * Created by Mark Levy on 08/02/2006. + * Copyright 2006 Centre for Digital Music, Queen Mary, University of London. All rights reserved. + * + */ + +void pca_project(double** data, int n, int m, int ncomponents); + + +#endif \ No newline at end of file
--- a/qm-dsp.pro Wed Nov 21 16:53:51 2007 +0000 +++ b/qm-dsp.pro Wed Jan 09 10:31:29 2008 +0000 @@ -24,11 +24,6 @@ dsp/chromagram/ChromaProcess.h \ dsp/chromagram/ConstantQ.h \ dsp/keydetection/GetKeyMode.h \ - dsp/maths/Correlation.h \ - dsp/maths/Histogram.h \ - dsp/maths/MathAliases.h \ - dsp/maths/MathUtilities.h \ - dsp/maths/Polyfit.h \ dsp/onsets/DetectionFunction.h \ dsp/onsets/PeakPicking.h \ dsp/phasevocoder/PhaseVocoder.h \ @@ -41,14 +36,17 @@ dsp/tonal/ChangeDetectionFunction.h \ dsp/tonal/TCSgram.h \ dsp/tonal/TonalEstimator.h \ - dsp/transforms/FFT.h + dsp/transforms/FFT.h \ + maths/Correlation.h \ + maths/Histogram.h \ + maths/MathAliases.h \ + maths/MathUtilities.h \ + maths/Polyfit.h SOURCES += base/Pitch.cpp \ dsp/chromagram/Chromagram.cpp \ dsp/chromagram/ChromaProcess.cpp \ dsp/chromagram/ConstantQ.cpp \ dsp/keydetection/GetKeyMode.cpp \ - dsp/maths/Correlation.cpp \ - dsp/maths/MathUtilities.cpp \ dsp/onsets/DetectionFunction.cpp \ dsp/onsets/PeakPicking.cpp \ dsp/phasevocoder/PhaseVocoder.cpp \ @@ -61,4 +59,6 @@ dsp/tonal/ChangeDetectionFunction.cpp \ dsp/tonal/TCSgram.cpp \ dsp/tonal/TonalEstimator.cpp \ - dsp/transforms/FFT.cpp + dsp/transforms/FFT.cpp \ + maths/Correlation.cpp \ + maths/MathUtilities.cpp