update qm-dsp library

11 files changed, 313 insertions, 119 deletions

diff --git a/libs/qm-dsp/maths/Correlation.cpp b/libs/qm-dsp/maths/Correlation.cpp
index fd36ac3473..17ee28f749 100644
--- a/libs/qm-dsp/maths/Correlation.cpp
+++ b/libs/qm-dsp/maths/Correlation.cpp
@@ -40,7 +40,7 @@ void Correlation::doAutoUnBiased(double *src, double *dst, unsigned int length)
     {
 	for( j = i; j < length; j++)
 	{
-	    tmp += src[ j-i ] * src[ j ];
+	    tmp += src[ j-i ] * src[ j ]; 
 	}
 
 
diff --git a/libs/qm-dsp/maths/Correlation.h b/libs/qm-dsp/maths/Correlation.h
index df2f8e17d4..85fcc7316b 100644
--- a/libs/qm-dsp/maths/Correlation.h
+++ b/libs/qm-dsp/maths/Correlation.h
@@ -18,7 +18,7 @@
 
 #define  EPS  2.2204e-016
 
-class Correlation
+class Correlation  
 {
 public:
     void doAutoUnBiased( double* src, double* dst, unsigned int length );
@@ -27,4 +27,4 @@ public:
 
 };
 
-#endif //
+#endif // 
diff --git a/libs/qm-dsp/maths/CosineDistance.cpp b/libs/qm-dsp/maths/CosineDistance.cpp
index fd2aeef3fe..13ab9ce0e8 100644
--- a/libs/qm-dsp/maths/CosineDistance.cpp
+++ b/libs/qm-dsp/maths/CosineDistance.cpp
@@ -34,7 +34,7 @@ double CosineDistance::distance(const vector<double> &v1,
     }
     else
     {
-        for(unsigned int i=0; i<v1.size(); i++)
+        for(int i=0; i<v1.size(); i++)
         {
             dSum1 += v1[i]*v2[i];
             dDen1 += v1[i]*v1[i];
diff --git a/libs/qm-dsp/maths/KLDivergence.cpp b/libs/qm-dsp/maths/KLDivergence.cpp
index 3502afb5cc..3c3cb134ca 100644
--- a/libs/qm-dsp/maths/KLDivergence.cpp
+++ b/libs/qm-dsp/maths/KLDivergence.cpp
@@ -50,7 +50,7 @@ double KLDivergence::distanceDistribution(const vector<double> &d1,
 
     double d = 0;
     double small = 1e-20;
-
+    
     for (int i = 0; i < sz; ++i) {
         d += d1[i] * log10((d1[i] + small) / (d2[i] + small));
     }
diff --git a/libs/qm-dsp/maths/MathAliases.h b/libs/qm-dsp/maths/MathAliases.h
index f597edc392..8660129cbb 100644
--- a/libs/qm-dsp/maths/MathAliases.h
+++ b/libs/qm-dsp/maths/MathAliases.h
@@ -24,7 +24,7 @@ typedef complex<double> ComplexData;
 
 
 #ifndef PI
-#define PI (3.14159265358979323846)
+#define PI (3.14159265358979232846)
 #endif
 
 #define TWO_PI 		(2. * PI)
diff --git a/libs/qm-dsp/maths/MathUtilities.cpp b/libs/qm-dsp/maths/MathUtilities.cpp
index a55252b365..cb35a2a7e0 100644
--- a/libs/qm-dsp/maths/MathUtilities.cpp
+++ b/libs/qm-dsp/maths/MathUtilities.cpp
@@ -16,6 +16,8 @@
 #include "MathUtilities.h"
 
 #include <iostream>
+#include <algorithm>
+#include <vector>
 #include <cmath>
 
 
@@ -41,11 +43,11 @@ void MathUtilities::getAlphaNorm(const double *data, unsigned int len, unsigned
     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;
@@ -60,11 +62,11 @@ double MathUtilities::getAlphaNorm( const std::vector <double> &data, unsigned 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;
@@ -75,54 +77,27 @@ double MathUtilities::getAlphaNorm( const std::vector <double> &data, unsigned i
 
 double MathUtilities::round(double x)
 {
-    double val = (double)floor(x + 0.5);
-
-    return val;
+    if (x < 0) {
+        return -floor(-x + 0.5);
+    } else {
+        return floor(x + 0.5);
+    }
 }
 
 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];
+    if (len == 0) return 0;
+    
+    std::vector<double> scratch;
+    for (int i = 0; i < len; ++i) scratch.push_back(src[i]);
+    std::sort(scratch.begin(), scratch.end());
+
+    int middle = len/2;
+    if (len % 2 == 0) {
+        return (scratch[middle] + scratch[middle - 1]) / 2;
+    } else {
+        return scratch[middle];
     }
-
-    medianVal = tempMedian;
-
-    delete [] scratch;
-    return medianVal;
 }
 
 double MathUtilities::sum(const double *src, unsigned int len)
@@ -142,8 +117,10 @@ double MathUtilities::mean(const double *src, unsigned int len)
 {
     double retVal =0.0;
 
-    double s = sum( src, len );
+    if (len == 0) return 0;
 
+    double s = sum( src, len );
+	
     retVal =  s  / (double)len;
 
     return retVal;
@@ -154,8 +131,10 @@ double MathUtilities::mean(const std::vector<double> &src,
                            unsigned int count)
 {
     double sum = 0.;
-
-    for (unsigned int i = 0; i < count; ++i)
+	
+    if (count == 0) return 0;
+    
+    for (int i = 0; i < (int)count; ++i)
     {
         sum += src[start + i];
     }
@@ -172,7 +151,7 @@ void MathUtilities::getFrameMinMax(const double *data, unsigned int len, double
         *min = *max = 0;
         return;
     }
-
+	
     *min = data[0];
     *max = data[0];
 
@@ -188,7 +167,7 @@ void MathUtilities::getFrameMinMax(const double *data, unsigned int len, double
 	{
 	    *max =  temp ;
 	}
-
+		
     }
 }
 
@@ -197,7 +176,7 @@ 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++)
@@ -209,7 +188,7 @@ int MathUtilities::getMax( double* pData, unsigned int Length, double* pMax )
 			max =  temp ;
 			index = i;
 		}
-
+		
    	}
 
 	if (pMax) *pMax = max;
@@ -223,7 +202,7 @@ 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++)
@@ -235,7 +214,7 @@ int MathUtilities::getMax( const std::vector<double> & data, double* pMax )
 			max =  temp ;
 			index = i;
 		}
-
+		
    	}
 
 	if (pMax) *pMax = max;
@@ -265,7 +244,7 @@ void MathUtilities::circShift( double* pData, int length, int shift)
 
 int MathUtilities::compareInt (const void * a, const void * b)
 {
-  return ( *(const int*)a - *(const int*)b );
+  return ( *(int*)a - *(int*)b );
 }
 
 void MathUtilities::normalise(double *data, int length, NormaliseType type)
@@ -316,9 +295,9 @@ void MathUtilities::normalise(std::vector<double> &data, NormaliseType type)
     case NormaliseUnitSum:
     {
         double sum = 0.0;
-        for (unsigned int i = 0; i < data.size(); ++i) sum += data[i];
+        for (int i = 0; i < (int)data.size(); ++i) sum += data[i];
         if (sum != 0.0) {
-            for (unsigned int i = 0; i < data.size(); ++i) data[i] /= sum;
+            for (int i = 0; i < (int)data.size(); ++i) data[i] /= sum;
         }
     }
     break;
@@ -326,11 +305,11 @@ void MathUtilities::normalise(std::vector<double> &data, NormaliseType type)
     case NormaliseUnitMax:
     {
         double max = 0.0;
-        for (unsigned int i = 0; i < data.size(); ++i) {
+        for (int i = 0; i < (int)data.size(); ++i) {
             if (fabs(data[i]) > max) max = fabs(data[i]);
         }
         if (max != 0.0) {
-            for (unsigned int i = 0; i < data.size(); ++i) data[i] /= max;
+            for (int i = 0; i < (int)data.size(); ++i) data[i] /= max;
         }
     }
     break;
@@ -344,7 +323,7 @@ void MathUtilities::adaptiveThreshold(std::vector<double> &data)
     if (sz == 0) return;
 
     std::vector<double> smoothed(sz);
-
+	
     int p_pre = 8;
     int p_post = 7;
 
@@ -365,7 +344,7 @@ void MathUtilities::adaptiveThreshold(std::vector<double> &data)
 bool
 MathUtilities::isPowerOfTwo(int x)
 {
-    if (x < 2) return false;
+    if (x < 1) return false;
     if (x & (x-1)) return false;
     return true;
 }
@@ -374,6 +353,7 @@ int
 MathUtilities::nextPowerOfTwo(int x)
 {
     if (isPowerOfTwo(x)) return x;
+    if (x < 1) return 1;
     int n = 1;
     while (x) { x >>= 1; n <<= 1; }
     return n;
@@ -383,6 +363,7 @@ int
 MathUtilities::previousPowerOfTwo(int x)
 {
     if (isPowerOfTwo(x)) return x;
+    if (x < 1) return 1;
     int n = 1;
     x >>= 1;
     while (x) { x >>= 1; n <<= 1; }
@@ -393,8 +374,30 @@ int
 MathUtilities::nearestPowerOfTwo(int x)
 {
     if (isPowerOfTwo(x)) return x;
-    int n0 = previousPowerOfTwo(x), n1 = nearestPowerOfTwo(x);
+    int n0 = previousPowerOfTwo(x), n1 = nextPowerOfTwo(x);
     if (x - n0 < n1 - x) return n0;
     else return n1;
 }
 
+double
+MathUtilities::factorial(int x)
+{
+    if (x < 0) return 0;
+    double f = 1;
+    for (int i = 1; i <= x; ++i) {
+	f = f * i;
+    }
+    return f;
+}
+
+int
+MathUtilities::gcd(int a, int b)
+{
+    int c = a % b;
+    if (c == 0) {
+        return b;
+    } else {
+        return gcd(b, c);
+    }
+}
+
diff --git a/libs/qm-dsp/maths/MathUtilities.h b/libs/qm-dsp/maths/MathUtilities.h
index 85774756d4..fac710af9a 100644
--- a/libs/qm-dsp/maths/MathUtilities.h
+++ b/libs/qm-dsp/maths/MathUtilities.h
@@ -20,20 +20,57 @@
 
 #include "nan-inf.h"
 
-class MathUtilities
+/**
+ * Static helper functions for simple mathematical calculations.
+ */
+class MathUtilities  
 {
-public:
+public:	
+    /**
+     * Round x to the nearest integer.
+     */
     static double round( double x );
 
+    /**
+     * Return through min and max pointers the highest and lowest
+     * values in the given array of the given length.
+     */
     static void	  getFrameMinMax( const double* data, unsigned int len,  double* min, double* max );
 
+    /**
+     * Return the mean of the given array of the given length.
+     */
     static double mean( const double* src, unsigned int len );
+
+    /**
+     * Return the mean of the subset of the given vector identified by
+     * start and count.
+     */
     static double mean( const std::vector<double> &data,
                         unsigned int start, unsigned int count );
+    
+    /**
+     * Return the sum of the values in the given array of the given
+     * length.
+     */
     static double sum( const double* src, unsigned int len );
+
+    /**
+     * Return the median of the values in the given array of the given
+     * length. If the array is even in length, the returned value will
+     * be half-way between the two values adjacent to median.
+     */
     static double median( const double* src, unsigned int len );
 
+    /**
+     * The principle argument function. Map the phase angle ang into
+     * the range [-pi,pi).
+     */
     static double princarg( double ang );
+
+    /**
+     * Floating-point division modulus: return x % y.
+     */
     static double mod( double x, double y);
 
     static void	  getAlphaNorm(const double *data, unsigned int len, unsigned int alpha, double* ANorm);
@@ -51,19 +88,50 @@ public:
         NormaliseUnitMax
     };
 
-    static void   normalise(double *data, int length,
-                            NormaliseType n = NormaliseUnitMax);
+    static void normalise(double *data, int length,
+                          NormaliseType n = NormaliseUnitMax);
 
-    static void   normalise(std::vector<double> &data,
-                            NormaliseType n = NormaliseUnitMax);
+    static void normalise(std::vector<double> &data,
+                          NormaliseType n = NormaliseUnitMax);
 
-    // moving mean threshholding:
+    /**
+     * Threshold the input/output vector data against a moving-mean
+     * average filter.
+     */
     static void adaptiveThreshold(std::vector<double> &data);
 
+    /** 
+     * Return true if x is 2^n for some integer n >= 0.
+     */
     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
+
+    /**
+     * Return the next higher integer power of two from x, e.g. 1300
+     * -> 2048, 2048 -> 2048.
+     */
+    static int nextPowerOfTwo(int x);
+
+    /**
+     * Return the next lower integer power of two from x, e.g. 1300 ->
+     * 1024, 2048 -> 2048.
+     */
+    static int previousPowerOfTwo(int x);
+
+    /**
+     * Return the nearest integer power of two to x, e.g. 1300 -> 1024,
+     * 12 -> 16 (not 8; if two are equidistant, the higher is returned).
+     */
+    static int nearestPowerOfTwo(int x);
+
+    /**
+     * Return x!
+     */
+    static double factorial(int x); // returns double in case it is large
+
+    /**
+     * Return the greatest common divisor of natural numbers a and b.
+     */
+    static int gcd(int a, int b);
 };
 
 #endif
diff --git a/libs/qm-dsp/maths/MedianFilter.h b/libs/qm-dsp/maths/MedianFilter.h
new file mode 100644
index 0000000000..6c01cc4ab5
--- /dev/null
+++ b/libs/qm-dsp/maths/MedianFilter.h
@@ -0,0 +1,131 @@
+/* -*- 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 2010 Chris Cannam.
+
+    This program is free software; you can redistribute it and/or
+    modify it under the terms of the GNU General Public License as
+    published by the Free Software Foundation; either version 2 of the
+    License, or (at your option) any later version.  See the file
+    COPYING included with this distribution for more information.
+*/
+
+#ifndef MEDIAN_FILTER_H
+#define MEDIAN_FILTER_H
+
+#include <algorithm>
+#include <cassert>
+#include <cmath>
+#include <iostream>
+#include <vector>
+
+template <typename T>
+class MedianFilter
+{
+public:
+    MedianFilter(int size, float percentile = 50.f) :
+	m_size(size),
+	m_frame(new T[size]),
+	m_sorted(new T[size]),
+	m_sortend(m_sorted + size - 1) {
+        setPercentile(percentile);
+	reset();
+    }
+
+    ~MedianFilter() { 
+	delete[] m_frame;
+	delete[] m_sorted;
+    }
+
+    void setPercentile(float p) {
+        m_index = int((m_size * p) / 100.f);
+        if (m_index >= m_size) m_index = m_size-1;
+        if (m_index < 0) m_index = 0;
+    }
+
+    void push(T value) {
+        if (value != value) {
+            std::cerr << "WARNING: MedianFilter::push: attempt to push NaN, pushing zero instead" << std::endl;
+            // we do need to push something, to maintain the filter length
+            value = T();
+        }
+	drop(m_frame[0]);
+	const int sz1 = m_size-1;
+	for (int i = 0; i < sz1; ++i) m_frame[i] = m_frame[i+1];
+	m_frame[m_size-1] = value;
+	put(value);
+    }
+
+    T get() const {
+	return m_sorted[m_index];
+    }
+
+    int getSize() const {
+        return m_size; 
+    }
+
+    T getAt(float percentile) {
+	int ix = int((m_size * percentile) / 100.f);
+        if (ix >= m_size) ix = m_size-1;
+        if (ix < 0) ix = 0;
+	return m_sorted[ix];
+    }
+
+    void reset() {
+	for (int i = 0; i < m_size; ++i) m_frame[i] = 0;
+	for (int i = 0; i < m_size; ++i) m_sorted[i] = 0;
+    }
+
+    static std::vector<T> filter(int size, const std::vector<T> &in) {
+        std::vector<T> out;
+        MedianFilter<T> f(size);
+        for (int i = 0; i < int(in.size()); ++i) {
+            f.push(in[i]);
+            T median = f.get();
+            if (i >= size/2) out.push_back(median);
+        }
+        while (out.size() < in.size()) {
+            f.push(T());
+            out.push_back(f.get());
+        }
+        return out;
+    }
+
+private:
+    const int m_size;
+    T *const m_frame;
+    T *const m_sorted;
+    T *const m_sortend;
+    int m_index;
+
+    void put(T value) {
+	// precondition: m_sorted contains m_size-1 values, packed at start
+	// postcondition: m_sorted contains m_size values, one of which is value
+	T *point = std::lower_bound(m_sorted, m_sortend, value);
+	const int n = m_sortend - point;
+	for (int i = n; i > 0; --i) point[i] = point[i-1];
+	*point = value;
+    }
+
+    void drop(T value) {
+	// precondition: m_sorted contains m_size values, one of which is value
+	// postcondition: m_sorted contains m_size-1 values, packed at start
+	T *point = std::lower_bound(m_sorted, m_sortend + 1, value);
+	if (*point != value) {
+	    std::cerr << "WARNING: MedianFilter::drop: *point is " << *point
+		      << ", expected " << value << std::endl;
+	}
+	const int n = m_sortend - point;
+	for (int i = 0; i < n; ++i) point[i] = point[i+1];
+	*m_sortend = T(0);
+    }
+
+    MedianFilter(const MedianFilter &); // not provided
+    MedianFilter &operator=(const MedianFilter &); // not provided
+};
+
+#endif
+
diff --git a/libs/qm-dsp/maths/Polyfit.h b/libs/qm-dsp/maths/Polyfit.h
index a5c2ca2af4..5cf97d9df0 100644
--- a/libs/qm-dsp/maths/Polyfit.h
+++ b/libs/qm-dsp/maths/Polyfit.h
@@ -53,13 +53,13 @@ public:
 			    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
@@ -105,13 +105,13 @@ double TPolyFit::PolyFit2 (const vector<double> &x,
 // nterms = coefs.size()
 // npoints = x.size()
 {
-    unsigned int i, j;
+    int i, j;
     double xi, yi, yc, srs, sum_y, sum_y2;
     Matrix xmatr;        // Data matrix
     Matrix a;
     vector<double> g;      // Constant vector
-    const unsigned int npoints(x.size());
-    const unsigned int nterms(coefs.size());
+    const int npoints(x.size());
+    const int nterms(coefs.size());
     double correl_coef;
     zeroise(g, nterms);
     zeroise(a, nterms, nterms);
@@ -124,7 +124,7 @@ double TPolyFit::PolyFit2 (const vector<double> &x,
         std::cerr << "ERROR: PolyFit called with less than two points" << std::endl;
         return 0;
     }
-    if(npoints != y.size()) {
+    if(npoints != (int)y.size()) {
         std::cerr << "ERROR: PolyFit called with x and y of unequal size" << std::endl;
         return 0;
     }
@@ -260,8 +260,8 @@ bool TPolyFit::GaussJordan (Matrix &b,
 
     for( int i = 0; i < ncol; ++i)
 	coef[i] = w[i][0];
-
-
+ 
+ 
     return true;
 }   // end;	{ procedure GaussJordan }
 //----------------------------------------------------------------------------------------------
@@ -274,12 +274,11 @@ bool TPolyFit::GaussJordan2(Matrix &b,
 {
     //GaussJordan2;         // first half of GaussJordan
     // actual start of gaussj
-
+ 
     double big, t;
     double pivot;
     double determ;
-    int irow = 0;
-    int icol = 0;
+    int irow, icol;
     int ncol(b.size());
     int nv = 1;                  // single constant vector
     for(int i = 0; i < ncol; ++i)
@@ -405,4 +404,4 @@ void NSUtility::zeroise(vector<vector<int> > &matrix, int m, int n)
 
 
 #endif
-
+ 
diff --git a/libs/qm-dsp/maths/nan-inf.h b/libs/qm-dsp/maths/nan-inf.h
index 8191713000..c6be51276a 100644
--- a/libs/qm-dsp/maths/nan-inf.h
+++ b/libs/qm-dsp/maths/nan-inf.h
@@ -2,19 +2,12 @@
 #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
+#define ISNAN(x) (sizeof(x) == sizeof(double) ? ISNANd(x) : ISNANf(x))
+static inline int ISNANf(float x) { return x != x; }
+static inline int ISNANd(double x) { return x != x; }
+          
+#define ISINF(x) (sizeof(x) == sizeof(double) ? ISINFd(x) : ISINFf(x))
+static inline int ISINFf(float x) { return !ISNANf(x) && ISNANf(x - x); }
+static inline int ISINFd(double x) { return !ISNANd(x) && ISNANd(x - x); }
 
 #endif
diff --git a/libs/qm-dsp/maths/pca/pca.c b/libs/qm-dsp/maths/pca/pca.c
index d2b2f11c47..9dadb8687d 100644
--- a/libs/qm-dsp/maths/pca/pca.c
+++ b/libs/qm-dsp/maths/pca/pca.c
@@ -15,8 +15,8 @@
    Earn/Bitnet:  fionn@dgaeso51,  fim@dgaipp1s,  murtagh@stsci
    Span:         esomc1::fionn
    Internet:     murtagh@scivax.stsci.edu
-
-   F. Murtagh, Munich, 6 June 1989                                   */
+   
+   F. Murtagh, Munich, 6 June 1989                                   */   
 /*********************************************************************/
 
 #include <stdio.h>
@@ -110,7 +110,7 @@ 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;
@@ -188,7 +188,7 @@ 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;
@@ -253,23 +253,23 @@ 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);
+    //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 */
@@ -278,7 +278,7 @@ void pca_project(double** data, int n, int m, int ncomponents)
     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. */
+columns of symmat now contain the associated eigenvectors. */	
 
 /*
 	printf("\nEigenvalues:\n");
@@ -289,7 +289,7 @@ columns of symmat now contain the associated eigenvectors. */
 	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++) {
@@ -310,7 +310,7 @@ for (i = 0; i < n; i++) {
         }
 }
 
-/*
+/*	
 printf("\nProjections of row-points on first 3 prin. comps.:\n");
  for (i = 0; i < n; i++) {
 	 for (j = 0; j < 3; j++)  {