add queen mary DSP library

git-svn-id: svn://localhost/ardour2/branches/3.0@9029 d708f5d6-7413-0410-9779-e7cbd77b26cf

14 files changed, 1667 insertions, 0 deletions

diff --git a/libs/qm-dsp/maths/Correlation.cpp b/libs/qm-dsp/maths/Correlation.cpp
new file mode 100644
index 0000000000..17ee28f749
--- /dev/null
+++ b/libs/qm-dsp/maths/Correlation.cpp
@@ -0,0 +1,56 @@
+/* -*- 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 2005-2006 Christian Landone.
+
+    This program is free software; you can redistribute it and/or
+    modify it under the terms of the GNU General Public License as
+    published by the Free Software Foundation; either version 2 of the
+    License, or (at your option) any later version.  See the file
+    COPYING included with this distribution for more information.
+*/
+
+#include "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;
+    }
+}
diff --git a/libs/qm-dsp/maths/Correlation.h b/libs/qm-dsp/maths/Correlation.h
new file mode 100644
index 0000000000..85fcc7316b
--- /dev/null
+++ b/libs/qm-dsp/maths/Correlation.h
@@ -0,0 +1,30 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file 2005-2006 Christian Landone.
+
+    This program is free software; you can redistribute it and/or
+    modify it under the terms of the GNU General Public License as
+    published by the Free Software Foundation; either version 2 of the
+    License, or (at your option) any later version.  See the file
+    COPYING included with this distribution for more information.
+*/
+
+#ifndef 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 // 
diff --git a/libs/qm-dsp/maths/CosineDistance.cpp b/libs/qm-dsp/maths/CosineDistance.cpp
new file mode 100644
index 0000000000..13ab9ce0e8
--- /dev/null
+++ b/libs/qm-dsp/maths/CosineDistance.cpp
@@ -0,0 +1,47 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2008 Kurt Jacobson.
+
+    This program is free software; you can redistribute it and/or
+    modify it under the terms of the GNU General Public License as
+    published by the Free Software Foundation; either version 2 of the
+    License, or (at your option) any later version.  See the file
+    COPYING included with this distribution for more information.
+*/
+
+#include "CosineDistance.h"
+
+#include <iostream>
+#include <limits>
+
+using std::cerr;
+
+double CosineDistance::distance(const vector<double> &v1,
+                                const vector<double> &v2)
+{
+    dist = 1.0; dDenTot = 0; dDen1 = 0; dDen2 = 0; dSum1 =0;
+    double small = 1e-20;
+
+    //check if v1, v2 same size
+    if (v1.size() != v2.size())
+    {
+        cerr << "CosineDistance::distance: ERROR: vectors not the same size\n";
+        return 1.0;
+    }
+    else
+    {
+        for(int i=0; i<v1.size(); i++)
+        {
+            dSum1 += v1[i]*v2[i];
+            dDen1 += v1[i]*v1[i];
+            dDen2 += v2[i]*v2[i];
+        }
+        dDenTot = sqrt(fabs(dDen1*dDen2)) + small;
+        dist = 1-((dSum1)/dDenTot);
+        return dist;
+    }
+}
diff --git a/libs/qm-dsp/maths/CosineDistance.h b/libs/qm-dsp/maths/CosineDistance.h
new file mode 100644
index 0000000000..4cdfd4f521
--- /dev/null
+++ b/libs/qm-dsp/maths/CosineDistance.h
@@ -0,0 +1,37 @@
+/* -*- 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 2008 Kurt Jacobson.
+
+    This program is free software; you can redistribute it and/or
+    modify it under the terms of the GNU General Public License as
+    published by the Free Software Foundation; either version 2 of the
+    License, or (at your option) any later version.  See the file
+    COPYING included with this distribution for more information.
+*/
+
+#ifndef COSINEDISTANCE_H
+#define COSINEDISTANCE_H
+
+#include <vector>
+#include <math.h>
+
+using std::vector;
+
+class CosineDistance
+{
+public:
+    CosineDistance() { }
+    ~CosineDistance() { }
+
+    double distance(const vector<double> &v1, const vector<double> &v2);
+
+protected:
+    double dist, dDenTot, dDen1, dDen2, dSum1;
+};
+
+#endif
+
diff --git a/libs/qm-dsp/maths/KLDivergence.cpp b/libs/qm-dsp/maths/KLDivergence.cpp
new file mode 100644
index 0000000000..3c3cb134ca
--- /dev/null
+++ b/libs/qm-dsp/maths/KLDivergence.cpp
@@ -0,0 +1,64 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file copyright 2008 QMUL
+
+    This program is free software; you can redistribute it and/or
+    modify it under the terms of the GNU General Public License as
+    published by the Free Software Foundation; either version 2 of the
+    License, or (at your option) any later version.  See the file
+    COPYING included with this distribution for more information.
+*/
+
+#include "KLDivergence.h"
+
+#include <cmath>
+
+double KLDivergence::distanceGaussian(const vector<double> &m1,
+                                      const vector<double> &v1,
+                                      const vector<double> &m2,
+                                      const vector<double> &v2)
+{
+    int sz = m1.size();
+
+    double d = -2.0 * sz;
+    double small = 1e-20;
+
+    for (int k = 0; k < sz; ++k) {
+
+        double kv1 = v1[k] + small;
+        double kv2 = v2[k] + small;
+        double km = (m1[k] - m2[k]) + small;
+
+        d += kv1 / kv2 + kv2 / kv1;
+        d += km * (1.0 / kv1 + 1.0 / kv2) * km;
+    }
+
+    d /= 2.0;
+
+    return d;
+}
+
+double KLDivergence::distanceDistribution(const vector<double> &d1,
+                                          const vector<double> &d2,
+                                          bool symmetrised)
+{
+    int sz = d1.size();
+
+    double d = 0;
+    double small = 1e-20;
+    
+    for (int i = 0; i < sz; ++i) {
+        d += d1[i] * log10((d1[i] + small) / (d2[i] + small));
+    }
+
+    if (symmetrised) {
+        d += distanceDistribution(d2, d1, false);
+    }
+
+    return d;
+}
+
diff --git a/libs/qm-dsp/maths/KLDivergence.h b/libs/qm-dsp/maths/KLDivergence.h
new file mode 100644
index 0000000000..6f0342955d
--- /dev/null
+++ b/libs/qm-dsp/maths/KLDivergence.h
@@ -0,0 +1,54 @@
+/* -*- 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 2008 QMUL.
+
+    This program is free software; you can redistribute it and/or
+    modify it under the terms of the GNU General Public License as
+    published by the Free Software Foundation; either version 2 of the
+    License, or (at your option) any later version.  See the file
+    COPYING included with this distribution for more information.
+*/
+
+#ifndef KLDIVERGENCE_H
+#define KLDIVERGENCE_H
+
+#include <vector>
+
+using std::vector;
+
+/**
+ * Helper methods for calculating Kullback-Leibler divergences.
+ */
+class KLDivergence
+{
+public:
+    KLDivergence() { }
+    ~KLDivergence() { }
+
+    /**
+     * Calculate a symmetrised Kullback-Leibler divergence of Gaussian
+     * models based on mean and variance vectors.  All input vectors
+     * must be of equal size.
+     */
+    double distanceGaussian(const vector<double> &means1,
+                            const vector<double> &variances1,
+                            const vector<double> &means2,
+                            const vector<double> &variances2);
+
+    /**
+     * Calculate a Kullback-Leibler divergence of two probability
+     * distributions.  Input vectors must be of equal size.  If
+     * symmetrised is true, the result will be the symmetrised
+     * distance (equal to KL(d1, d2) + KL(d2, d1)).
+     */
+    double distanceDistribution(const vector<double> &d1,
+                                const vector<double> &d2,
+                                bool symmetrised);
+};
+
+#endif
+
diff --git a/libs/qm-dsp/maths/MathAliases.h b/libs/qm-dsp/maths/MathAliases.h
new file mode 100644
index 0000000000..8660129cbb
--- /dev/null
+++ b/libs/qm-dsp/maths/MathAliases.h
@@ -0,0 +1,60 @@
+/* -*- 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 2005-2006 Christian Landone.
+
+    This program is free software; you can redistribute it and/or
+    modify it under the terms of the GNU General Public License as
+    published by the Free Software Foundation; either version 2 of the
+    License, or (at your option) any later version.  See the file
+    COPYING included with this distribution for more information.
+*/
+
+#ifndef 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
diff --git a/libs/qm-dsp/maths/MathUtilities.cpp b/libs/qm-dsp/maths/MathUtilities.cpp
new file mode 100644
index 0000000000..809874121e
--- /dev/null
+++ b/libs/qm-dsp/maths/MathUtilities.cpp
@@ -0,0 +1,401 @@
+/* -*- 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 2005-2006 Christian Landone.
+
+    This program is free software; you can redistribute it and/or
+    modify it under the terms of the GNU General Public License as
+    published by the Free Software Foundation; either version 2 of the
+    License, or (at your option) any later version.  See the file
+    COPYING included with this distribution for more information.
+*/
+
+#include "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), double(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), double(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;
+}
+
+double MathUtilities::mean(const std::vector<double> &src,
+                           unsigned int start,
+                           unsigned int count)
+{
+    double sum = 0.;
+	
+    for (int i = 0; i < count; ++i)
+    {
+        sum += src[start + i];
+    }
+
+    return sum / count;
+}
+
+void MathUtilities::getFrameMinMax(const double *data, unsigned int len, double *min, double *max)
+{
+    unsigned int i;
+    double temp = 0.0;
+    double a=0.0;
+
+    if (len == 0) {
+        *min = *max = 0;
+        return;
+    }
+	
+    *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;
+		}
+		
+   	}
+
+	if (pMax) *pMax = max;
+
+
+	return index;
+}
+
+int MathUtilities::getMax( const std::vector<double> & data, double* pMax )
+{
+	unsigned int index = 0;
+	unsigned int i;
+	double temp = 0.0;
+	
+	double max = data[0];
+
+	for( i = 0; i < data.size(); i++)
+	{
+		temp = data[ i ];
+
+		if( temp > max )
+		{
+			max =  temp ;
+			index = i;
+		}
+		
+   	}
+
+	if (pMax) *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 );
+}
+
+void MathUtilities::normalise(double *data, int length, NormaliseType type)
+{
+    switch (type) {
+
+    case NormaliseNone: return;
+
+    case NormaliseUnitSum:
+    {
+        double sum = 0.0;
+        for (int i = 0; i < length; ++i) {
+            sum += data[i];
+        }
+        if (sum != 0.0) {
+            for (int i = 0; i < length; ++i) {
+                data[i] /= sum;
+            }
+        }
+    }
+    break;
+
+    case NormaliseUnitMax:
+    {
+        double max = 0.0;
+        for (int i = 0; i < length; ++i) {
+            if (fabs(data[i]) > max) {
+                max = fabs(data[i]);
+            }
+        }
+        if (max != 0.0) {
+            for (int i = 0; i < length; ++i) {
+                data[i] /= max;
+            }
+        }
+    }
+    break;
+
+    }
+}
+
+void MathUtilities::normalise(std::vector<double> &data, NormaliseType type)
+{
+    switch (type) {
+
+    case NormaliseNone: return;
+
+    case NormaliseUnitSum:
+    {
+        double sum = 0.0;
+        for (int i = 0; i < data.size(); ++i) sum += data[i];
+        if (sum != 0.0) {
+            for (int i = 0; i < data.size(); ++i) data[i] /= sum;
+        }
+    }
+    break;
+
+    case NormaliseUnitMax:
+    {
+        double max = 0.0;
+        for (int i = 0; i < data.size(); ++i) {
+            if (fabs(data[i]) > max) max = fabs(data[i]);
+        }
+        if (max != 0.0) {
+            for (int i = 0; i < data.size(); ++i) data[i] /= max;
+        }
+    }
+    break;
+
+    }
+}
+
+void MathUtilities::adaptiveThreshold(std::vector<double> &data)
+{
+    int sz = int(data.size());
+    if (sz == 0) return;
+
+    std::vector<double> smoothed(sz);
+	
+    int p_pre = 8;
+    int p_post = 7;
+
+    for (int i = 0; i < sz; ++i) {
+
+        int first = std::max(0,      i - p_pre);
+        int last  = std::min(sz - 1, i + p_post);
+
+        smoothed[i] = mean(data, first, last - first + 1);
+    }
+
+    for (int i = 0; i < sz; i++) {
+        data[i] -= smoothed[i];
+        if (data[i] < 0.0) data[i] = 0.0;
+    }
+}
+
+bool
+MathUtilities::isPowerOfTwo(int x)
+{
+    if (x < 2) return false;
+    if (x & (x-1)) return false;
+    return true;
+}
+
+int
+MathUtilities::nextPowerOfTwo(int x)
+{
+    if (isPowerOfTwo(x)) return x;
+    int n = 1;
+    while (x) { x >>= 1; n <<= 1; }
+    return n;
+}
+
+int
+MathUtilities::previousPowerOfTwo(int x)
+{
+    if (isPowerOfTwo(x)) return x;
+    int n = 1;
+    x >>= 1;
+    while (x) { x >>= 1; n <<= 1; }
+    return n;
+}
+
+int
+MathUtilities::nearestPowerOfTwo(int x)
+{
+    if (isPowerOfTwo(x)) return x;
+    int n0 = previousPowerOfTwo(x), n1 = nearestPowerOfTwo(x);
+    if (x - n0 < n1 - x) return n0;
+    else return n1;
+}
+
diff --git a/libs/qm-dsp/maths/MathUtilities.h b/libs/qm-dsp/maths/MathUtilities.h
new file mode 100644
index 0000000000..4b271bf698
--- /dev/null
+++ b/libs/qm-dsp/maths/MathUtilities.h
@@ -0,0 +1,69 @@
+/* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
+
+/*
+    QM DSP Library
+
+    Centre for Digital Music, Queen Mary, University of London.
+    This file 2005-2006 Christian Landone.
+
+    This program is free software; you can redistribute it and/or
+    modify it under the terms of the GNU General Public License as
+    published by the Free Software Foundation; either version 2 of the
+    License, or (at your option) any later version.  See the file
+    COPYING included with this distribution for more information.
+*/
+
+#ifndef MATHUTILITIES_H
+#define MATHUTILITIES_H
+
+#include <vector>
+
+#include "nan-inf.h"
+
+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 mean( const std::vector<double> &data,
+                        unsigned int start, unsigned int count );
+    static double sum( const double* src, unsigned int len );
+    static double median( const double* src, unsigned int len );
+
+    static double princarg( double ang );
+    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 = 0 );
+    static int	  getMax( const std::vector<double> &data, double* max = 0 );
+    static int    compareInt(const void * a, const void * b);
+
+    enum NormaliseType {
+        NormaliseNone,
+        NormaliseUnitSum,
+        NormaliseUnitMax
+    };
+
+    static void   normalise(double *data, int length,
+                            NormaliseType n = NormaliseUnitMax);
+
+    static void   normalise(std::vector<double> &data,
+                            NormaliseType n = NormaliseUnitMax);
+
+    // moving mean threshholding:
+    static void adaptiveThreshold(std::vector<double> &data);
+
+    static bool isPowerOfTwo(int x);
+    static int nextPowerOfTwo(int x); // e.g. 1300 -> 2048, 2048 -> 2048
+    static int previousPowerOfTwo(int x); // e.g. 1300 -> 1024, 2048 -> 2048
+    static int nearestPowerOfTwo(int x); // e.g. 1300 -> 1024, 1700 -> 2048
+};
+
+#endif
diff --git a/libs/qm-dsp/maths/Polyfit.h b/libs/qm-dsp/maths/Polyfit.h
new file mode 100644
index 0000000000..86bb64cb1e
--- /dev/null
+++ b/libs/qm-dsp/maths/Polyfit.h
@@ -0,0 +1,407 @@
+/* -*- 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
+ ///////////////////////////////////////////////////////////////////////////////
+
+#include <iostream>
+
+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) {
+        std::cerr << "ERROR: PolyFit called with less than one term" << std::endl;
+        return 0;
+    }
+    if(npoints < 2) {
+        std::cerr << "ERROR: PolyFit called with less than two points" << std::endl;
+        return 0;
+    }
+    if(npoints != y.size()) {
+        std::cerr << "ERROR: PolyFit called with x and y of unequal size" << std::endl;
+        return 0;
+    }
+    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)
+	{
+            std::cerr << "ERROR: Error in PolyFit::GaussJordan: matrix is singular" << std::endl;
+            return false;
+	}
+    }
+
+    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) {
+                        std::cerr << "ERROR: Error in PolyFit::GaussJordan2: matrix is singular" << std::endl;
+                        return false;
+                    }
+
+		    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
+ 
diff --git a/libs/qm-dsp/maths/nan-inf.h b/libs/qm-dsp/maths/nan-inf.h
new file mode 100644
index 0000000000..8191713000
--- /dev/null
+++ b/libs/qm-dsp/maths/nan-inf.h
@@ -0,0 +1,20 @@
+
+#ifndef NAN_INF_H
+#define NAN_INF_H
+
+#include <math.h>
+
+#ifdef sun
+
+#include <ieeefp.h>
+#define ISNAN(x) ((sizeof(x)==sizeof(float))?isnanf(x):isnand(x))
+#define ISINF(x) (!finite(x))
+
+#else
+
+#define ISNAN(x) isnan(x)
+#define ISINF(x) isinf(x)
+
+#endif
+
+#endif
diff --git a/libs/qm-dsp/maths/pca/data.txt b/libs/qm-dsp/maths/pca/data.txt
new file mode 100644
index 0000000000..8dcfa61f2a
--- /dev/null
+++ b/libs/qm-dsp/maths/pca/data.txt
@@ -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
diff --git a/libs/qm-dsp/maths/pca/pca.c b/libs/qm-dsp/maths/pca/pca.c
new file mode 100644
index 0000000000..9dadb8687d
--- /dev/null
+++ b/libs/qm-dsp/maths/pca/pca.c
@@ -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);
+
+}
+
+
+
diff --git a/libs/qm-dsp/maths/pca/pca.h b/libs/qm-dsp/maths/pca/pca.h
new file mode 100644
index 0000000000..9c568867af
--- /dev/null
+++ b/libs/qm-dsp/maths/pca/pca.h
@@ -0,0 +1,30 @@
+#ifndef _PCA_H
+#define _PCA_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+/*
+ *  pca.h
+ *
+ *  Created by Mark Levy on 08/02/2006.
+ *  Copyright 2006 Centre for Digital Music, Queen Mary, University of London.
+
+    This program is free software; you can redistribute it and/or
+    modify it under the terms of the GNU General Public License as
+    published by the Free Software Foundation; either version 2 of the
+    License, or (at your option) any later version.  See the file
+    COPYING included with this distribution for more information.
+ *
+ */
+
+void pca_project(double** data, int n, int m, int ncomponents);
+
+#ifdef __cplusplus
+}
+#endif
+
+
+#endif
+