Update qm-dsp library (v1.7.1-20-g4d15479)

5 files changed, 134 insertions, 111 deletions

diff --git a/libs/qm-dsp/maths/CosineDistance.cpp b/libs/qm-dsp/maths/CosineDistance.cpp
index 13ab9ce0e8..4b06a7ad2a 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(int i=0; i<v1.size(); i++)
+        for(int i=0; i<int(v1.size()); i++)
         {
             dSum1 += v1[i]*v2[i];
             dDen1 += v1[i]*v1[i];
diff --git a/libs/qm-dsp/maths/MathUtilities.cpp b/libs/qm-dsp/maths/MathUtilities.cpp
index cb35a2a7e0..a9dabf372d 100644
--- a/libs/qm-dsp/maths/MathUtilities.cpp
+++ b/libs/qm-dsp/maths/MathUtilities.cpp
@@ -20,6 +20,7 @@
 #include <vector>
 #include <cmath>
 
+using namespace std;
 
 double MathUtilities::mod(double x, double y)
 {
@@ -38,9 +39,9 @@ double MathUtilities::princarg(double ang)
     return ValOut;
 }
 
-void MathUtilities::getAlphaNorm(const double *data, unsigned int len, unsigned int alpha, double* ANorm)
+void MathUtilities::getAlphaNorm(const double *data, int len, int alpha, double* ANorm)
 {
-    unsigned int i;
+    int i;
     double temp = 0.0;
     double a=0.0;
 	
@@ -56,17 +57,16 @@ void MathUtilities::getAlphaNorm(const double *data, unsigned int len, unsigned
     *ANorm = a;
 }
 
-double MathUtilities::getAlphaNorm( const std::vector <double> &data, unsigned int alpha )
+double MathUtilities::getAlphaNorm( const vector <double> &data, int alpha )
 {
-    unsigned int i;
-    unsigned int len = data.size();
+    int i;
+    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;
@@ -84,13 +84,13 @@ double MathUtilities::round(double x)
     }
 }
 
-double MathUtilities::median(const double *src, unsigned int len)
+double MathUtilities::median(const double *src, int len)
 {
     if (len == 0) return 0;
     
-    std::vector<double> scratch;
+    vector<double> scratch;
     for (int i = 0; i < len; ++i) scratch.push_back(src[i]);
-    std::sort(scratch.begin(), scratch.end());
+    sort(scratch.begin(), scratch.end());
 
     int middle = len/2;
     if (len % 2 == 0) {
@@ -100,9 +100,9 @@ double MathUtilities::median(const double *src, unsigned int len)
     }
 }
 
-double MathUtilities::sum(const double *src, unsigned int len)
+double MathUtilities::sum(const double *src, int len)
 {
-    unsigned int i ;
+    int i ;
     double retVal =0.0;
 
     for(  i = 0; i < len; i++)
@@ -113,7 +113,7 @@ double MathUtilities::sum(const double *src, unsigned int len)
     return retVal;
 }
 
-double MathUtilities::mean(const double *src, unsigned int len)
+double MathUtilities::mean(const double *src, int len)
 {
     double retVal =0.0;
 
@@ -126,9 +126,9 @@ double MathUtilities::mean(const double *src, unsigned int len)
     return retVal;
 }
 
-double MathUtilities::mean(const std::vector<double> &src,
-                           unsigned int start,
-                           unsigned int count)
+double MathUtilities::mean(const vector<double> &src,
+                           int start,
+                           int count)
 {
     double sum = 0.;
 	
@@ -142,9 +142,9 @@ double MathUtilities::mean(const std::vector<double> &src,
     return sum / count;
 }
 
-void MathUtilities::getFrameMinMax(const double *data, unsigned int len, double *min, double *max)
+void MathUtilities::getFrameMinMax(const double *data, int len, double *min, double *max)
 {
-    unsigned int i;
+    int i;
     double temp = 0.0;
 
     if (len == 0) {
@@ -171,10 +171,10 @@ void MathUtilities::getFrameMinMax(const double *data, unsigned int len, double
     }
 }
 
-int MathUtilities::getMax( double* pData, unsigned int Length, double* pMax )
+int MathUtilities::getMax( double* pData, int Length, double* pMax )
 {
-	unsigned int index = 0;
-	unsigned int i;
+	int index = 0;
+	int i;
 	double temp = 0.0;
 	
 	double max = pData[0];
@@ -197,15 +197,15 @@ int MathUtilities::getMax( double* pData, unsigned int Length, double* pMax )
 	return index;
 }
 
-int MathUtilities::getMax( const std::vector<double> & data, double* pMax )
+int MathUtilities::getMax( const vector<double> & data, double* pMax )
 {
-	unsigned int index = 0;
-	unsigned int i;
+	int index = 0;
+	int i;
 	double temp = 0.0;
 	
 	double max = data[0];
 
-	for( i = 0; i < data.size(); i++)
+	for( i = 0; i < int(data.size()); i++)
 	{
 		temp = data[ i ];
 
@@ -286,7 +286,7 @@ void MathUtilities::normalise(double *data, int length, NormaliseType type)
     }
 }
 
-void MathUtilities::normalise(std::vector<double> &data, NormaliseType type)
+void MathUtilities::normalise(vector<double> &data, NormaliseType type)
 {
     switch (type) {
 
@@ -317,20 +317,46 @@ void MathUtilities::normalise(std::vector<double> &data, NormaliseType type)
     }
 }
 
-void MathUtilities::adaptiveThreshold(std::vector<double> &data)
+double MathUtilities::getLpNorm(const vector<double> &data, int p)
+{
+    double tot = 0.0;
+    for (int i = 0; i < int(data.size()); ++i) {
+        tot += abs(pow(data[i], p));
+    }
+    return pow(tot, 1.0 / p);
+}
+
+vector<double> MathUtilities::normaliseLp(const vector<double> &data,
+                                               int p,
+                                               double threshold)
+{
+    int n = int(data.size());
+    if (n == 0 || p == 0) return data;
+    double norm = getLpNorm(data, p);
+    if (norm < threshold) {
+        return vector<double>(n, 1.0 / pow(n, 1.0 / p)); // unit vector
+    }
+    vector<double> out(n);
+    for (int i = 0; i < n; ++i) {
+        out[i] = data[i] / norm;
+    }
+    return out;
+}
+    
+void MathUtilities::adaptiveThreshold(vector<double> &data)
 {
     int sz = int(data.size());
     if (sz == 0) return;
 
-    std::vector<double> smoothed(sz);
+    vector<double> smoothed(sz);
 	
     int p_pre = 8;
     int p_post = 7;
 
     for (int i = 0; i < sz; ++i) {
 
-        int first = std::max(0,      i - p_pre);
-        int last  = std::min(sz - 1, i + p_post);
+        int first = max(0,      i - p_pre);
+        int last  = min(sz - 1, i + p_post);
 
         smoothed[i] = mean(data, first, last - first + 1);
     }
diff --git a/libs/qm-dsp/maths/MathUtilities.h b/libs/qm-dsp/maths/MathUtilities.h
index fac710af9a..ec86efcf2f 100644
--- a/libs/qm-dsp/maths/MathUtilities.h
+++ b/libs/qm-dsp/maths/MathUtilities.h
@@ -35,32 +35,32 @@ public:
      * 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 );
+    static void	  getFrameMinMax( const double* data, 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 );
+    static double mean( const double* src, 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 );
+                        int start, 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 );
+    static double sum( const double* src, 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 );
+    static double median( const double* src, int len );
 
     /**
      * The principle argument function. Map the phase angle ang into
@@ -73,14 +73,21 @@ public:
      */
     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);
+    /**
+     * The alpha norm is the alpha'th root of the mean alpha'th power
+     * magnitude. For example if alpha = 2 this corresponds to the RMS
+     * of the input data, and when alpha = 1 this is the mean
+     * magnitude.
+     */
+    static void	  getAlphaNorm(const double *data, int len, int alpha, double* ANorm);
 
-    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);
+    /**
+     * The alpha norm is the alpha'th root of the mean alpha'th power
+     * magnitude. For example if alpha = 2 this corresponds to the RMS
+     * of the input data, and when alpha = 1 this is the mean
+     * magnitude.
+     */
+    static double getAlphaNorm(const std::vector <double> &data, int alpha );
 
     enum NormaliseType {
         NormaliseNone,
@@ -95,11 +102,34 @@ public:
                           NormaliseType n = NormaliseUnitMax);
 
     /**
+     * Calculate the L^p norm of a vector. Equivalent to MATLAB's
+     * norm(data, p).
+     */
+    static double getLpNorm(const std::vector<double> &data,
+                            int p);
+
+    /**
+     * Normalise a vector by dividing through by its L^p norm. If the
+     * norm is below the given threshold, the unit vector for that
+     * norm is returned. p may be 0, in which case no normalisation
+     * happens and the data is returned unchanged.
+     */
+    static std::vector<double> normaliseLp(const std::vector<double> &data,
+                                           int p,
+                                           double threshold = 1e-6);
+    
+    /**
      * Threshold the input/output vector data against a moving-mean
      * average filter.
      */
     static void adaptiveThreshold(std::vector<double> &data);
 
+    static void   circShift( double* data, int length, int shift);
+
+    static int	  getMax( double* data, int length, double* max = 0 );
+    static int	  getMax( const std::vector<double> &data, double* max = 0 );
+    static int    compareInt(const void * a, const void * b);
+
     /** 
      * Return true if x is 2^n for some integer n >= 0.
      */
diff --git a/libs/qm-dsp/maths/Polyfit.h b/libs/qm-dsp/maths/Polyfit.h
index 5cf97d9df0..3405b5cb70 100644
--- a/libs/qm-dsp/maths/Polyfit.h
+++ b/libs/qm-dsp/maths/Polyfit.h
@@ -80,22 +80,36 @@ private:
 
 // some utility functions
 
-namespace NSUtility
+struct 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;}
+    static void swap(double &a, double &b) {double t = a; a = b; b = t;}
+    // fills a vector with zeros.
+    static void zeroise(vector<double> &array, int n) { 
+        array.clear();
+        for(int j = 0; j < n; ++j) array.push_back(0);
+    }
+    // fills a vector with zeros.
+    static void 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.
+    static void 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.
+    static void 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);
+    }
+    static double sqr(const double &x) {return x * x;}
 };
 
-//---------------------------------------------------------------------------
-// Implementation
-//---------------------------------------------------------------------------
-using namespace NSUtility;
-//------------------------------------------------------------------------------------------
-
 
 // main PolyFit routine
 
@@ -113,9 +127,9 @@ double TPolyFit::PolyFit2 (const vector<double> &x,
     const int npoints(x.size());
     const int nterms(coefs.size());
     double correl_coef;
-    zeroise(g, nterms);
-    zeroise(a, nterms, nterms);
-    zeroise(xmatr, npoints, nterms);
+    NSUtility::zeroise(g, nterms);
+    NSUtility::zeroise(a, nterms, nterms);
+    NSUtility::zeroise(xmatr, npoints, nterms);
     if (nterms < 1) {
         std::cerr << "ERROR: PolyFit called with less than one term" << std::endl;
         return 0;
@@ -148,13 +162,13 @@ double TPolyFit::PolyFit2 (const vector<double> &x,
 	yc = 0.0;
 	for(j = 0; j < nterms; ++j)
 	    yc += coefs [j] * xmatr [i][j];
-	srs += sqr (yc - yi);
+	srs += NSUtility::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;
+    correl_coef = sum_y2 - NSUtility::sqr (sum_y) / npoints;
     // Either return 0 or the correct value of correlation coefficient
     if (correl_coef != 0)
 	correl_coef = srs / correl_coef;
@@ -229,8 +243,8 @@ bool TPolyFit::GaussJordan (Matrix &b,
     vector<vector<int> >index;
     Matrix w;
 
-    zeroise(w, ncol, ncol);
-    zeroise(index, ncol, 3);
+    NSUtility::zeroise(w, ncol, ncol);
+    NSUtility::zeroise(index, ncol, 3);
 
     if(!GaussJordan2(b, y, w, index))
 	return false;
@@ -278,7 +292,7 @@ bool TPolyFit::GaussJordan2(Matrix &b,
     double big, t;
     double pivot;
     double determ;
-    int irow, icol;
+    int irow = 0, icol = 0;
     int ncol(b.size());
     int nv = 1;                  // single constant vector
     for(int i = 0; i < ncol; ++i)
@@ -355,53 +369,6 @@ bool TPolyFit::GaussJordan2(Matrix &b,
     } // { 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/pca/pca.c b/libs/qm-dsp/maths/pca/pca.c
index 9dadb8687d..1a7bfdd549 100644
--- a/libs/qm-dsp/maths/pca/pca.c
+++ b/libs/qm-dsp/maths/pca/pca.c
@@ -252,7 +252,7 @@ void tqli(double* d, double* e, int n, double** z)
 void pca_project(double** data, int n, int m, int ncomponents)
 {
 	int  i, j, k, k2;
-	double  **symmat, **symmat2, *evals, *interm;
+	double  **symmat, /* **symmat2, */ *evals, *interm;
 	
 	//TODO: assert ncomponents < m