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

25 files changed, 373 insertions, 339 deletions

diff --git a/libs/qm-dsp/README.txt b/libs/qm-dsp/README.txt
index 557ac23294..5856500f49 100644
--- a/libs/qm-dsp/README.txt
+++ b/libs/qm-dsp/README.txt
@@ -1,7 +1,6 @@
 This is a stripped down version of the Queen Mary DSP library
-Version 1.7.1 from Sept 2015.
 
-https://code.soundsoftware.ac.uk/attachments/download/1582/qm-dsp-1.7.1.tar.gz
+https://github.com/c4dm/qm-dsp  -- see gitrev.txt for version
 
 ---
 
diff --git a/libs/qm-dsp/base/KaiserWindow.h b/libs/qm-dsp/base/KaiserWindow.h
index 0253d6d4c4..f16a4b6c16 100644
--- a/libs/qm-dsp/base/KaiserWindow.h
+++ b/libs/qm-dsp/base/KaiserWindow.h
@@ -81,7 +81,7 @@ public:
     }
 
     const double *getWindow() const { 
-	return &m_window[0];
+	return m_window.data();
     }
 
     void cut(double *src) const { 
diff --git a/libs/qm-dsp/base/SincWindow.h b/libs/qm-dsp/base/SincWindow.h
index 02de928679..bb35d90c20 100644
--- a/libs/qm-dsp/base/SincWindow.h
+++ b/libs/qm-dsp/base/SincWindow.h
@@ -37,7 +37,7 @@ public:
     }
 
     const double *getWindow() const { 
-	return &m_window[0];
+	return m_window.data();
     }
 
     void cut(double *src) const { 
diff --git a/libs/qm-dsp/dsp/chromagram/Chromagram.cpp b/libs/qm-dsp/dsp/chromagram/Chromagram.cpp
index a8597a5ddd..f42a4d2aad 100644
--- a/libs/qm-dsp/dsp/chromagram/Chromagram.cpp
+++ b/libs/qm-dsp/dsp/chromagram/Chromagram.cpp
@@ -34,7 +34,7 @@ int Chromagram::initialise( ChromaConfig Config )
     m_normalise = Config.normalise;     // if frame normalisation is required
 
     // No. of constant Q bins
-    m_uK = ( unsigned int ) ceil( m_BPO * log(m_FMax/m_FMin)/log(2.0));	
+    m_uK = (int) ceil( m_BPO * log(m_FMax/m_FMin)/log(2.0));	
 
     // Create array for chroma result
     m_chromadata = new double[ m_BPO ];
@@ -112,7 +112,7 @@ void Chromagram::unityNormalise(double *src)
 
     MathUtilities::getFrameMinMax( src, m_BPO, & min, &max );
 
-    for( unsigned int i = 0; i < m_BPO; i++ )
+    for (int i = 0; i < m_BPO; i++)
     {
 	val = src[ i ] / max;
 
@@ -153,19 +153,17 @@ double* Chromagram::process( const double *real, const double *imag )
     }
 
     // initialise chromadata to 0
-    for (unsigned i = 0; i < m_BPO; i++) m_chromadata[i] = 0;
+    for (int i = 0; i < m_BPO; i++) m_chromadata[i] = 0;
 
-    double cmax = 0.0;
-    double cval = 0;
     // Calculate ConstantQ frame
     m_ConstantQ->process( real, imag, m_CQRe, m_CQIm );
 	
     // add each octave of cq data into Chromagram
-    const unsigned octaves = (int)floor(double( m_uK/m_BPO))-1;
-    for (unsigned octave = 0; octave <= octaves; octave++) 
+    const int octaves = (int)floor(double( m_uK/m_BPO))-1;
+    for (int octave = 0; octave <= octaves; octave++) 
     {
-	unsigned firstBin = octave*m_BPO;
-	for (unsigned i = 0; i < m_BPO; i++) 
+	int firstBin = octave*m_BPO;
+	for (int i = 0; i < m_BPO; i++) 
 	{
 	    m_chromadata[i] += kabs( m_CQRe[ firstBin + i ], m_CQIm[ firstBin + i ]);
 	}
diff --git a/libs/qm-dsp/dsp/chromagram/Chromagram.h b/libs/qm-dsp/dsp/chromagram/Chromagram.h
index bd928f5d48..790abd8e23 100644
--- a/libs/qm-dsp/dsp/chromagram/Chromagram.h
+++ b/libs/qm-dsp/dsp/chromagram/Chromagram.h
@@ -21,10 +21,10 @@
 #include "ConstantQ.h"
 
 struct ChromaConfig{
-    unsigned int FS;
+    int FS;
     double min;
     double max;
-    unsigned int BPO;
+    int BPO;
     double CQThresh;
     MathUtilities::NormaliseType normalise;
 };
@@ -44,10 +44,10 @@ public:
     double kabs( double real, double imag );
 	
     // Results
-    unsigned int getK() { return m_uK;}
-    unsigned int getFrameSize() { return m_frameSize; }
-    unsigned int getHopSize()   { return m_hopSize; }
-
+    int getK() { return m_uK;}
+    int getFrameSize() { return m_frameSize; }
+    int getHopSize()   { return m_hopSize; }
+    
 private:
     int initialise( ChromaConfig Config );
     int deInitialise();
@@ -58,13 +58,13 @@ private:
     double* m_chromadata;
     double m_FMin;
     double m_FMax;
-    unsigned int m_BPO;
-    unsigned int m_uK;
+    int m_BPO;
+    int m_uK;
 
     MathUtilities::NormaliseType m_normalise;
 
-    unsigned int m_frameSize;
-    unsigned int m_hopSize;
+    int m_frameSize;
+    int m_hopSize;
 
     FFTReal* m_FFT;
     ConstantQ* m_ConstantQ;
diff --git a/libs/qm-dsp/dsp/rateconversion/DecimatorB.cpp b/libs/qm-dsp/dsp/rateconversion/DecimatorB.cpp
index 55df1e9fc0..5f27130c49 100644
--- a/libs/qm-dsp/dsp/rateconversion/DecimatorB.cpp
+++ b/libs/qm-dsp/dsp/rateconversion/DecimatorB.cpp
@@ -112,7 +112,6 @@ void DecimatorB::doProcess()
 {
     int filteridx = 0;
     int factorDone = 1;
-    int factorRemaining = m_decFactor;
 
     while (factorDone < m_decFactor) {
 
diff --git a/libs/qm-dsp/dsp/signalconditioning/DFProcess.cpp b/libs/qm-dsp/dsp/signalconditioning/DFProcess.cpp
index 52443fb384..7aa9649e6c 100644
--- a/libs/qm-dsp/dsp/signalconditioning/DFProcess.cpp
+++ b/libs/qm-dsp/dsp/signalconditioning/DFProcess.cpp
@@ -59,13 +59,11 @@ void DFProcess::initialise( DFProcConfig Config )
     filtSrc = new double[ m_length ];
     filtDst = new double[ m_length ];
 
-	
-    //Low Pass Smoothing Filter Config
-    m_FilterConfigParams.ord = Config.LPOrd;
-    m_FilterConfigParams.ACoeffs = Config.LPACoeffs;
-    m_FilterConfigParams.BCoeffs = Config.LPBCoeffs;
-	
-    m_FiltFilt = new FiltFilt( m_FilterConfigParams );
+    Filter::Parameters params;
+    params.a = std::vector<double>(Config.LPACoeffs, Config.LPACoeffs + Config.LPOrd + 1);
+    params.b = std::vector<double>(Config.LPBCoeffs, Config.LPBCoeffs + Config.LPOrd + 1);
+    
+    m_FiltFilt = new FiltFilt(params);
 	
     //add delta threshold
     m_delta = Config.delta;
@@ -193,7 +191,7 @@ void DFProcess::removeDCNormalize( double *src, double*dst )
 
     MathUtilities::getAlphaNorm( src, m_length, m_alphaNormParam, &DFAlphaNorm );
 
-    for( unsigned int i = 0; i< m_length; i++)
+    for (int i = 0; i < m_length; i++)
     {
 	dst[ i ] = ( src[ i ] - DFMin ) / DFAlphaNorm; 
     }
diff --git a/libs/qm-dsp/dsp/signalconditioning/DFProcess.h b/libs/qm-dsp/dsp/signalconditioning/DFProcess.h
index b93b535518..b256eca16e 100644
--- a/libs/qm-dsp/dsp/signalconditioning/DFProcess.h
+++ b/libs/qm-dsp/dsp/signalconditioning/DFProcess.h
@@ -81,8 +81,6 @@ private:
     double* m_filtScratchIn;
     double* m_filtScratchOut;
 
-    FilterConfig m_FilterConfigParams;
-
     FiltFilt* m_FiltFilt;
 
     bool m_isMedianPositive;
diff --git a/libs/qm-dsp/dsp/signalconditioning/FiltFilt.cpp b/libs/qm-dsp/dsp/signalconditioning/FiltFilt.cpp
index 03b21138a4..c88641de7b 100644
--- a/libs/qm-dsp/dsp/signalconditioning/FiltFilt.cpp
+++ b/libs/qm-dsp/dsp/signalconditioning/FiltFilt.cpp
@@ -13,43 +13,22 @@
     COPYING included with this distribution for more information.
 */
 
-#include <stdio.h>
 #include "FiltFilt.h"
 
 //////////////////////////////////////////////////////////////////////
 // Construction/Destruction
 //////////////////////////////////////////////////////////////////////
 
-FiltFilt::FiltFilt( FilterConfig Config )
+FiltFilt::FiltFilt(Filter::Parameters parameters) :
+    m_filter(parameters)
 {
-    m_filtScratchIn = NULL;
-    m_filtScratchOut = NULL;
-    m_ord = 0;
-	
-    initialise( Config );
+    m_ord = m_filter.getOrder();
 }
 
 FiltFilt::~FiltFilt()
 {
-    deInitialise();
-}
-
-void FiltFilt::initialise( FilterConfig Config )
-{
-    m_ord = Config.ord;
-    m_filterConfig.ord = Config.ord;
-    m_filterConfig.ACoeffs = Config.ACoeffs;
-    m_filterConfig.BCoeffs = Config.BCoeffs;
-	
-    m_filter = new Filter( m_filterConfig );
-}
-
-void FiltFilt::deInitialise()
-{
-    delete m_filter;
 }
 
-
 void FiltFilt::process(double *src, double *dst, unsigned int length)
 {	
     unsigned int i;
@@ -57,7 +36,6 @@ void FiltFilt::process(double *src, double *dst, unsigned int length)
     if (length == 0) return;
 
     if (length < 2) {
-	fprintf (stderr, "FiltFilt::process called for %d samples\n", length);
 	for( i = 0; i < length; i++ ) {
 	    dst[i] = src [i];
 	}
@@ -68,14 +46,13 @@ void FiltFilt::process(double *src, double *dst, unsigned int length)
     unsigned int nFact = 3 * ( nFilt - 1);
     unsigned int nExt	= length + 2 * nFact;
 
-    m_filtScratchIn = new double[ nExt ];
-    m_filtScratchOut = new double[ nExt ];
-
+    double *filtScratchIn = new double[ nExt ];
+    double *filtScratchOut = new double[ nExt ];
 	
     for( i = 0; i< nExt; i++ ) 
     {
-	m_filtScratchIn[ i ] = 0.0;
-	m_filtScratchOut[ i ] = 0.0;
+	filtScratchIn[ i ] = 0.0;
+	filtScratchOut[ i ] = 0.0;
     }
 
     // Edge transients reflection
@@ -85,56 +62,56 @@ void FiltFilt::process(double *src, double *dst, unsigned int length)
     unsigned int index = 0;
     for( i = nFact; i > 0; i-- )
     {
-	m_filtScratchIn[ index++ ] = sample0 - src[ i ];
+	filtScratchIn[ index++ ] = sample0 - src[ i ];
     }
     index = 0;
     for( i = 0; i < nFact && i + 2 < length; i++ )
     {
-	m_filtScratchIn[ (nExt - nFact) + index++ ] = sampleN - src[ (length - 2) - i ];
+	filtScratchIn[ (nExt - nFact) + index++ ] = sampleN - src[ (length - 2) - i ];
     }
 
     for(; i < nFact; i++ )
     {
-	m_filtScratchIn[ (nExt - nFact) + index++ ] = 0;
+	filtScratchIn[ (nExt - nFact) + index++ ] = 0;
     }
 
     index = 0;
     for( i = 0; i < length; i++ )
     {
-	m_filtScratchIn[ i + nFact ] = src[ i ];
+	filtScratchIn[ i + nFact ] = src[ i ];
     }
 	
     ////////////////////////////////
     // Do  0Ph filtering
-    m_filter->process( m_filtScratchIn, m_filtScratchOut, nExt);
+    m_filter.process( filtScratchIn, filtScratchOut, nExt);
 	
     // reverse the series for FILTFILT 
     for ( i = 0; i < nExt; i++)
     { 
-	m_filtScratchIn[ i ] = m_filtScratchOut[ nExt - i - 1];
+	filtScratchIn[ i ] = filtScratchOut[ nExt - i - 1];
     }
 
     // do FILTER again 
-    m_filter->process( m_filtScratchIn, m_filtScratchOut, nExt);
+    m_filter.process( filtScratchIn, filtScratchOut, nExt);
 	
     // reverse the series back 
     for ( i = 0; i < nExt; i++)
     {
-	m_filtScratchIn[ i ] = m_filtScratchOut[ nExt - i - 1 ];
+	filtScratchIn[ i ] = filtScratchOut[ nExt - i - 1 ];
     }
     for ( i = 0;i < nExt; i++)
     {
-	m_filtScratchOut[ i ] = m_filtScratchIn[ i ];
+	filtScratchOut[ i ] = filtScratchIn[ i ];
     }
 
     index = 0;
     for( i = 0; i < length; i++ )
     {
-	dst[ index++ ] = m_filtScratchOut[ i + nFact ];
+	dst[ index++ ] = filtScratchOut[ i + nFact ];
     }	
 
-    delete [] m_filtScratchIn;
-    delete [] m_filtScratchOut;
+    delete [] filtScratchIn;
+    delete [] filtScratchOut;
 
 }
 
diff --git a/libs/qm-dsp/dsp/signalconditioning/FiltFilt.h b/libs/qm-dsp/dsp/signalconditioning/FiltFilt.h
index e5a38124cf..9fc4395068 100644
--- a/libs/qm-dsp/dsp/signalconditioning/FiltFilt.h
+++ b/libs/qm-dsp/dsp/signalconditioning/FiltFilt.h
@@ -20,30 +20,21 @@
 
 /**
  * Zero-phase digital filter, implemented by processing the data
- * through a filter specified by the given FilterConfig structure (see
+ * through a filter specified by the given filter parameters (see
  * Filter) and then processing it again in reverse.
  */
 class FiltFilt  
 {
 public:
-    FiltFilt( FilterConfig Config );
+    FiltFilt(Filter::Parameters);
     virtual ~FiltFilt();
 
     void reset();
     void process( double* src, double* dst, unsigned int length );
 
 private:
-    void initialise( FilterConfig Config );
-    void deInitialise();
-
-    unsigned int m_ord;
-
-    Filter* m_filter;
-
-    double* m_filtScratchIn;
-    double* m_filtScratchOut;
-
-    FilterConfig m_filterConfig;
+    Filter m_filter;
+    int m_ord;
 };
 
 #endif
diff --git a/libs/qm-dsp/dsp/signalconditioning/Filter.cpp b/libs/qm-dsp/dsp/signalconditioning/Filter.cpp
index fcc12e590a..e9523e229d 100644
--- a/libs/qm-dsp/dsp/signalconditioning/Filter.cpp
+++ b/libs/qm-dsp/dsp/signalconditioning/Filter.cpp
@@ -4,7 +4,6 @@
     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
@@ -15,73 +14,112 @@
 
 #include "Filter.h"
 
-//////////////////////////////////////////////////////////////////////
-// Construction/Destruction
-//////////////////////////////////////////////////////////////////////
+#include <stdexcept>
 
-Filter::Filter( FilterConfig Config )
-{
-    m_ord = 0;
-    m_outBuffer = NULL;
-    m_inBuffer = NULL;
+using namespace std;
 
-    initialise( Config );
+Filter::Filter(Parameters params)
+{
+    if (params.a.empty()) {
+        m_fir = true;
+        if (params.b.empty()) {
+            throw logic_error("Filter must have at least one pair of coefficients");
+        }
+    } else {
+        m_fir = false;
+        if (params.a.size() != params.b.size()) {
+            throw logic_error("Inconsistent numbers of filter coefficients");
+        }
+    }
+    
+    m_sz = int(params.b.size());
+    m_order = m_sz - 1;
+
+    m_a = params.a;
+    m_b = params.b;
+    
+    // We keep some empty space at the start of the buffer, and
+    // encroach gradually into it as we add individual sample
+    // calculations at the start. Then when we run out of space, we
+    // move the buffer back to the end and begin again. This is
+    // significantly faster than moving the whole buffer along in
+    // 1-sample steps every time.
+
+    m_offmax = 20;
+    m_offa = m_offmax;
+    m_offb = m_offmax;
+
+    if (!m_fir) {
+        m_bufa.resize(m_order + m_offmax);
+    }
+
+    m_bufb.resize(m_sz + m_offmax);
 }
 
 Filter::~Filter()
 {
-    deInitialise();
 }
 
-void Filter::initialise( FilterConfig Config )
+void
+Filter::reset()
 {
-    m_ord = Config.ord;
-    m_ACoeffs = Config.ACoeffs;
-    m_BCoeffs = Config.BCoeffs;
+    m_offb = m_offmax;
+    m_offa = m_offmax;
 
-    m_inBuffer = new double[ m_ord + 1 ];
-    m_outBuffer = new double[ m_ord + 1 ];
-
-    reset();
-}
+    if (!m_fir) {
+        m_bufa.assign(m_bufa.size(), 0.0);
+    }
 
-void Filter::deInitialise()
-{
-    delete[] m_inBuffer;
-    delete[] m_outBuffer;
+    m_bufb.assign(m_bufb.size(), 0.0);
 }
 
-void Filter::reset()
+void
+Filter::process(const double *const __restrict__ in,
+		double *const __restrict__ out,
+		const int n)
 {
-    for( unsigned int i = 0; i < m_ord+1; i++ ){ m_inBuffer[ i ] = 0.0; }
-    for(unsigned int  i = 0; i < m_ord+1; i++ ){ m_outBuffer[ i ] = 0.0; }
+    for (int s = 0; s < n; ++s) {
+
+        if (m_offb > 0) --m_offb;
+        else {
+            for (int i = m_sz - 2; i >= 0; --i) {
+                m_bufb[i + m_offmax + 1] = m_bufb[i];
+            }
+            m_offb = m_offmax;
+        }
+        m_bufb[m_offb] = in[s];
+
+        double b_sum = 0.0;
+        for (int i = 0; i < m_sz; ++i) {
+            b_sum += m_b[i] * m_bufb[i + m_offb];
+        }
+
+        double outval;
+
+        if (m_fir) {
+
+            outval = b_sum;
+
+        } else {
+
+            double a_sum = 0.0;
+            for (int i = 0; i < m_order; ++i) {
+                a_sum += m_a[i + 1] * m_bufa[i + m_offa];
+            }
+
+            outval = b_sum - a_sum;
+
+            if (m_offa > 0) --m_offa;
+            else {
+                for (int i = m_order - 2; i >= 0; --i) {
+                    m_bufa[i + m_offmax + 1] = m_bufa[i];
+                }
+                m_offa = m_offmax;
+            }
+            m_bufa[m_offa] = outval;
+        }
+        
+	out[s] = outval;
+    }
 }
 
-void Filter::process( double *src, double *dst, unsigned int length )
-{
-    unsigned int SP,i,j;
-
-    double xin,xout;
-
-    for (SP=0;SP<length;SP++)
-    {
-        xin=src[SP];
-        /* move buffer */
-        for ( i = 0; i < m_ord; i++) {m_inBuffer[ m_ord - i ]=m_inBuffer[ m_ord - i - 1 ];}
-        m_inBuffer[0]=xin;
-
-        xout=0.0;
-        for (j=0;j< m_ord + 1; j++)
-	    xout = xout + m_BCoeffs[ j ] * m_inBuffer[ j ];
-        for (j = 0; j < m_ord; j++)
-	    xout= xout - m_ACoeffs[ j + 1 ] * m_outBuffer[ j ];
-
-        dst[ SP ] = xout;
-        for ( i = 0; i < m_ord - 1; i++ ) { m_outBuffer[ m_ord - i - 1 ] = m_outBuffer[ m_ord - i - 2 ];}
-        m_outBuffer[0]=xout;
-
-    } /* end of SP loop */
-}
-
-
-
diff --git a/libs/qm-dsp/dsp/signalconditioning/Filter.h b/libs/qm-dsp/dsp/signalconditioning/Filter.h
index 9f25945fc1..05f79e9723 100644
--- a/libs/qm-dsp/dsp/signalconditioning/Filter.h
+++ b/libs/qm-dsp/dsp/signalconditioning/Filter.h
@@ -4,7 +4,6 @@
     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
@@ -16,46 +15,53 @@
 #ifndef FILTER_H
 #define FILTER_H
 
-#ifndef NULL
-#define NULL 0
-#endif
-
-/**
- * Filter specification. For a filter of order ord, the ACoeffs and
- * BCoeffs arrays must point to ord+1 values each. ACoeffs provides
- * the denominator and BCoeffs the numerator coefficients of the
- * filter.
- */
-struct FilterConfig{
-    unsigned int ord;
-    double* ACoeffs;
-    double* BCoeffs;
-};
+#include <vector>
 
-/**
- * Digital filter specified through FilterConfig structure.
- */
-class Filter  
+class Filter
 {
 public:
-    Filter( FilterConfig Config );
-    virtual ~Filter();
+    struct Parameters {
+        std::vector<double> a;
+        std::vector<double> b;
+    };
+
+    /**
+     * Construct an IIR filter with numerators b and denominators
+     * a. The filter will have order b.size()-1. To make an FIR
+     * filter, leave the vector a in the param struct empty.
+     * Otherwise, a and b must have the same number of values.
+     */
+    Filter(Parameters params);
+    
+    ~Filter();
 
     void reset();
 
-    void process( double *src, double *dst, unsigned int length );
+    /**
+     * Filter the input sequence \arg in of length \arg n samples, and
+     * write the resulting \arg n samples into \arg out. There must be
+     * enough room in \arg out for \arg n samples to be written.
+     */
+    void process(const double *const __restrict__ in,
+                 double *const __restrict__ out,
+                 const int n);
 
+    int getOrder() const { return m_order; }
+    
 private:
-    void initialise( FilterConfig Config );
-    void deInitialise();
+    int m_order;
+    int m_sz;
+    std::vector<double> m_a;
+    std::vector<double> m_b;
+    std::vector<double> m_bufa;
+    std::vector<double> m_bufb;
+    int m_offa;
+    int m_offb;
+    int m_offmax;
+    bool m_fir;
 
-    unsigned int m_ord;
-
-    double* m_inBuffer;
-    double* m_outBuffer;
-
-    double* m_ACoeffs;
-    double* m_BCoeffs;
+    Filter(const Filter &); // not supplied
+    Filter &operator=(const Filter &); // not supplied
 };
-
+    
 #endif
diff --git a/libs/qm-dsp/dsp/tempotracking/TempoTrack.cpp b/libs/qm-dsp/dsp/tempotracking/TempoTrack.cpp
index 389403edaa..7129b37966 100644
--- a/libs/qm-dsp/dsp/tempotracking/TempoTrack.cpp
+++ b/libs/qm-dsp/dsp/tempotracking/TempoTrack.cpp
@@ -70,10 +70,6 @@ void TempoTrack::initialise( TTParams Params )
     m_tempoScratch = new double[ m_lagLength ];
 	m_smoothRCF = new double[ m_lagLength ];
 
-
-    unsigned int winPre = Params.WinT.pre;
-    unsigned int winPost = Params.WinT.post;
-
     m_DFFramer.configure( m_winLength, m_lagLength );
 	
     m_DFPParams.length = m_winLength;
@@ -120,9 +116,9 @@ void TempoTrack::deInitialise()
 
 }
 
-void TempoTrack::createCombFilter(double* Filter, unsigned int winLength, unsigned int TSig, double beatLag)
+void TempoTrack::createCombFilter(double* Filter, int winLength, int /* TSig */, double beatLag)
 {
-    unsigned int i;
+    int i;
 
     if( beatLag == 0 )
     {
@@ -147,15 +143,15 @@ double TempoTrack::tempoMM(double* ACF, double* weight, int tsig)
 
     double period = 0;
     double maxValRCF = 0.0;
-    unsigned int maxIndexRCF = 0;
+    int maxIndexRCF = 0;
 
     double* pdPeaks;
 
-    unsigned int maxIndexTemp;
-    double	maxValTemp;
-    unsigned int count; 
+    int maxIndexTemp;
+    double maxValTemp;
+    int count; 
 	
-    unsigned int numelem,i,j;
+    int numelem,i,j;
     int a, b;
 
     for( i = 0; i < m_lagLength; i++ )
@@ -476,7 +472,7 @@ void TempoTrack::constDetect( double* periodP, int currentIdx, int* flag )
     }
 }
 
-int TempoTrack::findMeter(double *ACF, unsigned int len, double period)
+int TempoTrack::findMeter(double *ACF, int len, double period)
 {
     int i;
     int p = (int)MathUtilities::round( period );
@@ -491,7 +487,7 @@ int TempoTrack::findMeter(double *ACF, unsigned int len, double period)
     double temp4B = 0.0;
 
     double* dbf = new double[ len ]; int t = 0;
-    for( unsigned int u = 0; u < len; u++ ){ dbf[ u ] = 0.0; }
+    for( int u = 0; u < len; u++ ){ dbf[ u ] = 0.0; }
 
     if( (double)len < 6 * p + 2 )
     {
@@ -548,7 +544,7 @@ int TempoTrack::findMeter(double *ACF, unsigned int len, double period)
     return tsig;
 }
 
-void TempoTrack::createPhaseExtractor(double *Filter, unsigned int winLength, double period, unsigned int fsp, unsigned int lastBeat)
+void TempoTrack::createPhaseExtractor(double *Filter, int /* winLength */, double period, int fsp, int lastBeat)
 {	
     int p = (int)MathUtilities::round( period );
     int predictedOffset = 0;
@@ -584,7 +580,7 @@ void TempoTrack::createPhaseExtractor(double *Filter, unsigned int winLength, do
 	double sigma = (double)p/8;
 	double PhaseMin = 0.0;
 	double PhaseMax = 0.0;
-	unsigned int scratchLength = p*2;
+	int scratchLength = p*2;
 	double temp = 0.0;
 
 	for(  int i = 0; i < scratchLength; i++ )
@@ -604,7 +600,7 @@ void TempoTrack::createPhaseExtractor(double *Filter, unsigned int winLength, do
         std::cerr << "predictedOffset = " << predictedOffset << std::endl;
 #endif
 
-	unsigned int index = 0;
+	int index = 0;
 	for (int i = p - ( predictedOffset - 1); i < p + ( p - predictedOffset) + 1; i++)
 	{
 #ifdef DEBUG_TEMPO_TRACK
@@ -624,7 +620,7 @@ void TempoTrack::createPhaseExtractor(double *Filter, unsigned int winLength, do
     delete [] phaseScratch;
 }
 
-int TempoTrack::phaseMM(double *DF, double *weighting, unsigned int winLength, double period)
+int TempoTrack::phaseMM(double *DF, double *weighting, int winLength, double period)
 {
     int alignment = 0;
     int p = (int)MathUtilities::round( period );
@@ -667,7 +663,7 @@ int TempoTrack::phaseMM(double *DF, double *weighting, unsigned int winLength, d
     return alignment;
 }
 
-int TempoTrack::beatPredict(unsigned int FSP0, double alignment, double period, unsigned int step )
+int TempoTrack::beatPredict(int FSP0, double alignment, double period, int step )
 {
     int beat = 0;
 
@@ -712,39 +708,39 @@ vector<int> TempoTrack::process( vector <double> DF,
     causalDF = DF;
 
     //Prepare Causal Extension DFData
-    unsigned int DFCLength = m_dataLength + m_winLength;
+//    int DFCLength = m_dataLength + m_winLength;
 	
-    for( unsigned int j = 0; j < m_winLength; j++ )
+    for( int j = 0; j < m_winLength; j++ )
     {
 	causalDF.push_back( 0 );
     }
 	
 	
     double* RW = new double[ m_lagLength ];
-    for( unsigned int clear = 0; clear < m_lagLength; clear++){ RW[ clear ] = 0.0;}
+    for (int clear = 0; clear < m_lagLength; clear++){ RW[ clear ] = 0.0;}
 
     double* GW = new double[ m_lagLength ];
-    for(unsigned int clear = 0; clear < m_lagLength; clear++){ GW[ clear ] = 0.0;}
+    for (int clear = 0; clear < m_lagLength; clear++){ GW[ clear ] = 0.0;}
 
     double* PW = new double[ m_lagLength ];
-    for(unsigned clear = 0; clear < m_lagLength; clear++){ PW[ clear ] = 0.0;}
+    for(int clear = 0; clear < m_lagLength; clear++){ PW[ clear ] = 0.0;}
 
     m_DFFramer.setSource( &causalDF[0], m_dataLength );
 
-    unsigned int TTFrames = m_DFFramer.getMaxNoFrames();
+    int TTFrames = m_DFFramer.getMaxNoFrames();
 
 #ifdef DEBUG_TEMPO_TRACK
     std::cerr << "TTFrames = " << TTFrames << std::endl;
 #endif
 	
     double* periodP = new double[ TTFrames ];
-    for(unsigned clear = 0; clear < TTFrames; clear++){ periodP[ clear ] = 0.0;}
+    for(int clear = 0; clear < TTFrames; clear++){ periodP[ clear ] = 0.0;}
 	
     double* periodG = new double[ TTFrames ];
-    for(unsigned clear = 0; clear < TTFrames; clear++){ periodG[ clear ] = 0.0;}
+    for(int clear = 0; clear < TTFrames; clear++){ periodG[ clear ] = 0.0;}
 	
     double* alignment = new double[ TTFrames ];
-    for(unsigned clear = 0; clear < TTFrames; clear++){ alignment[ clear ] = 0.0;}
+    for(int clear = 0; clear < TTFrames; clear++){ alignment[ clear ] = 0.0;}
 
     m_beats.clear();
 
@@ -752,7 +748,7 @@ vector<int> TempoTrack::process( vector <double> DF,
 
     int TTLoopIndex = 0;
 
-    for( unsigned int i = 0; i < TTFrames; i++ )
+    for( int i = 0; i < TTFrames; i++ )
     {
 	m_DFFramer.getFrame( m_rawDFFrame );
 
diff --git a/libs/qm-dsp/dsp/tempotracking/TempoTrack.h b/libs/qm-dsp/dsp/tempotracking/TempoTrack.h
index 0c315717ba..72c2f756ef 100644
--- a/libs/qm-dsp/dsp/tempotracking/TempoTrack.h
+++ b/libs/qm-dsp/dsp/tempotracking/TempoTrack.h
@@ -30,16 +30,16 @@ using std::vector;
 

 struct WinThresh

 {

-    unsigned int pre;

-    unsigned int post;

+    int pre;

+    int post;

 };

 

 struct TTParams

 {

-    unsigned int winLength; //Analysis window length

-    unsigned int lagLength; //Lag & Stride size

-    unsigned int alpha; //alpha-norm parameter

-    unsigned int LPOrd; // low-pass Filter order

+    int winLength; //Analysis window length

+    int lagLength; //Lag & Stride size

+    int alpha; //alpha-norm parameter

+    int LPOrd; // low-pass Filter order

     double* LPACoeffs; //low pass Filter den coefficients

     double* LPBCoeffs; //low pass Filter num coefficients

     WinThresh WinT;//window size in frames for adaptive thresholding [pre post]:

@@ -59,22 +59,22 @@ private:
     void initialise( TTParams Params );

     void deInitialise();

 

-    int beatPredict( unsigned int FSP, double alignment, double period, unsigned int step);

-    int phaseMM( double* DF, double* weighting, unsigned int winLength, double period );

-    void createPhaseExtractor( double* Filter, unsigned int winLength,  double period,  unsigned int fsp, unsigned int lastBeat );

-    int findMeter( double* ACF,  unsigned int len, double period );

+    int beatPredict( int FSP, double alignment, double period, int step);

+    int phaseMM( double* DF, double* weighting, int winLength, double period );

+    void createPhaseExtractor( double* Filter, int winLength,  double period,  int fsp, int lastBeat );

+    int findMeter( double* ACF,  int len, double period );

     void constDetect( double* periodP, int currentIdx, int* flag );

     void stepDetect( double* periodP, double* periodG, int currentIdx, int* flag );

-    void createCombFilter( double* Filter, unsigned int winLength, unsigned int TSig, double beatLag );

+    void createCombFilter( double* Filter, int winLength, int TSig, double beatLag );

     double tempoMM( double* ACF, double* weight, int sig );

 	

-    unsigned int m_dataLength;

-    unsigned int m_winLength;

-    unsigned int m_lagLength;

+    int m_dataLength;

+    int m_winLength;

+    int m_lagLength;

 

-    double		 m_rayparam;

-    double		 m_sigma;

-    double		 m_DFWVNnorm;

+    double m_rayparam;

+    double m_sigma;

+    double m_DFWVNnorm;

 

     vector<int>	 m_beats; // Vector of detected beats

 

diff --git a/libs/qm-dsp/dsp/tonal/TCSgram.cpp b/libs/qm-dsp/dsp/tonal/TCSgram.cpp
index 8840872576..8730dcf110 100644
--- a/libs/qm-dsp/dsp/tonal/TCSgram.cpp
+++ b/libs/qm-dsp/dsp/tonal/TCSgram.cpp
@@ -36,7 +36,7 @@ void TCSGram::getTCSVector(int iPosition, TCSVector& rTCSVector) const
 {
 	if (iPosition < 0) 
 		rTCSVector = TCSVector();
-	else if (iPosition >= m_VectorList.size())
+	else if (iPosition >= int(m_VectorList.size()))
 		rTCSVector = TCSVector();
 	else
 		rTCSVector = m_VectorList[iPosition].second;
diff --git a/libs/qm-dsp/dsp/tonal/TonalEstimator.h b/libs/qm-dsp/dsp/tonal/TonalEstimator.h
index f555ef9e3b..b3c9c7b235 100644
--- a/libs/qm-dsp/dsp/tonal/TonalEstimator.h
+++ b/libs/qm-dsp/dsp/tonal/TonalEstimator.h
@@ -32,7 +32,7 @@ public:
 	
 	void printDebug()
 	{
-		for (int i = 0; i < size(); i++)
+		for (int i = 0; i < int(size()); i++)
 		{
 			std::cout <<  (*this)[i] << ";";
 		}
@@ -68,7 +68,7 @@ public:
 
 	void printDebug()
 	{
-		for (int i = 0; i < size(); i++)
+		for (int i = 0; i < int(size()); i++)
 		{
 			std::cout <<  (*this)[i] << ";";
 		}
diff --git a/libs/qm-dsp/dsp/transforms/FFT.h b/libs/qm-dsp/dsp/transforms/FFT.h
index c8956ac7f2..73af5616e7 100644
--- a/libs/qm-dsp/dsp/transforms/FFT.h
+++ b/libs/qm-dsp/dsp/transforms/FFT.h
@@ -4,6 +4,12 @@
     QM DSP Library
 
     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.
 */
 
 #ifndef FFT_H
diff --git a/libs/qm-dsp/dsp/wavelet/Wavelet.cpp b/libs/qm-dsp/dsp/wavelet/Wavelet.cpp
index 764c84b24a..9c29b73cf5 100644
--- a/libs/qm-dsp/dsp/wavelet/Wavelet.cpp
+++ b/libs/qm-dsp/dsp/wavelet/Wavelet.cpp
@@ -1843,7 +1843,10 @@ Wavelet::createDecompositionFilters(Type wavelet,
         break;
     }
 
-    assert(flength == lpd.size());
-    assert(flength == hpd.size());
+    // avoid compiler warning for unused value if assert is not compiled in:
+    (void)flength;
+
+    assert(flength == int(lpd.size()));
+    assert(flength == int(hpd.size()));
 }
 
diff --git a/libs/qm-dsp/ext/kissfft/kiss_fft.c b/libs/qm-dsp/ext/kissfft/kiss_fft.c
index 465d6c97a0..7824d34521 100644
--- a/libs/qm-dsp/ext/kissfft/kiss_fft.c
+++ b/libs/qm-dsp/ext/kissfft/kiss_fft.c
@@ -275,7 +275,8 @@ void kf_work(
 
     if (m==1) {
         do{
-            *Fout = *f;
+            Fout->r = f->r;
+            Fout->i = f->i;
             f += fstride*in_stride;
         }while(++Fout != Fout_end );
     }else{
diff --git a/libs/qm-dsp/gitrev.txt b/libs/qm-dsp/gitrev.txt
new file mode 100644
index 0000000000..dde4075d05
--- /dev/null
+++ b/libs/qm-dsp/gitrev.txt
@@ -0,0 +1 @@
+qm-vamp-plugins-v1.7.1-20-g4d15479
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