changeset 483:fdaa63607c15

Untabify, indent, tidy
author Chris Cannam <cannam@all-day-breakfast.com>
date Fri, 31 May 2019 11:54:32 +0100
parents cbe668c7d724
children d48276a3ae24
files base/KaiserWindow.cpp base/KaiserWindow.h base/Pitch.cpp base/Pitch.h base/SincWindow.cpp base/SincWindow.h base/Window.h dsp/chromagram/Chromagram.cpp dsp/chromagram/Chromagram.h dsp/chromagram/ConstantQ.cpp dsp/chromagram/ConstantQ.h dsp/onsets/DetectionFunction.cpp dsp/onsets/DetectionFunction.h dsp/onsets/PeakPicking.cpp dsp/onsets/PeakPicking.h dsp/phasevocoder/PhaseVocoder.cpp dsp/rateconversion/Decimator.cpp dsp/rateconversion/Decimator.h dsp/rateconversion/DecimatorB.cpp dsp/rateconversion/DecimatorB.h dsp/rateconversion/Resampler.cpp dsp/signalconditioning/DFProcess.cpp dsp/signalconditioning/DFProcess.h dsp/signalconditioning/FiltFilt.cpp dsp/signalconditioning/Filter.cpp dsp/signalconditioning/Framer.cpp dsp/signalconditioning/Framer.h dsp/transforms/DCT.cpp dsp/wavelet/Wavelet.cpp dsp/wavelet/Wavelet.h hmm/hmm.c hmm/hmm.h maths/pca/pca.c
diffstat 33 files changed, 1427 insertions(+), 1616 deletions(-) [+]
line wrap: on
line diff
--- a/base/KaiserWindow.cpp	Fri May 31 11:02:28 2019 +0100
+++ b/base/KaiserWindow.cpp	Fri May 31 11:54:32 2019 +0100
@@ -17,27 +17,27 @@
 
 KaiserWindow::Parameters
 KaiserWindow::parametersForTransitionWidth(double attenuation,
-					   double transition)
+                                           double transition)
 {
     Parameters p;
     p.length = 1 + (attenuation > 21.0 ?
-		    ceil((attenuation - 7.95) / (2.285 * transition)) :
-		    ceil(5.79 / transition));
+                    ceil((attenuation - 7.95) / (2.285 * transition)) :
+                    ceil(5.79 / transition));
     p.beta = (attenuation > 50.0 ? 
-	      0.1102 * (attenuation - 8.7) :
-	      attenuation > 21.0 ? 
-	      0.5842 * pow(attenuation - 21.0, 0.4) + 0.07886 * (attenuation - 21.0) :
-	      0);
+              0.1102 * (attenuation - 8.7) :
+              attenuation > 21.0 ? 
+              0.5842 * pow(attenuation - 21.0, 0.4) + 0.07886 * (attenuation - 21.0) :
+              0);
     return p;
 }
 
 static double besselTerm(double x, int i)
 {
     if (i == 0) {
-	return 1;
+        return 1;
     } else {
-	double f = MathUtilities::factorial(i);
-	return pow(x/2, i*2) / (f*f);
+        double f = MathUtilities::factorial(i);
+        return pow(x/2, i*2) / (f*f);
     }
 }
 
@@ -45,7 +45,7 @@
 {
     double b = 0.0;
     for (int i = 0; i < 20; ++i) {
-	b += besselTerm(x, i);
+        b += besselTerm(x, i);
     }
     return b;
 }
@@ -56,8 +56,8 @@
     double denominator = bessel0(m_beta);
     bool even = (m_length % 2 == 0);
     for (int i = 0; i < (even ? m_length/2 : (m_length+1)/2); ++i) {
-	double k = double(2*i) / double(m_length-1) - 1.0;
-	m_window.push_back(bessel0(m_beta * sqrt(1.0 - k*k)) / denominator);
+        double k = double(2*i) / double(m_length-1) - 1.0;
+        m_window.push_back(bessel0(m_beta * sqrt(1.0 - k*k)) / denominator);
     }
     for (int i = 0; i < (even ? m_length/2 : (m_length-1)/2); ++i) {
         m_window.push_back(m_window[int(m_length/2) - i - 1]);
--- a/base/KaiserWindow.h	Fri May 31 11:02:28 2019 +0100
+++ b/base/KaiserWindow.h	Fri May 31 11:54:32 2019 +0100
@@ -26,8 +26,8 @@
 {
 public:
     struct Parameters {
-	int length;
-	double beta;
+        int length;
+        double beta;
     };
 
     /**
@@ -41,9 +41,9 @@
      * and transition width in samples.
      */
     static KaiserWindow byTransitionWidth(double attenuation,
-					  double transition) {
-	return KaiserWindow
-	    (parametersForTransitionWidth(attenuation, transition));
+                                          double transition) {
+        return KaiserWindow
+            (parametersForTransitionWidth(attenuation, transition));
     }
 
     /**
@@ -51,10 +51,10 @@
      * and transition bandwidth in Hz for the given samplerate.
      */
     static KaiserWindow byBandwidth(double attenuation,
-				    double bandwidth,
-				    double samplerate) {
-	return KaiserWindow
-	    (parametersForBandwidth(attenuation, bandwidth, samplerate));
+                                    double bandwidth,
+                                    double samplerate) {
+        return KaiserWindow
+            (parametersForBandwidth(attenuation, bandwidth, samplerate));
     }
 
     /**
@@ -62,7 +62,7 @@
      * given attenuation in dB and transition width in samples.
      */
     static Parameters parametersForTransitionWidth(double attenuation,
-						   double transition);
+                                                   double transition);
 
     /**
      * Obtain the parameters necessary for a Kaiser window of the
@@ -70,28 +70,28 @@
      * given samplerate.
      */
     static Parameters parametersForBandwidth(double attenuation,
-					     double bandwidth,
-					     double samplerate) {
-	return parametersForTransitionWidth
-	    (attenuation, (bandwidth * 2 * M_PI) / samplerate);
+                                             double bandwidth,
+                                             double samplerate) {
+        return parametersForTransitionWidth
+            (attenuation, (bandwidth * 2 * M_PI) / samplerate);
     } 
 
     int getLength() const {
-	return m_length;
+        return m_length;
     }
 
     const double *getWindow() const { 
-	return m_window.data();
+        return m_window.data();
     }
 
     void cut(double *src) const { 
-	cut(src, src); 
+        cut(src, src); 
     }
 
     void cut(const double *src, double *dst) const {
-	for (int i = 0; i < m_length; ++i) {
-	    dst[i] = src[i] * m_window[i];
-	}
+        for (int i = 0; i < m_length; ++i) {
+            dst[i] = src[i] * m_window[i];
+        }
     }
 
 private:
--- a/base/Pitch.cpp	Fri May 31 11:02:28 2019 +0100
+++ b/base/Pitch.cpp	Fri May 31 11:54:32 2019 +0100
@@ -18,8 +18,8 @@
 
 float
 Pitch::getFrequencyForPitch(int midiPitch,
-			    float centsOffset,
-			    float concertA)
+                            float centsOffset,
+                            float concertA)
 {
     float p = float(midiPitch) + (centsOffset / 100);
     return concertA * powf(2.0, (p - 69.0) / 12.0);
@@ -27,8 +27,8 @@
 
 int
 Pitch::getPitchForFrequency(float frequency,
-			    float *centsOffsetReturn,
-			    float concertA)
+                            float *centsOffsetReturn,
+                            float concertA)
 {
     float p = 12.0 * (log(frequency / (concertA / 2.0)) / log(2.0)) + 57.0;
 
@@ -36,8 +36,8 @@
     float centsOffset = (p - midiPitch) * 100.0;
 
     if (centsOffset >= 50.0) {
-	midiPitch = midiPitch + 1;
-	centsOffset = -(100.0 - centsOffset);
+        midiPitch = midiPitch + 1;
+        centsOffset = -(100.0 - centsOffset);
     }
     
     if (centsOffsetReturn) *centsOffsetReturn = centsOffset;
--- a/base/Pitch.h	Fri May 31 11:02:28 2019 +0100
+++ b/base/Pitch.h	Fri May 31 11:54:32 2019 +0100
@@ -23,12 +23,12 @@
 {
 public:
     static float getFrequencyForPitch(int midiPitch,
-				      float centsOffset = 0,
-				      float concertA = 440.0);
+                                      float centsOffset = 0,
+                                      float concertA = 440.0);
 
     static int getPitchForFrequency(float frequency,
-				    float *centsOffsetReturn = 0,
-				    float concertA = 440.0);
+                                    float *centsOffsetReturn = 0,
+                                    float concertA = 440.0);
 };
 
 
--- a/base/SincWindow.cpp	Fri May 31 11:02:28 2019 +0100
+++ b/base/SincWindow.cpp	Fri May 31 11:54:32 2019 +0100
@@ -19,27 +19,27 @@
 SincWindow::init()
 {
     if (m_length < 1) {
-	return;
+        return;
     } else if (m_length < 2) {
-	m_window.push_back(1);
-	return;
+        m_window.push_back(1);
+        return;
     } else {
 
-	int n0 = (m_length % 2 == 0 ? m_length/2 : (m_length - 1)/2);
-	int n1 = (m_length % 2 == 0 ? m_length/2 : (m_length + 1)/2);
-	double m = 2 * M_PI / m_p;
+        int n0 = (m_length % 2 == 0 ? m_length/2 : (m_length - 1)/2);
+        int n1 = (m_length % 2 == 0 ? m_length/2 : (m_length + 1)/2);
+        double m = 2 * M_PI / m_p;
 
-	for (int i = 0; i < n0; ++i) {
-	    double x = ((m_length / 2) - i) * m;
-	    m_window.push_back(sin(x) / x);
-	}
+        for (int i = 0; i < n0; ++i) {
+            double x = ((m_length / 2) - i) * m;
+            m_window.push_back(sin(x) / x);
+        }
 
-	m_window.push_back(1.0);
+        m_window.push_back(1.0);
 
-	for (int i = 1; i < n1; ++i) {
-	    double x = i * m;
-	    m_window.push_back(sin(x) / x);
-	}
+        for (int i = 1; i < n1; ++i) {
+            double x = i * m;
+            m_window.push_back(sin(x) / x);
+        }
     }
 }
 
--- a/base/SincWindow.h	Fri May 31 11:02:28 2019 +0100
+++ b/base/SincWindow.h	Fri May 31 11:54:32 2019 +0100
@@ -33,21 +33,21 @@
     SincWindow(int length, double p) : m_length(length), m_p(p) { init(); }
 
     int getLength() const {
-	return m_length;
+        return m_length;
     }
 
     const double *getWindow() const { 
-	return m_window.data();
+        return m_window.data();
     }
 
     void cut(double *src) const { 
-	cut(src, src); 
+        cut(src, src); 
     }
 
     void cut(const double *src, double *dst) const {
-	for (int i = 0; i < m_length; ++i) {
-	    dst[i] = src[i] * m_window[i];
-	}
+        for (int i = 0; i < m_length; ++i) {
+            dst[i] = src[i] * m_window[i];
+        }
     }
 
 private:
--- a/base/Window.h	Fri May 31 11:02:28 2019 +0100
+++ b/base/Window.h	Fri May 31 11:54:32 2019 +0100
@@ -50,17 +50,17 @@
     Window(WindowType type, int size) : m_type(type), m_size(size) { encache(); }
     Window(const Window &w) : m_type(w.m_type), m_size(w.m_size) { encache(); }
     Window &operator=(const Window &w) {
-	if (&w == this) return *this;
-	m_type = w.m_type;
-	m_size = w.m_size;
-	encache();
-	return *this;
+        if (&w == this) return *this;
+        m_type = w.m_type;
+        m_size = w.m_size;
+        encache();
+        return *this;
     }
     virtual ~Window() { delete[] m_cache; }
     
     void cut(T *src) const { cut(src, src); }
     void cut(const T *src, T *dst) const {
-	for (int i = 0; i < m_size; ++i) dst[i] = src[i] * m_cache[i];
+        for (int i = 0; i < m_size; ++i) dst[i] = src[i] * m_cache[i];
     }
 
     WindowType getType() const { return m_type; }
@@ -91,13 +91,13 @@
     for (i = 0; i < n; ++i) mult[i] = 1.0;
 
     switch (m_type) {
-		
+                
     case RectangularWindow:
         for (i = 0; i < n; ++i) {
             mult[i] = mult[i] * 0.5;
-	}
-	break;
-	    
+        }
+        break;
+            
     case BartlettWindow:
         if (n == 2) {
             mult[0] = mult[1] = 0; // "matlab compatible"
@@ -109,34 +109,34 @@
                 mult[i] = mult[i] * (i / T(n/2));
                 mult[i + n - n/2] = mult[i + n - n/2] * (1.0 - (i / T(n/2)));
             }
-	}
-	break;
-	    
+        }
+        break;
+            
     case HammingWindow:
         if (n > 1) {
             for (i = 0; i < n; ++i) {
                 mult[i] = mult[i] * (0.54 - 0.46 * cos(2 * M_PI * i / n));
             }
-	}
-	break;
-	    
+        }
+        break;
+            
     case HanningWindow:
         if (n > 1) {
             for (i = 0; i < n; ++i) {
                 mult[i] = mult[i] * (0.50 - 0.50 * cos(2 * M_PI * i / n));
             }
-	}
-	break;
-	    
+        }
+        break;
+            
     case BlackmanWindow:
         if (n > 1) {
             for (i = 0; i < n; ++i) {
                 mult[i] = mult[i] * (0.42 - 0.50 * cos(2 * M_PI * i / n)
                                      + 0.08 * cos(4 * M_PI * i / n));
             }
-	}
-	break;
-	    
+        }
+        break;
+            
     case BlackmanHarrisWindow:
         if (n > 1) {
             for (i = 0; i < n; ++i) {
@@ -145,10 +145,10 @@
                                      + 0.14128 * cos(4 * M_PI * i / n)
                                      - 0.01168 * cos(6 * M_PI * i / n));
             }
-	}
-	break;
+        }
+        break;
     }
-	   
+           
     m_cache = mult;
 }
 
--- a/dsp/chromagram/Chromagram.cpp	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/chromagram/Chromagram.cpp	Fri May 31 11:54:32 2019 +0100
@@ -1,5 +1,4 @@
 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
-
 /*
     QM DSP Library
 
@@ -27,10 +26,10 @@
 }
 
 int Chromagram::initialise( ChromaConfig Config )
-{	
-    m_FMin = Config.min;		// min freq
-    m_FMax = Config.max;		// max freq
-    m_BPO  = Config.BPO;		// bins per octave
+{       
+    m_FMin = Config.min;                // min freq
+    m_FMax = Config.max;                // max freq
+    m_BPO  = Config.BPO;                // bins per octave
     m_normalise = Config.normalise;     // if frame normalisation is required
 
     // Extend range to a full octave
@@ -50,7 +49,7 @@
     ConstantQConfig.max = m_FMax;
     ConstantQConfig.BPO = m_BPO;
     ConstantQConfig.CQThresh = Config.CQThresh;
-	
+        
     // Initialise ConstantQ operator
     m_ConstantQ = new ConstantQ( ConstantQConfig );
 
@@ -61,7 +60,7 @@
     m_frameSize = m_ConstantQ->getfftlength();
     m_hopSize = m_ConstantQ->gethop();
 
-    // Initialise FFT object	
+    // Initialise FFT object    
     m_FFT = new FFTReal(m_frameSize);
 
     m_FFTRe = new double[ m_frameSize ];
@@ -111,16 +110,13 @@
 void Chromagram::unityNormalise(double *src)
 {
     double min, max;
-
     double val = 0;
 
     MathUtilities::getFrameMinMax( src, m_BPO, & min, &max );
 
-    for (int i = 0; i < m_BPO; i++)
-    {
-	val = src[ i ] / max;
-
-	src[ i ] = val;
+    for (int i = 0; i < m_BPO; i++) {
+        val = src[ i ] / max;
+        src[ i ] = val;
     }
 }
 
@@ -169,16 +165,15 @@
 
     // Calculate ConstantQ frame
     m_ConstantQ->process( real, imag, m_CQRe, m_CQIm );
-	
+        
     // add each octave of cq data into Chromagram
     const int octaves = m_uK / m_BPO;
-    for (int octave = 0; octave < octaves; octave++) 
-    {
-	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 ]);
-	}
+    for (int octave = 0; octave < octaves; octave++) {
+        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 ]);
+        }
     }
 
     MathUtilities::normalise(m_chromadata, m_BPO, m_normalise);
--- a/dsp/chromagram/Chromagram.h	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/chromagram/Chromagram.h	Fri May 31 11:54:32 2019 +0100
@@ -1,5 +1,4 @@
 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
-
 /*
     QM DSP Library
 
@@ -31,8 +30,7 @@
 
 class Chromagram 
 {
-
-public:	
+public: 
     Chromagram( ChromaConfig Config );
     ~Chromagram();
 
@@ -67,7 +65,7 @@
 
     // Complex arithmetic
     double kabs( double real, double imag );
-	
+        
     // Results
     int getK() { return m_uK;}
     int getFrameSize() { return m_frameSize; }
@@ -79,7 +77,7 @@
 
     Window<double> *m_window;
     double *m_windowbuf;
-	
+        
     double* m_chromadata;
     double m_FMin;
     double m_FMax;
--- a/dsp/chromagram/ConstantQ.cpp	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/chromagram/ConstantQ.cpp	Fri May 31 11:54:32 2019 +0100
@@ -56,10 +56,9 @@
     double* transfHammingWindowRe = new double [ m_FFTLength ];
     double* transfHammingWindowIm = new double [ m_FFTLength ];
 
-    for (unsigned u=0; u < m_FFTLength; u++) 
-    {
-	hammingWindowRe[u] = 0;
-	hammingWindowIm[u] = 0;
+    for (unsigned u=0; u < m_FFTLength; u++) {
+        hammingWindowRe[u] = 0;
+        hammingWindowIm[u] = 0;
     }
 
     // Here, fftleng*2 is a guess of the number of sparse cells in the matrix
@@ -69,39 +68,36 @@
     sk->js.reserve( m_FFTLength*2 );
     sk->real.reserve( m_FFTLength*2 );
     sk->imag.reserve( m_FFTLength*2 );
-	
+        
     // for each bin value K, calculate temporal kernel, take its fft to
     //calculate the spectral kernel then threshold it to make it sparse and 
     //add it to the sparse kernels matrix
     double squareThreshold = m_CQThresh * m_CQThresh;
 
     FFT m_FFT(m_FFTLength);
-	
-    for (unsigned k = m_uK; k--; ) 
-    {
-        for (unsigned u=0; u < m_FFTLength; u++) 
-        {
+        
+    for (unsigned k = m_uK; k--; ) {
+        for (unsigned u=0; u < m_FFTLength; u++) {
             hammingWindowRe[u] = 0;
             hammingWindowIm[u] = 0;
         }
         
-	// Computing a hamming window
-	const unsigned hammingLength = (int) ceil( m_dQ * m_FS / ( m_FMin * pow(2,((double)(k))/(double)m_BPO)));
+        // Computing a hamming window
+        const unsigned hammingLength = (int) ceil( m_dQ * m_FS / ( m_FMin * pow(2,((double)(k))/(double)m_BPO)));
 
 //        cerr << "k = " << k << ", q = " << m_dQ << ", m_FMin = " << m_FMin << ", hammingLength = " << hammingLength << " (rounded up from " << (m_dQ * m_FS / ( m_FMin * pow(2,((double)(k))/(double)m_BPO))) << ")" << endl;
         
 
         unsigned origin = m_FFTLength/2 - hammingLength/2;
 
-	for (unsigned i=0; i<hammingLength; i++) 
-	{
-	    const double angle = 2*PI*m_dQ*i/hammingLength;
-	    const double real = cos(angle);
-	    const double imag = sin(angle);
-	    const double absol = hamming(hammingLength, i)/hammingLength;
-	    hammingWindowRe[ origin + i ] = absol*real;
-	    hammingWindowIm[ origin + i ] = absol*imag;
-	}
+        for (unsigned i=0; i<hammingLength; i++) {
+            const double angle = 2*PI*m_dQ*i/hammingLength;
+            const double real = cos(angle);
+            const double imag = sin(angle);
+            const double absol = hamming(hammingLength, i)/hammingLength;
+            hammingWindowRe[ origin + i ] = absol*real;
+            hammingWindowIm[ origin + i ] = absol*imag;
+        }
 
         for (unsigned i = 0; i < m_FFTLength/2; ++i) {
             double temp = hammingWindowRe[i];
@@ -112,25 +108,22 @@
             hammingWindowIm[i + m_FFTLength/2] = temp;
         }
     
-	//do fft of hammingWindow
-	m_FFT.process( 0, hammingWindowRe, hammingWindowIm, transfHammingWindowRe, transfHammingWindowIm );
+        //do fft of hammingWindow
+        m_FFT.process( 0, hammingWindowRe, hammingWindowIm, transfHammingWindowRe, transfHammingWindowIm );
+                
+        for (unsigned j=0; j<( m_FFTLength ); j++) {
+            // perform thresholding
+            const double squaredBin = squaredModule( transfHammingWindowRe[ j ], transfHammingWindowIm[ j ]);
+            if (squaredBin <= squareThreshold) continue;
+                
+            // Insert non-zero position indexes
+            sk->is.push_back(j);
+            sk->js.push_back(k);
 
-		
-	for (unsigned j=0; j<( m_FFTLength ); j++) 
-	{
-	    // perform thresholding
-	    const double squaredBin = squaredModule( transfHammingWindowRe[ j ], transfHammingWindowIm[ j ]);
-	    if (squaredBin <= squareThreshold) continue;
-		
-	    // Insert non-zero position indexes
-	    sk->is.push_back(j);
-	    sk->js.push_back(k);
-
-	    // take conjugate, normalise and add to array sparkernel
-	    sk->real.push_back( transfHammingWindowRe[ j ]/m_FFTLength);
-	    sk->imag.push_back(-transfHammingWindowIm[ j ]/m_FFTLength);
-	}
-
+            // take conjugate, normalise and add to array sparkernel
+            sk->real.push_back( transfHammingWindowRe[ j ]/m_FFTLength);
+            sk->imag.push_back(-transfHammingWindowIm[ j ]/m_FFTLength);
+        }
     }
 
     delete [] hammingWindowRe;
@@ -154,29 +147,27 @@
 
     SparseKernel *sk = m_sparseKernel;
 
-    for (unsigned row=0; row<2*m_uK; row++) 
-    {
-	m_CQdata[ row ] = 0;
-	m_CQdata[ row+1 ] = 0;
+    for (unsigned row=0; row<2*m_uK; row++) {
+        m_CQdata[ row ] = 0;
+        m_CQdata[ row+1 ] = 0;
     }
     const unsigned *fftbin = &(sk->is[0]);
     const unsigned *cqbin  = &(sk->js[0]);
     const double   *real   = &(sk->real[0]);
     const double   *imag   = &(sk->imag[0]);
     const unsigned int sparseCells = sk->real.size();
-	
-    for (unsigned i = 0; i<sparseCells; i++)
-    {
-	const unsigned row = cqbin[i];
-	const unsigned col = fftbin[i];
+        
+    for (unsigned i = 0; i<sparseCells; i++) {
+        const unsigned row = cqbin[i];
+        const unsigned col = fftbin[i];
         if (col == 0) continue;
-	const double & r1  = real[i];
-	const double & i1  = imag[i];
-	const double & r2  = fftdata[ (2*m_FFTLength) - 2*col - 2 ];
-	const double & i2  = fftdata[ (2*m_FFTLength) - 2*col - 2 + 1 ];
-	// add the multiplication
-	m_CQdata[ 2*row  ] += (r1*r2 - i1*i2);
-	m_CQdata[ 2*row+1] += (r1*i2 + i1*r2);
+        const double & r1  = real[i];
+        const double & i1  = imag[i];
+        const double & r2  = fftdata[ (2*m_FFTLength) - 2*col - 2 ];
+        const double & i2  = fftdata[ (2*m_FFTLength) - 2*col - 2 + 1 ];
+        // add the multiplication
+        m_CQdata[ 2*row  ] += (r1*r2 - i1*i2);
+        m_CQdata[ 2*row+1] += (r1*i2 + i1*r2);
     }
 
     return m_CQdata;
@@ -186,13 +177,13 @@
 void ConstantQ::initialise( CQConfig Config )
 {
     m_FS = Config.FS;
-    m_FMin = Config.min;		// min freq
-    m_FMax = Config.max;		// max freq
-    m_BPO = Config.BPO;		// bins per octave
+    m_FMin = Config.min;                // min freq
+    m_FMax = Config.max;                // max freq
+    m_BPO = Config.BPO;         // bins per octave
     m_CQThresh = Config.CQThresh;// ConstantQ threshold for kernel generation
 
-    m_dQ = 1/(pow(2,(1/(double)m_BPO))-1);	// Work out Q value for Filter bank
-    m_uK = (unsigned int) ceil(m_BPO * log(m_FMax/m_FMin)/log(2.0));	// No. of constant Q bins
+    m_dQ = 1/(pow(2,(1/(double)m_BPO))-1);      // Work out Q value for Filter bank
+    m_uK = (unsigned int) ceil(m_BPO * log(m_FMax/m_FMin)/log(2.0));    // No. of constant Q bins
 
 //    std::cerr << "ConstantQ::initialise: rate = " << m_FS << ", fmin = " << m_FMin << ", fmax = " << m_FMax << ", bpo = " << m_BPO << ", K = " << m_uK << ", Q = " << m_dQ << std::endl;
 
@@ -223,10 +214,9 @@
 
     SparseKernel *sk = m_sparseKernel;
 
-    for (unsigned row=0; row<m_uK; row++) 
-    {
-	CQRe[ row ] = 0;
-	CQIm[ row ] = 0;
+    for (unsigned row=0; row<m_uK; row++) {
+        CQRe[ row ] = 0;
+        CQIm[ row ] = 0;
     }
 
     const unsigned *fftbin = &(sk->is[0]);
@@ -234,18 +224,17 @@
     const double   *real   = &(sk->real[0]);
     const double   *imag   = &(sk->imag[0]);
     const unsigned int sparseCells = sk->real.size();
-	
-    for (unsigned i = 0; i<sparseCells; i++)
-    {
-	const unsigned row = cqbin[i];
-	const unsigned col = fftbin[i];
+        
+    for (unsigned i = 0; i<sparseCells; i++) {
+        const unsigned row = cqbin[i];
+        const unsigned col = fftbin[i];
         if (col == 0) continue;
-	const double & r1  = real[i];
-	const double & i1  = imag[i];
-	const double & r2  = FFTRe[ m_FFTLength - col ];
-	const double & i2  = FFTIm[ m_FFTLength - col ];
-	// add the multiplication
-	CQRe[ row ] += (r1*r2 - i1*i2);
-	CQIm[ row ] += (r1*i2 + i1*r2);
+        const double & r1  = real[i];
+        const double & i1  = imag[i];
+        const double & r2  = FFTRe[ m_FFTLength - col ];
+        const double & i2  = FFTIm[ m_FFTLength - col ];
+        // add the multiplication
+        CQRe[ row ] += (r1*r2 - i1*i2);
+        CQIm[ row ] += (r1*i2 + i1*r2);
     }
 }
--- a/dsp/chromagram/ConstantQ.h	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/chromagram/ConstantQ.h	Fri May 31 11:54:32 2019 +0100
@@ -28,9 +28,8 @@
     double CQThresh;   // threshold
 };
 
-class ConstantQ {
-	
-//public functions incl. sparsekernel so can keep out of loop in main
+class ConstantQ
+{
 public:
     void process( const double* FFTRe, const double* FFTIm,
                   double* CQRe, double* CQIm );
@@ -43,10 +42,10 @@
     void sparsekernel();
 
     double hamming(int len, int n) {
-	double out = 0.54 - 0.46*cos(2*PI*n/len);
-	return(out);
+        double out = 0.54 - 0.46*cos(2*PI*n/len);
+        return(out);
     }
-	
+        
     int getnumwin() { return m_numWin;}
     double getQ() { return m_dQ;}
     int getK() {return m_uK ;}
@@ -56,7 +55,7 @@
 private:
     void initialise( CQConfig Config );
     void deInitialise();
-	
+        
     double* m_CQdata;
     double m_FS;
     double m_FMin;
--- a/dsp/onsets/DetectionFunction.cpp	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/onsets/DetectionFunction.cpp	Fri May 31 11:54:32 2019 +0100
@@ -54,7 +54,7 @@
 
     m_magHistory = new double[ m_halfLength ];
     memset(m_magHistory,0, m_halfLength*sizeof(double));
-		
+                
     m_phaseHistory = new double[ m_halfLength ];
     memset(m_phaseHistory,0, m_halfLength*sizeof(double));
 
@@ -134,30 +134,30 @@
     switch( m_DFType )
     {
     case DF_HFC:
-	retVal = HFC( m_halfLength, m_magnitude);
-	break;
-	
+        retVal = HFC( m_halfLength, m_magnitude);
+        break;
+        
     case DF_SPECDIFF:
-	retVal = specDiff( m_halfLength, m_magnitude);
-	break;
-	
+        retVal = specDiff( m_halfLength, m_magnitude);
+        break;
+        
     case DF_PHASEDEV:
         // Using the instantaneous phases here actually provides the
         // same results (for these calculations) as if we had used
         // unwrapped phases, but without the possible accumulation of
         // phase error over time
-	retVal = phaseDev( m_halfLength, m_thetaAngle);
-	break;
-	
+        retVal = phaseDev( m_halfLength, m_thetaAngle);
+        break;
+        
     case DF_COMPLEXSD:
-	retVal = complexSD( m_halfLength, m_magnitude, m_thetaAngle);
-	break;
+        retVal = complexSD( m_halfLength, m_magnitude, m_thetaAngle);
+        break;
 
     case DF_BROADBAND:
         retVal = broadband( m_halfLength, m_magnitude);
         break;
     }
-	
+        
     return retVal;
 }
 
@@ -166,9 +166,8 @@
     unsigned int i;
     double val = 0;
 
-    for( i = 0; i < length; i++)
-    {
-	val += src[ i ] * ( i + 1);
+    for( i = 0; i < length; i++) {
+        val += src[ i ] * ( i + 1);
     }
     return val;
 }
@@ -180,17 +179,17 @@
     double temp = 0.0;
     double diff = 0.0;
 
-    for( i = 0; i < length; i++)
-    {
-	temp = fabs( (src[ i ] * src[ i ]) - (m_magHistory[ i ] * m_magHistory[ i ]) );
-		
-	diff= sqrt(temp);
+    for( i = 0; i < length; i++) {
+        
+        temp = fabs( (src[ i ] * src[ i ]) - (m_magHistory[ i ] * m_magHistory[ i ]) );
+                
+        diff= sqrt(temp);
 
         // (See note in phaseDev below.)
 
         val += diff;
 
-	m_magHistory[ i ] = src[ i ];
+        m_magHistory[ i ] = src[ i ];
     }
 
     return val;
@@ -206,10 +205,9 @@
 
     double dev = 0;
 
-    for( i = 0; i < length; i++)
-    {
-	tmpPhase = (srcPhase[ i ]- 2*m_phaseHistory[ i ]+m_phaseHistoryOld[ i ]);
-	dev = MathUtilities::princarg( tmpPhase );
+    for( i = 0; i < length; i++) {
+        tmpPhase = (srcPhase[ i ]- 2*m_phaseHistory[ i ]+m_phaseHistoryOld[ i ]);
+        dev = MathUtilities::princarg( tmpPhase );
 
         // A previous version of this code only counted the value here
         // if the magnitude exceeded 0.1.  My impression is that
@@ -218,14 +216,14 @@
         // does significantly damage its ability to work with quieter
         // music, so I'm removing it and counting the result always.
         // Same goes for the spectral difference measure above.
-		
+                
         tmpVal  = fabs(dev);
         val += tmpVal ;
 
-	m_phaseHistoryOld[ i ] = m_phaseHistory[ i ] ;
-	m_phaseHistory[ i ] = srcPhase[ i ];
+        m_phaseHistoryOld[ i ] = m_phaseHistory[ i ] ;
+        m_phaseHistory[ i ] = srcPhase[ i ];
     }
-	
+        
     return val;
 }
 
@@ -242,21 +240,21 @@
     ComplexData meas = ComplexData( 0, 0 );
     ComplexData j = ComplexData( 0, 1 );
 
-    for( i = 0; i < length; i++)
-    {
-	tmpPhase = (srcPhase[ i ]- 2*m_phaseHistory[ i ]+m_phaseHistoryOld[ i ]);
-	dev= MathUtilities::princarg( tmpPhase );
-		
-	meas = m_magHistory[i] - ( srcMagnitude[ i ] * exp( j * dev) );
+    for( i = 0; i < length; i++) {
+        
+        tmpPhase = (srcPhase[ i ]- 2*m_phaseHistory[ i ]+m_phaseHistoryOld[ i ]);
+        dev= MathUtilities::princarg( tmpPhase );
+                
+        meas = m_magHistory[i] - ( srcMagnitude[ i ] * exp( j * dev) );
 
-	tmpReal = real( meas );
-	tmpImag = imag( meas );
+        tmpReal = real( meas );
+        tmpImag = imag( meas );
 
-	val += sqrt( (tmpReal * tmpReal) + (tmpImag * tmpImag) );
-		
-	m_phaseHistoryOld[ i ] = m_phaseHistory[ i ] ;
-	m_phaseHistory[ i ] = srcPhase[ i ];
-	m_magHistory[ i ] = srcMagnitude[ i ];
+        val += sqrt( (tmpReal * tmpReal) + (tmpImag * tmpImag) );
+                
+        m_phaseHistoryOld[ i ] = m_phaseHistory[ i ] ;
+        m_phaseHistory[ i ] = srcPhase[ i ];
+        m_magHistory[ i ] = srcMagnitude[ i ];
     }
 
     return val;
--- a/dsp/onsets/DetectionFunction.h	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/onsets/DetectionFunction.h	Fri May 31 11:54:32 2019 +0100
@@ -65,7 +65,7 @@
     double phaseDev(unsigned int length, double *srcPhase);
     double complexSD(unsigned int length, double *srcMagnitude, double *srcPhase);
     double broadband(unsigned int length, double *srcMagnitude);
-	
+        
 private:
     void initialise( DFConfig Config );
     void deInitialise();
@@ -90,7 +90,7 @@
     double* m_unwrapped; // Unwrapped phase of analysis frame
 
     Window<double> *m_window;
-    PhaseVocoder* m_phaseVoc;	// Phase Vocoder
+    PhaseVocoder* m_phaseVoc;   // Phase Vocoder
 };
 
 #endif 
--- a/dsp/onsets/PeakPicking.cpp	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/onsets/PeakPicking.cpp	Fri May 31 11:54:32 2019 +0100
@@ -49,7 +49,7 @@
     Qfilta = Config.QuadThresh.a ;
     Qfiltb = Config.QuadThresh.b ;
     Qfiltc = Config.QuadThresh.c ;
-	
+        
     m_DFProcessingParams.length = m_DFLength; 
     m_DFProcessingParams.LPOrd = Config.LPOrd; 
     m_DFProcessingParams.LPACoeffs = Config.LPACoeffs; 
@@ -77,21 +77,19 @@
 {
     if (len < 4) return;
 
-    vector <double> m_maxima;	
+    vector <double> m_maxima;   
 
     // Signal conditioning 
     m_DFSmoothing->process( src, m_workBuffer );
-	
-    for( unsigned int u = 0; u < len; u++)
-    {
-	m_maxima.push_back( m_workBuffer[ u ] );		
+        
+    for( unsigned int u = 0; u < len; u++) {
+        m_maxima.push_back( m_workBuffer[ u ] );                
     }
-	
+        
     quadEval( m_maxima, onsets );
 
-    for( int b = 0; b <  (int)m_maxima.size(); b++)
-    {
-	src[ b ] = m_maxima[ b ];
+    for( int b = 0; b <  (int)m_maxima.size(); b++) {
+        src[ b ] = m_maxima[ b ];
     }
 }
 
@@ -101,7 +99,7 @@
 
     vector <int> m_maxIndex;
     vector <int> m_onsetPosition;
-	
+        
     vector <double> m_maxFit;
     vector <double> m_poly;
     vector <double> m_err;
@@ -110,42 +108,36 @@
     m_poly.push_back(0);
     m_poly.push_back(0);
 
-    for(  int t = -2; t < 3; t++)
-    {
-	m_err.push_back( (double)t );
+    for(  int t = -2; t < 3; t++) {
+        m_err.push_back( (double)t );
     }
-    for( unsigned int i = 2; i < src.size() - 2; i++)
-    {
-	if( (src[i] > src[i-1]) && (src[i] > src[i+1]) && (src[i] > 0) )
-	{
-//	    m_maxIndex.push_back(  i + 1 );
+
+    for( unsigned int i = 2; i < src.size() - 2; i++) {
+        if( (src[i] > src[i-1]) && (src[i] > src[i+1]) && (src[i] > 0) ) {
             m_maxIndex.push_back(i);
-	}
+        }
     }
 
     maxLength = m_maxIndex.size();
 
     double selMax = 0;
 
-    for( unsigned int j = 0; j < maxLength ; j++)
-    {
-        for (int k = -2; k <= 2; ++k)
-	{
-	    selMax = src[ m_maxIndex[j] + k ] ;
-	    m_maxFit.push_back(selMax);			
-	}
+    for( unsigned int j = 0; j < maxLength ; j++) {
+        for (int k = -2; k <= 2; ++k) {
+            selMax = src[ m_maxIndex[j] + k ] ;
+            m_maxFit.push_back(selMax);                 
+        }
 
-	TPolyFit::PolyFit2(m_err, m_maxFit, m_poly);
+        TPolyFit::PolyFit2(m_err, m_maxFit, m_poly);
 
-	double f = m_poly[0];
-	double h = m_poly[2];
+        double f = m_poly[0];
+        double h = m_poly[2];
 
-	if (h < -Qfilta || f > Qfiltc)
-	{
-	    idx.push_back(m_maxIndex[j]);
-	}
-		
-	m_maxFit.clear();
+        if (h < -Qfilta || f > Qfiltc) {
+            idx.push_back(m_maxIndex[j]);
+        }
+                
+        m_maxFit.clear();
     }
 
     return 1;
--- a/dsp/onsets/PeakPicking.h	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/onsets/PeakPicking.h	Fri May 31 11:54:32 2019 +0100
@@ -92,15 +92,14 @@
 public:
     PeakPicking( PPickParams Config );
     virtual ~PeakPicking();
-	
+        
     void process( double* src, unsigned int len, vector<int> &onsets  );
 
-
 private:
     void initialise( PPickParams Config  );
     void deInitialise();
     int  quadEval( vector<double> &src, vector<int> &idx );
-	
+        
     DFProcConfig m_DFProcessingParams;
 
     unsigned int m_DFLength ;
@@ -108,10 +107,9 @@
     double Qfiltb;
     double Qfiltc;
 
-
     double* m_workBuffer;
-	
-    DFProcess*	m_DFSmoothing;
+        
+    DFProcess*  m_DFSmoothing;
 };
 
 #endif
--- a/dsp/phasevocoder/PhaseVocoder.cpp	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/phasevocoder/PhaseVocoder.cpp	Fri May 31 11:54:32 2019 +0100
@@ -105,17 +105,17 @@
 }
 
 void PhaseVocoder::getMagnitudes(double *mag)
-{	
+{       
     for (int i = 0; i < m_n/2 + 1; i++) {
-	mag[i] = sqrt(m_real[i] * m_real[i] + m_imag[i] * m_imag[i]);
+        mag[i] = sqrt(m_real[i] * m_real[i] + m_imag[i] * m_imag[i]);
     }
 }
 
 void PhaseVocoder::getPhases(double *theta)
 {
     for (int i = 0; i < m_n/2 + 1; i++) {
-	theta[i] = atan2(m_imag[i], m_real[i]);
-    }	
+        theta[i] = atan2(m_imag[i], m_real[i]);
+    }   
 }
 
 void PhaseVocoder::unwrapPhases(double *theta, double *unwrapped)
--- a/dsp/rateconversion/Decimator.cpp	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/rateconversion/Decimator.cpp	Fri May 31 11:54:32 2019 +0100
@@ -23,7 +23,6 @@
 
 Decimator::Decimator( unsigned int inLength, unsigned int decFactor )
 {
-
     m_inputLength = 0;
     m_outputLength = 0;
     m_decFactor = 1;
@@ -47,8 +46,8 @@
     // If adding new factors here, add them to
     // getHighestSupportedFactor in the header as well
 
-    if(m_decFactor == 8)
-    {
+    if(m_decFactor == 8) {
+
         //////////////////////////////////////////////////
         b[0] = 0.060111378492136;
         b[1] = -0.257323420830598;
@@ -68,74 +67,74 @@
         a[6] = 2.577553446979888;
         a[7] = -0.326903916815751;
         //////////////////////////////////////////////////
-    }
-    else if( m_decFactor == 4 )
-    {
-	//////////////////////////////////////////////////
-	b[ 0 ] = 0.10133306904918619;
-	b[ 1 ] = -0.2447523353702363;
-	b[ 2 ] = 0.33622528590120965;
-	b[ 3 ] = -0.13936581560633518;
-	b[ 4 ] = -0.13936581560633382;
-	b[ 5 ] = 0.3362252859012087;
-	b[ 6 ] = -0.2447523353702358;
-	b[ 7 ] = 0.10133306904918594;
 
-	a[ 0 ] = 1;
-	a[ 1 ] = -3.9035590278139427;
-	a[ 2 ] = 7.5299379980621133;
-	a[ 3 ] = -8.6890803793177511;
-	a[ 4 ] = 6.4578667096099176;
-	a[ 5 ] = -3.0242979431223631;
-	a[ 6 ] = 0.83043385136748382;
-	a[ 7 ] = -0.094420800837809335;
-	//////////////////////////////////////////////////
-    }
-    else if( m_decFactor == 2 )
-    {
-	//////////////////////////////////////////////////
-	b[ 0 ] = 0.20898944260075727;
-	b[ 1 ] = 0.40011234879814367;
-	b[ 2 ] = 0.819741973072733;
-	b[ 3 ] = 1.0087419911682323;
-	b[ 4 ] = 1.0087419911682325;
-	b[ 5 ] = 0.81974197307273156;
-	b[ 6 ] = 0.40011234879814295;
-	b[ 7 ] = 0.20898944260075661;
+    } else if( m_decFactor == 4 ) {
+        
+        //////////////////////////////////////////////////
+        b[ 0 ] = 0.10133306904918619;
+        b[ 1 ] = -0.2447523353702363;
+        b[ 2 ] = 0.33622528590120965;
+        b[ 3 ] = -0.13936581560633518;
+        b[ 4 ] = -0.13936581560633382;
+        b[ 5 ] = 0.3362252859012087;
+        b[ 6 ] = -0.2447523353702358;
+        b[ 7 ] = 0.10133306904918594;
 
-	a[ 0 ] = 1;
-	a[ 1 ] = 0.0077331184208358217;
-	a[ 2 ] = 1.9853971155964376;
-	a[ 3 ] = 0.19296739275341004;
-	a[ 4 ] = 1.2330748872852182;
-	a[ 5 ] = 0.18705341389316466;
-	a[ 6 ] = 0.23659265908013868;
-	a[ 7 ] = 0.032352924250533946;
-    }
-    else
-    {
+        a[ 0 ] = 1;
+        a[ 1 ] = -3.9035590278139427;
+        a[ 2 ] = 7.5299379980621133;
+        a[ 3 ] = -8.6890803793177511;
+        a[ 4 ] = 6.4578667096099176;
+        a[ 5 ] = -3.0242979431223631;
+        a[ 6 ] = 0.83043385136748382;
+        a[ 7 ] = -0.094420800837809335;
+        //////////////////////////////////////////////////
+
+    } else if( m_decFactor == 2 ) {
+        
+        //////////////////////////////////////////////////
+        b[ 0 ] = 0.20898944260075727;
+        b[ 1 ] = 0.40011234879814367;
+        b[ 2 ] = 0.819741973072733;
+        b[ 3 ] = 1.0087419911682323;
+        b[ 4 ] = 1.0087419911682325;
+        b[ 5 ] = 0.81974197307273156;
+        b[ 6 ] = 0.40011234879814295;
+        b[ 7 ] = 0.20898944260075661;
+
+        a[ 0 ] = 1;
+        a[ 1 ] = 0.0077331184208358217;
+        a[ 2 ] = 1.9853971155964376;
+        a[ 3 ] = 0.19296739275341004;
+        a[ 4 ] = 1.2330748872852182;
+        a[ 5 ] = 0.18705341389316466;
+        a[ 6 ] = 0.23659265908013868;
+        a[ 7 ] = 0.032352924250533946;
+
+    } else {
+        
         if ( m_decFactor != 1 ) {
             std::cerr << "WARNING: Decimator::initialise: unsupported decimation factor " << m_decFactor << ", no antialiasing filter will be used" << std::endl;
         }
 
-	//////////////////////////////////////////////////
-	b[ 0 ] = 1;
-	b[ 1 ] = 0;
-	b[ 2 ] = 0;
-	b[ 3 ] = 0;
-	b[ 4 ] = 0;
-	b[ 5 ] = 0;
-	b[ 6 ] = 0;
-	b[ 7 ] = 0;
+        //////////////////////////////////////////////////
+        b[ 0 ] = 1;
+        b[ 1 ] = 0;
+        b[ 2 ] = 0;
+        b[ 3 ] = 0;
+        b[ 4 ] = 0;
+        b[ 5 ] = 0;
+        b[ 6 ] = 0;
+        b[ 7 ] = 0;
 
-	a[ 0 ] = 1;
-	a[ 1 ] = 0;
-	a[ 2 ] = 0;
-	a[ 3 ] = 0;
-	a[ 4 ] = 0;
-	a[ 5 ] = 0;
-	a[ 6 ] = 0;
-	a[ 7 ] = 0;
+        a[ 0 ] = 1;
+        a[ 1 ] = 0;
+        a[ 2 ] = 0;
+        a[ 3 ] = 0;
+        a[ 4 ] = 0;
+        a[ 5 ] = 0;
+        a[ 6 ] = 0;
+        a[ 7 ] = 0;
     }
 
     resetFilter();
@@ -155,46 +154,42 @@
 
 void Decimator::doAntiAlias(const double *src, double *dst, unsigned int length)
 {
+    for( unsigned int i = 0; i < length; i++ ) {
+        
+        Input = (double)src[ i ];
 
-    for( unsigned int i = 0; i < length; i++ )
-    {
-	Input = (double)src[ i ];
+        Output = Input * b[ 0 ] + o1;
 
-	Output = Input * b[ 0 ] + o1;
+        o1 = Input * b[ 1 ] - Output * a[ 1 ] + o2;
+        o2 = Input * b[ 2 ] - Output * a[ 2 ] + o3;
+        o3 = Input * b[ 3 ] - Output * a[ 3 ] + o4;
+        o4 = Input * b[ 4 ] - Output * a[ 4 ] + o5;
+        o5 = Input * b[ 5 ] - Output * a[ 5 ] + o6;
+        o6 = Input * b[ 6 ] - Output * a[ 6 ] + o7;
+        o7 = Input * b[ 7 ] - Output * a[ 7 ] ;
 
-	o1 = Input * b[ 1 ] - Output * a[ 1 ] + o2;
-	o2 = Input * b[ 2 ] - Output * a[ 2 ] + o3;
-	o3 = Input * b[ 3 ] - Output * a[ 3 ] + o4;
-	o4 = Input * b[ 4 ] - Output * a[ 4 ] + o5;
-	o5 = Input * b[ 5 ] - Output * a[ 5 ] + o6;
-	o6 = Input * b[ 6 ] - Output * a[ 6 ] + o7;
-	o7 = Input * b[ 7 ] - Output * a[ 7 ] ;
-
-	dst[ i ] = Output;
+        dst[ i ] = Output;
     }
-
 }
 
 void Decimator::doAntiAlias(const float *src, double *dst, unsigned int length)
 {
+    for( unsigned int i = 0; i < length; i++ ) {
+        
+        Input = (double)src[ i ];
 
-    for( unsigned int i = 0; i < length; i++ )
-    {
-	Input = (double)src[ i ];
+        Output = Input * b[ 0 ] + o1;
 
-	Output = Input * b[ 0 ] + o1;
+        o1 = Input * b[ 1 ] - Output * a[ 1 ] + o2;
+        o2 = Input * b[ 2 ] - Output * a[ 2 ] + o3;
+        o3 = Input * b[ 3 ] - Output * a[ 3 ] + o4;
+        o4 = Input * b[ 4 ] - Output * a[ 4 ] + o5;
+        o5 = Input * b[ 5 ] - Output * a[ 5 ] + o6;
+        o6 = Input * b[ 6 ] - Output * a[ 6 ] + o7;
+        o7 = Input * b[ 7 ] - Output * a[ 7 ] ;
 
-	o1 = Input * b[ 1 ] - Output * a[ 1 ] + o2;
-	o2 = Input * b[ 2 ] - Output * a[ 2 ] + o3;
-	o3 = Input * b[ 3 ] - Output * a[ 3 ] + o4;
-	o4 = Input * b[ 4 ] - Output * a[ 4 ] + o5;
-	o5 = Input * b[ 5 ] - Output * a[ 5 ] + o6;
-	o6 = Input * b[ 6 ] - Output * a[ 6 ] + o7;
-	o7 = Input * b[ 7 ] - Output * a[ 7 ] ;
-
-	dst[ i ] = Output;
+        dst[ i ] = Output;
     }
-
 }
 
 void Decimator::process(const double *src, double *dst)
@@ -210,9 +205,8 @@
 
     unsigned idx = 0;
 
-    for( unsigned int i = 0; i < m_outputLength; i++ )
-    {
-	dst[ idx++ ] = decBuffer[ m_decFactor * i ];
+    for( unsigned int i = 0; i < m_outputLength; i++ ) {
+        dst[ idx++ ] = decBuffer[ m_decFactor * i ];
     }
 }
 
@@ -229,8 +223,7 @@
 
     unsigned idx = 0;
 
-    for( unsigned int i = 0; i < m_outputLength; i++ )
-    {
-	dst[ idx++ ] = decBuffer[ m_decFactor * i ];
+    for( unsigned int i = 0; i < m_outputLength; i++ ) {
+        dst[ idx++ ] = decBuffer[ m_decFactor * i ];
     }
 }
--- a/dsp/rateconversion/Decimator.h	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/rateconversion/Decimator.h	Fri May 31 11:54:32 2019 +0100
@@ -75,7 +75,7 @@
 
     double a[ 9 ];
     double b[ 9 ];
-	
+        
     double* decBuffer;
 };
 
--- a/dsp/rateconversion/DecimatorB.cpp	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/rateconversion/DecimatorB.cpp	Fri May 31 11:54:32 2019 +0100
@@ -94,17 +94,17 @@
 
     for (int i = 0; i < length; i++) {
 
-	double input = src[i];
-	double output = input * m_b[0] + o[0];
+        double input = src[i];
+        double output = input * m_b[0] + o[0];
 
-	o[0] = input * m_b[1] - output * m_a[1] + o[1];
-	o[1] = input * m_b[2] - output * m_a[2] + o[2];
-	o[2] = input * m_b[3] - output * m_a[3] + o[3];
-	o[3] = input * m_b[4] - output * m_a[4] + o[4];
-	o[4] = input * m_b[5] - output * m_a[5] + o[5];
-	o[5] = input * m_b[6] - output * m_a[6];
+        o[0] = input * m_b[1] - output * m_a[1] + o[1];
+        o[1] = input * m_b[2] - output * m_a[2] + o[2];
+        o[2] = input * m_b[3] - output * m_a[3] + o[3];
+        o[3] = input * m_b[4] - output * m_a[4] + o[4];
+        o[4] = input * m_b[5] - output * m_a[5] + o[5];
+        o[5] = input * m_b[6] - output * m_a[6];
 
-	dst[i] = output;
+        dst[i] = output;
     }
 }
 
--- a/dsp/rateconversion/DecimatorB.h	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/rateconversion/DecimatorB.h	Fri May 31 11:54:32 2019 +0100
@@ -55,7 +55,7 @@
 
     double m_a[7];
     double m_b[7];
-	
+        
     double *m_aaBuffer;
     double *m_tmpBuffer;
 };
--- a/dsp/rateconversion/Resampler.cpp	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/rateconversion/Resampler.cpp	Fri May 31 11:54:32 2019 +0100
@@ -70,10 +70,10 @@
     }
 
     KaiserWindow::Parameters params =
-	KaiserWindow::parametersForBandwidth(snr, bandwidth, higher / m_gcd);
+        KaiserWindow::parametersForBandwidth(snr, bandwidth, higher / m_gcd);
 
     params.length =
-	(params.length % 2 == 0 ? params.length + 1 : params.length);
+        (params.length % 2 == 0 ? params.length + 1 : params.length);
     
     params.length =
         (params.length > 200001 ? 200001 : params.length);
@@ -174,23 +174,23 @@
 
     for (int phase = 0; phase < inputSpacing; ++phase) {
 
-	Phase p;
+        Phase p;
 
-	p.nextPhase = phase - outputSpacing;
-	while (p.nextPhase < 0) p.nextPhase += inputSpacing;
-	p.nextPhase %= inputSpacing;
-	
-	p.drop = int(ceil(std::max(0.0, double(outputSpacing - phase))
-			  / inputSpacing));
+        p.nextPhase = phase - outputSpacing;
+        while (p.nextPhase < 0) p.nextPhase += inputSpacing;
+        p.nextPhase %= inputSpacing;
+        
+        p.drop = int(ceil(std::max(0.0, double(outputSpacing - phase))
+                          / inputSpacing));
 
-	int filtZipLength = int(ceil(double(m_filterLength - phase)
-				     / inputSpacing));
+        int filtZipLength = int(ceil(double(m_filterLength - phase)
+                                     / inputSpacing));
 
-	for (int i = 0; i < filtZipLength; ++i) {
-	    p.filter.push_back(filter[i * inputSpacing + phase]);
-	}
+        for (int i = 0; i < filtZipLength; ++i) {
+            p.filter.push_back(filter[i * inputSpacing + phase]);
+        }
 
-	m_phaseData[phase] = p;
+        m_phaseData[phase] = p;
     }
 
 #ifdef DEBUG_RESAMPLER
@@ -268,7 +268,7 @@
 
 #ifdef DEBUG_RESAMPLER
     cerr << "initial phase " << m_phase << " (as " << (m_filterLength/2) << " % " << inputSpacing << ")"
-	      << ", latency " << m_latency << endl;
+              << ", latency " << m_latency << endl;
 #endif
 }
 
@@ -297,8 +297,8 @@
     const double *const R__ filt(pd.filter.data());
 
     for (int i = 0; i < n; ++i) {
-	// NB gcc can only vectorize this with -ffast-math
-	v += buf[i] * filt[i];
+        // NB gcc can only vectorize this with -ffast-math
+        v += buf[i] * filt[i];
     }
 
     m_bufferOrigin += pd.drop;
@@ -321,9 +321,9 @@
     double scaleFactor = (double(m_targetRate) / m_gcd) / m_peakToPole;
 
     while (outidx < maxout &&
-	   m_buffer.size() >= m_phaseData[m_phase].filter.size() + m_bufferOrigin) {
-	dst[outidx] = scaleFactor * reconstructOne();
-	outidx++;
+           m_buffer.size() >= m_phaseData[m_phase].filter.size() + m_bufferOrigin) {
+        dst[outidx] = scaleFactor * reconstructOne();
+        outidx++;
     }
 
     if (m_bufferOrigin > (int)m_buffer.size()) {
@@ -392,7 +392,7 @@
     int printN = 50;
     cerr << "first " << printN << " in:" << endl;
     for (int i = 0; i < printN && i < n; ++i) {
-	if (i % 5 == 0) cerr << endl << i << "... ";
+        if (i % 5 == 0) cerr << endl << i << "... ";
         cerr << data[i] << " ";
     }
     cerr << endl;
@@ -402,13 +402,13 @@
     if (toReturn > m) toReturn = m;
 
     vector<double> sliced(out.begin() + latency, 
-			  out.begin() + latency + toReturn);
+                          out.begin() + latency + toReturn);
 
 #ifdef DEBUG_RESAMPLER_VERBOSE
     cerr << "first " << printN << " out (after latency compensation), length " << sliced.size() << ":";
     for (int i = 0; i < printN && i < sliced.size(); ++i) {
-	if (i % 5 == 0) cerr << endl << i << "... ";
-	cerr << sliced[i] << " ";
+        if (i % 5 == 0) cerr << endl << i << "... ";
+        cerr << sliced[i] << " ";
     }
     cerr << endl;
 #endif
--- a/dsp/signalconditioning/DFProcess.cpp	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/signalconditioning/DFProcess.cpp	Fri May 31 11:54:32 2019 +0100
@@ -35,7 +35,7 @@
 DFProcess::DFProcess( DFProcConfig Config )
 {
     filtSrc = NULL;
-    filtDst = NULL;	
+    filtDst = NULL;     
     m_filtScratchIn = NULL;
     m_filtScratchOut = NULL;
 
@@ -62,11 +62,13 @@
     filtDst = new double[ m_length ];
 
     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);
+    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;
 }
@@ -74,13 +76,9 @@
 void DFProcess::deInitialise()
 {
     delete [] filtSrc;
-
     delete [] filtDst;
-
     delete [] m_filtScratchIn;
-
     delete [] m_filtScratchOut;
-
     delete m_FiltFilt;
 }
 
@@ -108,76 +106,71 @@
 
     double* scratch = new double[ m_length ];
 
-    for( i = 0; i < m_winPre; i++)
-    {
-        if (index >= m_length) break;
+    for( i = 0; i < m_winPre; i++) {
+        
+        if (index >= m_length) {
+            break;
+        }
 
-	k = i + m_winPost + 1;
+        k = i + m_winPost + 1;
 
-	for( j = 0; j < k; j++)
-	{
-	    y[ j ] = src[ j ];
-	}
-	scratch[ index ] = MathUtilities::median( y, k );
-	index++;
+        for( j = 0; j < k; j++) {
+            y[ j ] = src[ j ];
+        }
+        scratch[ index ] = MathUtilities::median( y, k );
+        index++;
     }
 
-    for(  i = 0; i + m_winPost + m_winPre < m_length; i ++)
-    {
-        if (index >= m_length) break;
+    for(  i = 0; i + m_winPost + m_winPre < m_length; i ++) {
+        
+        if (index >= m_length) {
+            break;
+        }
+                         
+        l = 0;
+        for(  j  = i; j < ( i + m_winPost + m_winPre + 1); j++) {
+            y[ l ] = src[ j ];
+            l++;
+        }
 
-			 
-	l = 0;
-	for(  j  = i; j < ( i + m_winPost + m_winPre + 1); j++)
-	{
-	    y[ l ] = src[ j ];
-	    l++;
-	}
-
-	scratch[ index++ ] = MathUtilities::median( y, (m_winPost + m_winPre + 1 ));
+        scratch[index] = MathUtilities::median( y, (m_winPost + m_winPre + 1 ));
+        index++;
     }
 
-    for( i = std::max( m_length - m_winPost, 1); i < m_length; i++)
-    {
-        if (index >= m_length) break;
+    for( i = std::max( m_length - m_winPost, 1); i < m_length; i++) {
+        
+        if (index >= m_length) {
+            break;
+        }
 
-	k = std::max( i - m_winPre, 1);
+        k = std::max( i - m_winPre, 1);
 
-	l = 0;
-	for( j = k; j < m_length; j++)
-	{
-	    y[ l ] = src[ j ];
-
-	    l++;
-	}
-		
-	scratch[ index++ ] = MathUtilities::median( y, l); 
+        l = 0;
+        for( j = k; j < m_length; j++) {
+            y[ l ] = src[ j ];
+            l++;
+        }
+                
+        scratch[index] = MathUtilities::median( y, l);
+        index++;
     }
 
-
-    for( i = 0; i < m_length; i++ )
-    {
-	//add a delta threshold used as an offset when computing the smoothed detection function
-	//(helps to discard noise when detecting peaks)	
-	val = src[ i ] - scratch[ i ] - m_delta;
-		
-	if( m_isMedianPositive )
-	{
-	    if( val > 0 )
-	    {
-		dst[ i ]  = val;
-	    }
-	    else
-	    {
-		dst[ i ]  = 0;
-	    }
-	}
-	else
-	{
-	    dst[ i ]  = val;
-	}
+    for( i = 0; i < m_length; i++ ) {
+        //add a delta threshold used as an offset when computing the smoothed detection function
+        //(helps to discard noise when detecting peaks) 
+        val = src[ i ] - scratch[ i ] - m_delta;
+                
+        if( m_isMedianPositive ) {
+            if( val > 0 ) {
+                dst[ i ]  = val;
+            } else {
+                dst[ i ]  = 0;
+            }
+        } else {
+            dst[ i ]  = val;
+        }
     }
-	
+        
     delete [] y;
     delete [] scratch;
 }
@@ -193,8 +186,7 @@
 
     MathUtilities::getAlphaNorm( src, m_length, m_alphaNormParam, &DFAlphaNorm );
 
-    for (int i = 0; i < m_length; i++)
-    {
-	dst[ i ] = ( src[ i ] - DFMin ) / DFAlphaNorm; 
+    for (int i = 0; i < m_length; i++) {
+        dst[ i ] = ( src[ i ] - DFMin ) / DFAlphaNorm; 
     }
 }
--- a/dsp/signalconditioning/DFProcess.h	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/signalconditioning/DFProcess.h	Fri May 31 11:54:32 2019 +0100
@@ -60,7 +60,7 @@
 
     void process( double* src, double* dst );
 
-	
+        
 private:
     void initialise( DFProcConfig Config );
     void deInitialise();
--- a/dsp/signalconditioning/FiltFilt.cpp	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/signalconditioning/FiltFilt.cpp	Fri May 31 11:54:32 2019 +0100
@@ -30,22 +30,21 @@
 }
 
 void FiltFilt::process(double *src, double *dst, unsigned int length)
-{	
+{       
     unsigned int i;
 
     if (length == 0) return;
 
     unsigned int nFilt = m_ord + 1;
     unsigned int nFact = 3 * ( nFilt - 1);
-    unsigned int nExt	= length + 2 * nFact;
+    unsigned int nExt   = length + 2 * nFact;
 
     double *filtScratchIn = new double[ nExt ];
     double *filtScratchOut = new double[ nExt ];
-	
-    for( i = 0; i< nExt; i++ ) 
-    {
-	filtScratchIn[ i ] = 0.0;
-	filtScratchOut[ i ] = 0.0;
+        
+    for( i = 0; i< nExt; i++ ) {
+        filtScratchIn[ i ] = 0.0;
+        filtScratchOut[ i ] = 0.0;
     }
 
     // Edge transients reflection
@@ -53,54 +52,46 @@
     double sampleN = 2 * src[ length - 1 ];
 
     unsigned int index = 0;
-    for( i = nFact; i > 0; i-- )
-    {
-	filtScratchIn[ index++ ] = sample0 - src[ i ];
+    for( i = nFact; i > 0; i-- ) {
+        filtScratchIn[ index++ ] = sample0 - src[ i ];
     }
     index = 0;
-    for( i = 0; i < nFact; i++ )
-    {
-	filtScratchIn[ (nExt - nFact) + index++ ] = sampleN - src[ (length - 2) - i ];
+    for( i = 0; i < nFact; i++ ) {
+        filtScratchIn[ (nExt - nFact) + index++ ] = sampleN - src[ (length - 2) - i ];
     }
 
     index = 0;
-    for( i = 0; i < length; i++ )
-    {
-	filtScratchIn[ i + nFact ] = src[ i ];
+    for( i = 0; i < length; i++ ) {
+        filtScratchIn[ i + nFact ] = src[ i ];
     }
-	
+        
     ////////////////////////////////
     // Do  0Ph filtering
     m_filter.process( filtScratchIn, filtScratchOut, nExt);
-	
+        
     // reverse the series for FILTFILT 
-    for ( i = 0; i < nExt; i++)
-    { 
-	filtScratchIn[ i ] = filtScratchOut[ nExt - i - 1];
+    for ( i = 0; i < nExt; i++) { 
+        filtScratchIn[ i ] = filtScratchOut[ nExt - i - 1];
     }
 
     // do FILTER again 
     m_filter.process( filtScratchIn, filtScratchOut, nExt);
-	
+        
     // reverse the series back 
-    for ( i = 0; i < nExt; i++)
-    {
-	filtScratchIn[ i ] = filtScratchOut[ nExt - i - 1 ];
+    for ( i = 0; i < nExt; i++) {
+        filtScratchIn[ i ] = filtScratchOut[ nExt - i - 1 ];
     }
-    for ( i = 0;i < nExt; i++)
-    {
-	filtScratchOut[ i ] = filtScratchIn[ i ];
+    for ( i = 0;i < nExt; i++) {
+        filtScratchOut[ i ] = filtScratchIn[ i ];
     }
 
     index = 0;
-    for( i = 0; i < length; i++ )
-    {
-	dst[ index++ ] = filtScratchOut[ i + nFact ];
-    }	
+    for( i = 0; i < length; i++ ) {
+        dst[ index++ ] = filtScratchOut[ i + nFact ];
+    }   
 
     delete [] filtScratchIn;
     delete [] filtScratchOut;
-
 }
 
 void FiltFilt::reset()
--- a/dsp/signalconditioning/Filter.cpp	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/signalconditioning/Filter.cpp	Fri May 31 11:54:32 2019 +0100
@@ -75,13 +75,14 @@
 
 void
 Filter::process(const double *const QM_R__ in,
-		double *const QM_R__ out,
-		const int n)
+                double *const QM_R__ out,
+                const int n)
 {
     for (int s = 0; s < n; ++s) {
 
-        if (m_offb > 0) --m_offb;
-        else {
+        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];
             }
@@ -109,8 +110,9 @@
 
             outval = b_sum - a_sum;
 
-            if (m_offa > 0) --m_offa;
-            else {
+            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];
                 }
@@ -119,7 +121,7 @@
             m_bufa[m_offa] = outval;
         }
         
-	out[s] = outval;
+        out[s] = outval;
     }
 }
 
--- a/dsp/signalconditioning/Framer.cpp	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/signalconditioning/Framer.cpp	Fri May 31 11:54:32 2019 +0100
@@ -28,11 +28,13 @@
 
 Framer::~Framer()
 {
-    if( m_dataFrame != NULL )
-	delete [] m_dataFrame;
+    if( m_dataFrame != NULL ) {
+        delete [] m_dataFrame;
+    }
 
-    if( m_strideFrame != NULL )
-	delete [] m_strideFrame;
+    if( m_strideFrame != NULL ) {
+        delete [] m_strideFrame;
+    }
 }
 
 void Framer::configure( unsigned int frameLength, unsigned int hop )
@@ -42,17 +44,15 @@
 
     resetCounters();
 
-    if( m_dataFrame != NULL )
-    {
-	delete [] m_dataFrame;	
-	m_dataFrame = NULL;
+    if( m_dataFrame != NULL ) {
+        delete [] m_dataFrame;  
+        m_dataFrame = NULL;
     }
     m_dataFrame = new double[ m_frameLength ];
 
-    if( m_strideFrame != NULL )
-    {
-	delete [] m_strideFrame;	
-	m_strideFrame = NULL;
+    if( m_strideFrame != NULL ) {
+        delete [] m_strideFrame;        
+        m_strideFrame = NULL;
     }
     m_strideFrame = new double[ m_stepSize ];
 }
@@ -60,30 +60,27 @@
 void Framer::getFrame(double *dst)
 {
 
-    if( (m_ulSrcIndex + ( m_frameLength) ) < m_ulSampleLen )
-    {
-	for( unsigned int u = 0; u < m_frameLength; u++)
-	{
-	    dst[ u ] = m_srcBuffer[ m_ulSrcIndex++ ]; 
-	}	
-	m_ulSrcIndex -= ( m_frameLength - m_stepSize );
-    }
-    else
-    {
-	unsigned int rem = (m_ulSampleLen - m_ulSrcIndex );
-	unsigned int zero = m_frameLength - rem;
+    if( (m_ulSrcIndex + ( m_frameLength) ) < m_ulSampleLen ) {
 
-	for( unsigned int u = 0; u < rem; u++ )
-	{
-	    dst[ u ] = m_srcBuffer[ m_ulSrcIndex++ ];
-	}
-		
-	for( unsigned int u = 0; u < zero; u++ )
-	{
-	    dst[ rem + u ] = 0;
-	}
+        for( unsigned int u = 0; u < m_frameLength; u++) {
+            dst[ u ] = m_srcBuffer[ m_ulSrcIndex++ ]; 
+        }       
+        m_ulSrcIndex -= ( m_frameLength - m_stepSize );
 
-	m_ulSrcIndex -= (( rem - m_stepSize ) );
+    } else {
+
+        unsigned int rem = (m_ulSampleLen - m_ulSrcIndex );
+        unsigned int zero = m_frameLength - rem;
+
+        for( unsigned int u = 0; u < rem; u++ ) {
+            dst[ u ] = m_srcBuffer[ m_ulSrcIndex++ ];
+        }
+                
+        for( unsigned int u = 0; u < zero; u++ ) {
+            dst[ rem + u ] = 0;
+        }
+
+        m_ulSrcIndex -= (( rem - m_stepSize ) );
     }
 
     m_framesRead++;
--- a/dsp/signalconditioning/Framer.h	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/signalconditioning/Framer.h	Fri May 31 11:54:32 2019 +0100
@@ -16,11 +16,6 @@
 #ifndef FRAMER_H
 #define FRAMER_H
 
-//#include <io.h>
-#include <fcntl.h>
-#include <stdio.h>
-
-
 class Framer  
 {
 public:
@@ -34,19 +29,18 @@
     void resetCounters();
 
 private:
+    unsigned long       m_ulSampleLen;          // DataLength (samples)
+    unsigned int        m_framesRead;           // Read Frames Index
 
-    unsigned long	m_ulSampleLen;		// DataLength (samples)
-    unsigned int	m_framesRead;		// Read Frames Index
+    double*             m_srcBuffer;
+    double*             m_dataFrame;            // Analysis Frame Buffer
+    double*             m_strideFrame;          // Stride Frame Buffer
+    unsigned int        m_frameLength;          // Analysis Frame Length
+    unsigned int        m_stepSize;             // Analysis Frame Stride
 
-    double*			m_srcBuffer;
-    double*			m_dataFrame;		// Analysis Frame Buffer
-    double*			m_strideFrame;		// Stride Frame Buffer
-    unsigned int	m_frameLength;		// Analysis Frame Length
-    unsigned int	m_stepSize;		// Analysis Frame Stride
+    unsigned int        m_maxFrames;
 
-    unsigned int	m_maxFrames;
-
-    unsigned long	m_ulSrcIndex;
+    unsigned long       m_ulSrcIndex;
 };
 
 #endif
--- a/dsp/transforms/DCT.cpp	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/transforms/DCT.cpp	Fri May 31 11:54:32 2019 +0100
@@ -36,14 +36,14 @@
 DCT::forward(const double *in, double *out)
 {
     for (int i = 0; i < m_n; ++i) {
-	m_time2[i*2 + 1] = in[i];
-	m_time2[m_n*4 - i*2 - 1] = in[i];
+        m_time2[i*2 + 1] = in[i];
+        m_time2[m_n*4 - i*2 - 1] = in[i];
     }
 
     m_fft.forward(m_time2.data(), m_freq2r.data(), m_freq2i.data());
 
     for (int i = 0; i < m_n; ++i) {
-	out[i] = m_freq2r[i];
+        out[i] = m_freq2r[i];
     }
 }
 
@@ -52,7 +52,7 @@
 {
     forward(in, out);
     for (int i = 0; i < m_n; ++i) {
-	out[i] /= m_scale;
+        out[i] /= m_scale;
     }
     out[0] /= sqrt(2.0);
 }
@@ -61,21 +61,21 @@
 DCT::inverse(const double *in, double *out)
 {
     for (int i = 0; i < m_n; ++i) {
-	m_freq2r[i] = in[i];
+        m_freq2r[i] = in[i];
     }
     for (int i = 0; i < m_n; ++i) {
-	m_freq2r[m_n*2 - i] = -in[i];
+        m_freq2r[m_n*2 - i] = -in[i];
     }
     m_freq2r[m_n] = 0.0;
 
     for (int i = 0; i <= m_n*2; ++i) {
-	m_freq2i[i] = 0.0;
+        m_freq2i[i] = 0.0;
     }
     
     m_fft.inverse(m_freq2r.data(), m_freq2i.data(), m_time2.data());
 
     for (int i = 0; i < m_n; ++i) {
-	out[i] = m_time2[i*2 + 1];
+        out[i] = m_time2[i*2 + 1];
     }
 }
 
@@ -83,7 +83,7 @@
 DCT::inverseUnitary(const double *in, double *out)
 {
     for (int i = 0; i < m_n; ++i) {
-	m_scaled[i] = in[i] * m_scale;
+        m_scaled[i] = in[i] * m_scale;
     }
     m_scaled[0] *= sqrt(2.0);
     inverse(m_scaled.data(), out);
--- a/dsp/wavelet/Wavelet.cpp	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/wavelet/Wavelet.cpp	Fri May 31 11:54:32 2019 +0100
@@ -33,37 +33,37 @@
         case Daubechies_10: return "Daubechies 10";
         case Daubechies_20: return "Daubechies 20";
         case Daubechies_40: return "Daubechies 40";
-	case Symlet_2: return "Symlet 2";
-	case Symlet_3: return "Symlet 3";
-	case Symlet_4: return "Symlet 4";
-	case Symlet_5: return "Symlet 5";
-	case Symlet_6: return "Symlet 6";
-	case Symlet_7: return "Symlet 7";
-	case Symlet_8: return "Symlet 8";
-	case Symlet_9: return "Symlet 9";
-	case Symlet_10: return "Symlet 10";
-	case Symlet_20: return "Symlet 20";
-	case Symlet_30: return "Symlet 30";
-	case Coiflet_1: return "Coiflet 1";
-	case Coiflet_2: return "Coiflet 2";
-	case Coiflet_3: return "Coiflet 3";
-	case Coiflet_4: return "Coiflet 4";
-	case Coiflet_5: return "Coiflet 5";
-	case Biorthogonal_1_3: return "Biorthogonal 1.3";
-	case Biorthogonal_1_5: return "Biorthogonal 1.5";
-	case Biorthogonal_2_2: return "Biorthogonal 2.2";
-	case Biorthogonal_2_4: return "Biorthogonal 2.4";
-	case Biorthogonal_2_6: return "Biorthogonal 2.6";
-	case Biorthogonal_2_8: return "Biorthogonal 2.8";
-	case Biorthogonal_3_1: return "Biorthogonal 3.1";
-	case Biorthogonal_3_3: return "Biorthogonal 3.3";
-	case Biorthogonal_3_5: return "Biorthogonal 3.5";
-	case Biorthogonal_3_7: return "Biorthogonal 3.7";
-	case Biorthogonal_3_9: return "Biorthogonal 3.9";
-	case Biorthogonal_4_4: return "Biorthogonal 4.4";
-	case Biorthogonal_5_5: return "Biorthogonal 5.5";
-	case Biorthogonal_6_8: return "Biorthogonal 6.8";
-	case Meyer: return "Meyer";
+        case Symlet_2: return "Symlet 2";
+        case Symlet_3: return "Symlet 3";
+        case Symlet_4: return "Symlet 4";
+        case Symlet_5: return "Symlet 5";
+        case Symlet_6: return "Symlet 6";
+        case Symlet_7: return "Symlet 7";
+        case Symlet_8: return "Symlet 8";
+        case Symlet_9: return "Symlet 9";
+        case Symlet_10: return "Symlet 10";
+        case Symlet_20: return "Symlet 20";
+        case Symlet_30: return "Symlet 30";
+        case Coiflet_1: return "Coiflet 1";
+        case Coiflet_2: return "Coiflet 2";
+        case Coiflet_3: return "Coiflet 3";
+        case Coiflet_4: return "Coiflet 4";
+        case Coiflet_5: return "Coiflet 5";
+        case Biorthogonal_1_3: return "Biorthogonal 1.3";
+        case Biorthogonal_1_5: return "Biorthogonal 1.5";
+        case Biorthogonal_2_2: return "Biorthogonal 2.2";
+        case Biorthogonal_2_4: return "Biorthogonal 2.4";
+        case Biorthogonal_2_6: return "Biorthogonal 2.6";
+        case Biorthogonal_2_8: return "Biorthogonal 2.8";
+        case Biorthogonal_3_1: return "Biorthogonal 3.1";
+        case Biorthogonal_3_3: return "Biorthogonal 3.3";
+        case Biorthogonal_3_5: return "Biorthogonal 3.5";
+        case Biorthogonal_3_7: return "Biorthogonal 3.7";
+        case Biorthogonal_3_9: return "Biorthogonal 3.9";
+        case Biorthogonal_4_4: return "Biorthogonal 4.4";
+        case Biorthogonal_5_5: return "Biorthogonal 5.5";
+        case Biorthogonal_6_8: return "Biorthogonal 6.8";
+        case Meyer: return "Meyer";
     }
 
     return "(unknown)";
@@ -78,7 +78,7 @@
     hpd.clear();
 
     int flength = 0;
-	
+        
     switch (wavelet) {
 
     case Haar: 
@@ -99,7 +99,7 @@
         hpd.push_back(-0.22414386804186);
         hpd.push_back(-0.12940952255092);
         flength = 4;
-        break;		
+        break;          
 
     case Daubechies_3:
         lpd.push_back(0.03522629188210);
@@ -588,7 +588,7 @@
         hpd.push_back(-0.00000000000000);
         flength = 80;
         break;
-			
+                        
     case Symlet_2:
         lpd.push_back(-0.12940952255092);
         lpd.push_back(0.22414386804186);
@@ -688,7 +688,7 @@
         hpd.push_back(0.01540410932703);
         flength = 12;
         break;
-			
+                        
     case Symlet_7:
         lpd.push_back(0.00268181456826);
         lpd.push_back(-0.00104738488868);
--- a/dsp/wavelet/Wavelet.h	Fri May 31 11:02:28 2019 +0100
+++ b/dsp/wavelet/Wavelet.h	Fri May 31 11:54:32 2019 +0100
@@ -1,5 +1,4 @@
 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
-
 /*
     QM DSP Library
 
@@ -35,37 +34,37 @@
         Daubechies_10,
         Daubechies_20,
         Daubechies_40,
-	Symlet_2,
-	Symlet_3,
-	Symlet_4,
-	Symlet_5,
-	Symlet_6,
-	Symlet_7,
-	Symlet_8,
-	Symlet_9,
-	Symlet_10,
-	Symlet_20,
-	Symlet_30,
-	Coiflet_1,
-	Coiflet_2,
-	Coiflet_3,
-	Coiflet_4,
-	Coiflet_5,
-	Biorthogonal_1_3,
-	Biorthogonal_1_5,
-	Biorthogonal_2_2,
-	Biorthogonal_2_4,
-	Biorthogonal_2_6,
-	Biorthogonal_2_8,
-	Biorthogonal_3_1,
-	Biorthogonal_3_3,
-	Biorthogonal_3_5,
-	Biorthogonal_3_7,
-	Biorthogonal_3_9,
-	Biorthogonal_4_4,
-	Biorthogonal_5_5,
-	Biorthogonal_6_8,
-	Meyer,
+        Symlet_2,
+        Symlet_3,
+        Symlet_4,
+        Symlet_5,
+        Symlet_6,
+        Symlet_7,
+        Symlet_8,
+        Symlet_9,
+        Symlet_10,
+        Symlet_20,
+        Symlet_30,
+        Coiflet_1,
+        Coiflet_2,
+        Coiflet_3,
+        Coiflet_4,
+        Coiflet_5,
+        Biorthogonal_1_3,
+        Biorthogonal_1_5,
+        Biorthogonal_2_2,
+        Biorthogonal_2_4,
+        Biorthogonal_2_6,
+        Biorthogonal_2_8,
+        Biorthogonal_3_1,
+        Biorthogonal_3_3,
+        Biorthogonal_3_5,
+        Biorthogonal_3_7,
+        Biorthogonal_3_9,
+        Biorthogonal_4_4,
+        Biorthogonal_5_5,
+        Biorthogonal_6_8,
+        Meyer,
 
         LastType = Meyer
     };
--- a/hmm/hmm.c	Fri May 31 11:02:28 2019 +0100
+++ b/hmm/hmm.c	Fri May 31 11:54:32 2019 +0100
@@ -16,9 +16,9 @@
 #include <math.h>
 #include <stdlib.h>
 #include <float.h>
-#include <time.h>				/* to seed random number generator */
+#include <time.h>               /* to seed random number generator */
 
-#include <clapack.h>		/* LAPACK for matrix inversion */
+#include <clapack.h>            /* LAPACK for matrix inversion */
 
 #include "maths/nan-inf.h"
 
@@ -37,788 +37,652 @@
 #ifdef _MAC_OS_X
 #include <vecLib/cblas.h>
 #else
-#include <cblas.h>		/* BLAS for matrix multiplication */
+#include <cblas.h>              /* BLAS for matrix multiplication */
 #endif
 
 #include "hmm.h"
 
 model_t* hmm_init(double** x, int T, int L, int N)
 {
-	int i, j, d, e, t;
-	double s, ss;
-	
-	model_t* model;
-	model = (model_t*) malloc(sizeof(model_t));
-	model->N = N;
-	model->L = L;	
-	model->p0 = (double*) malloc(N*sizeof(double));
-	model->a = (double**) malloc(N*sizeof(double*));
-	model->mu = (double**) malloc(N*sizeof(double*));
-	for (i = 0; i < N; i++)
-	{
-		model->a[i] = (double*) malloc(N*sizeof(double));
-		model->mu[i] = (double*) malloc(L*sizeof(double));
-	}
-	model->cov = (double**) malloc(L*sizeof(double*));
-	for (i = 0; i < L; i++)
-		model->cov[i] = (double*) malloc(L*sizeof(double));
-	
-	srand(time(0));
-	double* global_mean = (double*) malloc(L*sizeof(double));
-	
-	/* find global mean */
-	for (d = 0; d < L; d++)
-	{
-		global_mean[d] = 0;
-		for (t = 0; t < T; t++)
-			global_mean[d] += x[t][d];
-		global_mean[d] /= T;
-	}
-	
-	/* calculate global diagonal covariance */
-	for (d = 0; d < L; d++)
-	{
-		for (e = 0; e < L; e++)
-			model->cov[d][e] = 0;
-		for (t = 0; t < T; t++)
-			model->cov[d][d] += (x[t][d] - global_mean[d]) * (x[t][d] - global_mean[d]);
-		model->cov[d][d] /= T-1;
-	}
-	
-	/* set all means close to global mean */
-	for (i = 0; i < N; i++)
-	{
-		for (d = 0; d < L; d++)
-		{
-			/* add some random noise related to covariance */
-			/* ideally the random number would be Gaussian(0,1), as a hack we make it uniform on [-0.25,0.25] */
-			model->mu[i][d] = global_mean[d] + (0.5 * rand() / (double) RAND_MAX - 0.25) * sqrt(model->cov[d][d]);
-		}
-	}	
-	
-	/* random intial and transition probs */
-	s = 0;
-	for (i = 0; i < N; i++)
-	{
-		model->p0[i] = 1 + rand() / (double) RAND_MAX;
-		s += model->p0[i];
-		ss = 0;
-		for (j = 0; j < N; j++)
-		{
-			model->a[i][j] = 1 + rand() / (double) RAND_MAX;
-			ss += model->a[i][j];
-		}
-		for (j = 0; j < N; j++)
-		{
-			model->a[i][j] /= ss;
-		}
-	}
-	for (i = 0; i < N; i++)
-		model->p0[i] /= s;
-	
-	free(global_mean);
-	
-	return model;
+    int i, j, d, e, t;
+    double s, ss;
+        
+    model_t* model;
+    model = (model_t*) malloc(sizeof(model_t));
+    model->N = N;
+    model->L = L;   
+    model->p0 = (double*) malloc(N*sizeof(double));
+    model->a = (double**) malloc(N*sizeof(double*));
+    model->mu = (double**) malloc(N*sizeof(double*));
+    for (i = 0; i < N; i++) {
+        model->a[i] = (double*) malloc(N*sizeof(double));
+        model->mu[i] = (double*) malloc(L*sizeof(double));
+    }
+    model->cov = (double**) malloc(L*sizeof(double*));
+    for (i = 0; i < L; i++) {
+        model->cov[i] = (double*) malloc(L*sizeof(double));
+    }
+        
+    srand(time(0));
+    double* global_mean = (double*) malloc(L*sizeof(double));
+        
+    /* find global mean */
+    for (d = 0; d < L; d++) {
+        global_mean[d] = 0;
+        for (t = 0; t < T; t++) {
+            global_mean[d] += x[t][d];
+        }
+        global_mean[d] /= T;
+    }
+        
+    /* calculate global diagonal covariance */
+    for (d = 0; d < L; d++) {
+        for (e = 0; e < L; e++) {
+            model->cov[d][e] = 0;
+        }
+        for (t = 0; t < T; t++) {
+            model->cov[d][d] +=
+                (x[t][d] - global_mean[d]) * (x[t][d] - global_mean[d]);
+        }
+        model->cov[d][d] /= T-1;
+    }
+        
+    /* set all means close to global mean */
+    for (i = 0; i < N; i++) {
+        for (d = 0; d < L; d++) {
+            /* add some random noise related to covariance */
+            /* ideally the random number would be Gaussian(0,1),
+               as a hack we make it uniform on [-0.25,0.25] */
+            model->mu[i][d] = global_mean[d] +
+                (0.5 * rand() / (double) RAND_MAX - 0.25)
+                * sqrt(model->cov[d][d]);
+        }
+    }       
+        
+    /* random initial and transition probs */
+    s = 0;
+    for (i = 0; i < N; i++) {
+        model->p0[i] = 1 + rand() / (double) RAND_MAX;
+        s += model->p0[i];
+        ss = 0;
+        for (j = 0; j < N; j++) {
+            model->a[i][j] = 1 + rand() / (double) RAND_MAX;
+            ss += model->a[i][j];
+        }
+        for (j = 0; j < N; j++) {
+            model->a[i][j] /= ss;
+        }
+    }
+    for (i = 0; i < N; i++) {
+        model->p0[i] /= s;
+    }
+        
+    free(global_mean);
+        
+    return model;
 }
 
 void hmm_close(model_t* model)
 {
-	int i;
-	
-	for (i = 0; i < model->N; i++)
-	{
-		free(model->a[i]);
-		free(model->mu[i]);
-	}
-	free(model->a);
-	free(model->mu);
-	for (i = 0; i < model->L; i++)
-		free(model->cov[i]);
-	free(model->cov);	  
-	free(model);	
+    int i;
+        
+    for (i = 0; i < model->N; i++) {
+        free(model->a[i]);
+        free(model->mu[i]);
+    }
+    free(model->a);
+    free(model->mu);
+    for (i = 0; i < model->L; i++) {
+        free(model->cov[i]);
+    }
+    free(model->cov);         
+    free(model);    
 }
 
 void hmm_train(double** x, int T, model_t* model)
 {
-	int i, t;
-	double loglik;	/* overall log-likelihood at each iteration */
-	
-	int N = model->N;
-	int L = model->L;
-	double* p0 = model->p0;
-	double** a = model->a;
-	double** mu = model->mu;
-	double** cov = model->cov;
-	
-	/* allocate memory */	
-	double** gamma = (double**) malloc(T*sizeof(double*));
-	double*** xi = (double***) malloc(T*sizeof(double**));
-	for (t = 0; t < T; t++)
-	{
-		gamma[t] = (double*) malloc(N*sizeof(double));
-		xi[t] = (double**) malloc(N*sizeof(double*));
-		for (i = 0; i < N; i++)
-			xi[t][i] = (double*) malloc(N*sizeof(double));
-	}
-	
-	/* temporary memory */
-	double* gauss_y = (double*) malloc(L*sizeof(double));
-	double* gauss_z = (double*) malloc(L*sizeof(double));	
-			
-	/* obs probs P(j|{x}) */
-	double** b = (double**) malloc(T*sizeof(double*));
-	for (t = 0; t < T; t++)
-		b[t] = (double*) malloc(N*sizeof(double));
-	
-	/* inverse covariance and its determinant */
-	double** icov = (double**) malloc(L*sizeof(double*));
-	for (i = 0; i < L; i++)
-		icov[i] = (double*) malloc(L*sizeof(double));
-	double detcov;
-	
-	double thresh = 0.0001;
-	int niter = 50;	
-	int iter = 0;
-	double loglik1, loglik2;
-	int foundnan = 0;
+    int i, t;
+    double loglik;  /* overall log-likelihood at each iteration */
+        
+    int N = model->N;
+    int L = model->L;
+    double* p0 = model->p0;
+    double** a = model->a;
+    double** mu = model->mu;
+    double** cov = model->cov;
+        
+    /* allocate memory */   
+    double** gamma = (double**) malloc(T*sizeof(double*));
+    double*** xi = (double***) malloc(T*sizeof(double**));
+    for (t = 0; t < T; t++) {
+        gamma[t] = (double*) malloc(N*sizeof(double));
+        xi[t] = (double**) malloc(N*sizeof(double*));
+        for (i = 0; i < N; i++) {
+            xi[t][i] = (double*) malloc(N*sizeof(double));
+        }
+    }
+        
+    /* temporary memory */
+    double* gauss_y = (double*) malloc(L*sizeof(double));
+    double* gauss_z = (double*) malloc(L*sizeof(double));   
+                        
+    /* obs probs P(j|{x}) */
+    double** b = (double**) malloc(T*sizeof(double*));
+    for (t = 0; t < T; t++) {
+        b[t] = (double*) malloc(N*sizeof(double));
+    }
+        
+    /* inverse covariance and its determinant */
+    double** icov = (double**) malloc(L*sizeof(double*));
+    for (i = 0; i < L; i++) {
+        icov[i] = (double*) malloc(L*sizeof(double));
+    }
+    double detcov;
+        
+    double thresh = 0.0001;
+    int niter = 50; 
+    int iter = 0;
+    double loglik1, loglik2;
+    int foundnan = 0;
 
-	while (iter < niter && !foundnan && !(iter > 1 && (loglik - loglik1) < thresh * (loglik1 - loglik2)))	
-	{
-		++iter;
-/*		
-		fprintf(stderr, "calculating obsprobs...\n");
-		fflush(stderr);
-*/		
-		/* precalculate obs probs */
-		invert(cov, L, icov, &detcov);
-		
-		for (t = 0; t < T; t++)
-		{
-			//int allzero = 1;
-			for (i = 0; i < N; i++)
-			{
-				b[t][i] = exp(loggauss(x[t], L, mu[i], icov, detcov, gauss_y, gauss_z));
-		
-				//if (b[t][i] != 0)
-				//	allzero = 0;
-			}
-			/*
-			if (allzero)
-			{
-				printf("all the b[t][i] were zero for t = %d, correcting...\n", t);
-				for (i = 0; i < N; i++)
-				{
-					b[t][i] = 0.00001;
-				}
-			}
-			*/
-		}
-/*		
-		fprintf(stderr, "forwards-backwards...\n");
-		fflush(stderr);
-*/		
-		forward_backwards(xi, gamma, &loglik, &loglik1, &loglik2, iter, N, T, p0, a, b);
-/*		
-		fprintf(stderr, "iteration %d: loglik = %f\n", iter, loglik);		
-		fprintf(stderr, "re-estimation...\n");
-		fflush(stderr);
-*/
-		if (ISNAN(loglik)) {
-		    foundnan = 1;
-		    continue;
-		}
-		
-		baum_welch(p0, a, mu, cov, N, T, L, x, xi, gamma);
-			
-		/*
-		printf("a:\n");
-		for (i = 0; i < model->N; i++)
-		{
-			for (j = 0; j < model->N; j++)
-				printf("%f ", model->a[i][j]);
-			printf("\n");
-		}
-		printf("\n\n");
-		 */
-		//hmm_print(model);
-	}
-	
-	/* deallocate memory */
-	for (t = 0; t < T; t++)
-	{
-		free(gamma[t]);
-		free(b[t]);
-		for (i = 0; i < N; i++)
-			free(xi[t][i]);
-		free(xi[t]);
-	}
-	free(gamma);
-	free(xi);
-	free(b);	
-	
-	for (i = 0; i < L; i++)
-		free(icov[i]);
-	free(icov);
-	
-	free(gauss_y);
-	free(gauss_z);
+    while (iter < niter && !foundnan &&
+           !(iter > 1 && (loglik - loglik1) < thresh * (loglik1 - loglik2))) {
+        
+        ++iter;
+
+        /* precalculate obs probs */
+        invert(cov, L, icov, &detcov);
+                
+        for (t = 0; t < T; t++) {
+            for (i = 0; i < N; i++) {
+                b[t][i] = exp(loggauss(x[t], L, mu[i], icov, detcov, gauss_y, gauss_z));
+            }
+        }
+        forward_backwards(xi, gamma, &loglik, &loglik1, &loglik2,
+                          iter, N, T, p0, a, b);
+        if (ISNAN(loglik)) {
+            foundnan = 1;
+            continue;
+        }
+                
+        baum_welch(p0, a, mu, cov, N, T, L, x, xi, gamma);
+    }
+        
+    /* deallocate memory */
+    for (t = 0; t < T; t++) {
+        free(gamma[t]);
+        free(b[t]);
+        for (i = 0; i < N; i++) {
+            free(xi[t][i]);
+        }
+        free(xi[t]);
+    }
+    free(gamma);
+    free(xi);
+    free(b);        
+        
+    for (i = 0; i < L; i++) {
+        free(icov[i]);
+    }
+    free(icov);
+        
+    free(gauss_y);
+    free(gauss_z);
 }
 
-void baum_welch(double* p0, double** a, double** mu, double** cov, int N, int T, int L, double** x, double*** xi, double** gamma)
+void baum_welch(double* p0, double** a, double** mu, double** cov,
+                int N, int T, int L, double** x, double*** xi, double** gamma)
 {
-	int i, j, t;
-	
-	double* sum_gamma = (double*) malloc(N*sizeof(double));
-	
-	/* temporary memory */
-	double* u = (double*) malloc(L*L*sizeof(double));
-	double* yy = (double*) malloc(T*L*sizeof(double));
-	double* yy2 = (double*) malloc(T*L*sizeof(double));	
-	
-	/* re-estimate transition probs */
-	for (i = 0; i < N; i++)
-	{
-		sum_gamma[i] = 0;
-		for (t = 0; t < T-1; t++)
-			sum_gamma[i] += gamma[t][i];
-	}
-	
-	for (i = 0; i < N; i++)
-	{
-		if (sum_gamma[i] == 0)
-		{
-/*			fprintf(stderr, "sum_gamma[%d] was zero...\n", i); */
-		}
-		//double s = 0;
-		for (j = 0; j < N; j++)
-		{
-			a[i][j] = 0;
-			if (sum_gamma[i] == 0.) continue;
-			for (t = 0; t < T-1; t++)
-				a[i][j] += xi[t][i][j];
-			//s += a[i][j];
-			a[i][j] /= sum_gamma[i];	
-		}
-		/*
-		 for (j = 0; j < N; j++)
-		 {
-			 a[i][j] /= s;
-		 }
-		 */
-	}
-	
-	/* NB: now we need to sum gamma over all t */
-	for (i = 0; i < N; i++)
-		sum_gamma[i] += gamma[T-1][i];
-	
-	/* re-estimate initial probs */
-	for (i = 0; i < N; i++)
-		p0[i] = gamma[0][i];
-	
-	/* re-estimate covariance */
-	int d, e;
-	double sum_sum_gamma = 0;
-	for (i = 0; i < N; i++)
-		sum_sum_gamma += sum_gamma[i];		
-	
-	/*
-	 for (d = 0; d < L; d++)
-	 {
-		 for (e = d; e < L; e++)
-		 {
-			 cov[d][e] = 0;
-			 for (t = 0; t < T; t++)
-				 for (j = 0; j < N; j++)
-					 cov[d][e] += gamma[t][j] * (x[t][d] - mu[j][d]) * (x[t][e] - mu[j][e]);
-			 
-			 cov[d][e] /= sum_sum_gamma;
-			 
-			 if (ISNAN(cov[d][e]))
-			 {
-				 printf("cov[%d][%d] was nan\n", d, e);
-				 for (j = 0; j < N; j++)
-					 for (i = 0; i < L; i++)
-						 if (ISNAN(mu[j][i]))
-							 printf("mu[%d][%d] was nan\n", j, i);
-				 for (t = 0; t < T; t++)
-					 for (j = 0; j < N; j++)
-						 if (ISNAN(gamma[t][j]))
-							 printf("gamma[%d][%d] was nan\n", t, j);
-				 exit(-1);
-			 }
-		 }
-	 }
-	 for (d = 0; d < L; d++)
-	 for (e = 0; e < d; e++)
-	 cov[d][e] = cov[e][d];
-	 */
-	
-	/* using BLAS */
-	for (d = 0; d < L; d++)
-		for (e = 0; e < L; e++)
-			cov[d][e] = 0;
-	
-	for (j = 0; j < N; j++)
-	{
-		for (d = 0; d < L; d++)
-			for (t = 0; t < T; t++)
-			{
-				yy[d*T+t] = x[t][d] - mu[j][d];
-				yy2[d*T+t] = gamma[t][j] * (x[t][d] - mu[j][d]);
-			}
-				
-				cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans, L, L, T, 1.0, yy, T, yy2, T, 0, u, L);
-		
-		for (e = 0; e < L; e++)
-			for (d = 0; d < L; d++)
-				cov[d][e] += u[e*L+d];
-	}
-	
-	for (d = 0; d < L; d++)
-		for (e = 0; e < L; e++)
-			cov[d][e] /= T; /* sum_sum_gamma; */			
-	
-	//printf("sum_sum_gamma = %f\n", sum_sum_gamma); /* fine, = T IS THIS ALWAYS TRUE with pooled cov?? */
-	
-	/* re-estimate means */
-	for (j = 0; j < N; j++)
-	{
-		for (d = 0; d < L; d++)
-		{
-			mu[j][d] = 0;
-			for (t = 0; t < T; t++)
-				mu[j][d] += gamma[t][j] * x[t][d];
-			mu[j][d] /= sum_gamma[j];
-		}
-	}
-	
-	/* deallocate memory */
-	free(sum_gamma);
-	free(yy);
-	free(yy2);
-	free(u);
+    int i, j, t;
+        
+    double* sum_gamma = (double*) malloc(N*sizeof(double));
+        
+    /* temporary memory */
+    double* u = (double*) malloc(L*L*sizeof(double));
+    double* yy = (double*) malloc(T*L*sizeof(double));
+    double* yy2 = (double*) malloc(T*L*sizeof(double));     
+        
+    /* re-estimate transition probs */
+    for (i = 0; i < N; i++) {
+        sum_gamma[i] = 0;
+        for (t = 0; t < T-1; t++) {
+            sum_gamma[i] += gamma[t][i];
+        }
+    }
+        
+    for (i = 0; i < N; i++) {
+        for (j = 0; j < N; j++) {
+            a[i][j] = 0;
+            if (sum_gamma[i] == 0.) {
+                continue;
+            }
+            for (t = 0; t < T-1; t++) {
+                a[i][j] += xi[t][i][j];
+            }
+            a[i][j] /= sum_gamma[i];        
+        }
+    }
+        
+    /* NB: now we need to sum gamma over all t */
+    for (i = 0; i < N; i++) {
+        sum_gamma[i] += gamma[T-1][i];
+    }
+        
+    /* re-estimate initial probs */
+    for (i = 0; i < N; i++) {
+        p0[i] = gamma[0][i];
+    }
+        
+    /* re-estimate covariance */
+    int d, e;
+    double sum_sum_gamma = 0;
+    for (i = 0; i < N; i++) {
+        sum_sum_gamma += sum_gamma[i];
+    }
+        
+    /* using BLAS */
+    for (d = 0; d < L; d++) {
+        for (e = 0; e < L; e++) {
+            cov[d][e] = 0;
+        }
+    }
+        
+    for (j = 0; j < N; j++) {
+        
+        for (d = 0; d < L; d++) {
+            for (t = 0; t < T; t++) {
+                yy[d*T+t] = x[t][d] - mu[j][d];
+                yy2[d*T+t] = gamma[t][j] * (x[t][d] - mu[j][d]);
+            }
+        }
+                                
+        cblas_dgemm(CblasColMajor, CblasTrans, CblasNoTrans,
+                    L, L, T, 1.0, yy, T, yy2, T, 0, u, L);
+                
+        for (e = 0; e < L; e++) {
+            for (d = 0; d < L; d++) {
+                cov[d][e] += u[e*L+d];
+            }
+        }
+    }
+        
+    for (d = 0; d < L; d++) {
+        for (e = 0; e < L; e++) {
+            cov[d][e] /= T; /* sum_sum_gamma; */
+        }
+    }
+    
+    //printf("sum_sum_gamma = %f\n", sum_sum_gamma); /* fine, = T IS THIS ALWAYS TRUE with pooled cov?? */
+        
+    /* re-estimate means */
+    for (j = 0; j < N; j++) {
+        for (d = 0; d < L; d++) {
+            mu[j][d] = 0;
+            for (t = 0; t < T; t++)
+                mu[j][d] += gamma[t][j] * x[t][d];
+            mu[j][d] /= sum_gamma[j];
+        }
+    }
+        
+    /* deallocate memory */
+    free(sum_gamma);
+    free(yy);
+    free(yy2);
+    free(u);
 }
 
-void forward_backwards(double*** xi, double** gamma, double* loglik, double* loglik1, double* loglik2, int iter, int N, int T, double* p0, double** a, double** b)
+void forward_backwards(double*** xi, double** gamma,
+                       double* loglik, double* loglik1, double* loglik2,
+                       int iter, int N, int T,
+                       double* p0, double** a, double** b)
 {
-	/* forwards-backwards with scaling */
-	int i, j, t;
+    /* forwards-backwards with scaling */
+    int i, j, t;
+        
+    double** alpha = (double**) malloc(T*sizeof(double*));
+    double** beta = (double**) malloc(T*sizeof(double*));
+    for (t = 0; t < T; t++) {
+        alpha[t] = (double*) malloc(N*sizeof(double));
+        beta[t] = (double*) malloc(N*sizeof(double));
+    }
+        
+    /* scaling coefficients */
+    double* c = (double*) malloc(T*sizeof(double));
+        
+    /* calculate forward probs and scale coefficients */
+    c[0] = 0;
+    for (i = 0; i < N; i++) {
+        alpha[0][i] = p0[i] * b[0][i];
+        c[0] += alpha[0][i];
+    }
+    c[0] = 1 / c[0];
+    for (i = 0; i < N; i++) {
+        alpha[0][i] *= c[0];            
+    }
+        
+    *loglik1 = *loglik;
+    *loglik = -log(c[0]);
+    if (iter == 2) {
+        *loglik2 = *loglik;
+    }
+        
+    for (t = 1; t < T; t++) {
+        
+        c[t] = 0;
+
+        for (j = 0; j < N; j++) {
+            alpha[t][j] = 0;
+            for (i = 0; i < N; i++) {
+                alpha[t][j] += alpha[t-1][i] * a[i][j];
+            }
+            alpha[t][j] *= b[t][j];
+                        
+            c[t] += alpha[t][j];
+        }
+                
+        c[t] = 1 / c[t];
+        for (j = 0; j < N; j++) {
+            alpha[t][j] *= c[t];
+        }
+                
+        *loglik -= log(c[t]);
+    }
+        
+    /* calculate backwards probs using same coefficients */
+    for (i = 0; i < N; i++) {
+        beta[T-1][i] = 1;
+    }
+    t = T - 1;
+
+    while (1) {
+        for (i = 0; i < N; i++) {
+            beta[t][i] *= c[t];
+        }
+                
+        if (t == 0) {
+            break;
+        }
+                
+        for (i = 0; i < N; i++) {
+            beta[t-1][i] = 0;
+            for (j = 0; j < N; j++) {
+                beta[t-1][i] += a[i][j] * b[t][j] * beta[t][j];
+            }
+        }
+                
+        t--;
+    }
 	
-	double** alpha = (double**) malloc(T*sizeof(double*));
-	double** beta = (double**) malloc(T*sizeof(double*));
-	for (t = 0; t < T; t++)
-	{
-		alpha[t] = (double*) malloc(N*sizeof(double));
-		beta[t] = (double*) malloc(N*sizeof(double));
-	}
+    /* calculate posterior probs */
+    double tot;
+    for (t = 0; t < T; t++) {
+        tot = 0;
+        for (i = 0; i < N; i++) {
+            gamma[t][i] = alpha[t][i] * beta[t][i];
+            tot += gamma[t][i];
+        }
+        for (i = 0; i < N; i++) {
+            gamma[t][i] /= tot;				
+        }
+    }		
 	
-	/* scaling coefficients */
-	double* c = (double*) malloc(T*sizeof(double));
+    for (t = 0; t < T-1; t++) {
+        tot = 0;
+        for (i = 0; i < N; i++) {
+            for (j = 0; j < N; j++) {
+                xi[t][i][j] = alpha[t][i] * a[i][j] * b[t+1][j] * beta[t+1][j];
+                tot += xi[t][i][j];
+            }
+        }
+        for (i = 0; i < N; i++) {
+            for (j = 0; j < N; j++) {
+                xi[t][i][j] /= tot;
+            }
+        }
+    }
 	
-	/* calculate forward probs and scale coefficients */
-	c[0] = 0;
-	for (i = 0; i < N; i++)
-	{
-		alpha[0][i] = p0[i] * b[0][i];
-		c[0] += alpha[0][i];
-		
-		//printf("p0[%d] = %f, b[0][%d] = %f\n", i, p0[i], i, b[0][i]);
-	}
-	c[0] = 1 / c[0];
-	for (i = 0; i < N; i++)
-	{
-		alpha[0][i] *= c[0];		
-		
-		//printf("alpha[0][%d] = %f\n", i, alpha[0][i]);	/* OK agrees with Matlab */
-	}
-	
-	*loglik1 = *loglik;
-	*loglik = -log(c[0]);
-	if (iter == 2)
-		*loglik2 = *loglik;
-	
-	for (t = 1; t < T; t++)
-	{			
-		c[t] = 0;
-		for (j = 0; j < N; j++)
-		{
-			alpha[t][j] = 0;
-			for (i = 0; i < N; i++)
-				alpha[t][j] += alpha[t-1][i] * a[i][j];
-			alpha[t][j] *= b[t][j];
-			
-			c[t] += alpha[t][j];
-		}
-		
-		/*
-		 if (c[t] == 0)
-		 {
-			 printf("c[%d] = 0, going to blow up so exiting\n", t);
-			 for (i = 0; i < N; i++)
-				 if (b[t][i] == 0)
-					 fprintf(stderr, "b[%d][%d] was zero\n", t, i);
-			 fprintf(stderr, "x[t] was \n");
-			 for (i = 0; i < L; i++)
-				 fprintf(stderr, "%f ", x[t][i]);
-			 fprintf(stderr, "\n\n");
-			 exit(-1);
-		 }
-		 */
-		
-		c[t] = 1 / c[t];
-		for (j = 0; j < N; j++)
-			alpha[t][j] *= c[t];
-		
-		//printf("c[%d] = %e\n", t, c[t]);
-		
-		*loglik -= log(c[t]);
-	}
-	
-	/* calculate backwards probs using same coefficients */
-	for (i = 0; i < N; i++)
-		beta[T-1][i] = 1;
-	t = T - 1;
-	while (1)
-	{
-		for (i = 0; i < N; i++)
-			beta[t][i] *= c[t];
-		
-		if (t == 0)
-			break;
-		
-		for (i = 0; i < N; i++)
-		{
-			beta[t-1][i] = 0;
-			for (j = 0; j < N; j++)
-				beta[t-1][i] += a[i][j] * b[t][j] * beta[t][j];
-		}
-		
-		t--;
-	}
-	
-	/*
-	 printf("alpha:\n");
-	 for (t = 0; t < T; t++)
-	 {
-		 for (i = 0; i < N; i++)
-			 printf("%4.4e\t\t", alpha[t][i]);
-		 printf("\n");
-	 }
-	 printf("\n\n");printf("beta:\n");
-	 for (t = 0; t < T; t++)
-	 {
-		 for (i = 0; i < N; i++)
-			 printf("%4.4e\t\t", beta[t][i]);
-		 printf("\n");
-	 }
-	 printf("\n\n");
-	 */
-	
-	/* calculate posterior probs */
-	double tot;
-	for (t = 0; t < T; t++)
-	{
-		tot = 0;
-		for (i = 0; i < N; i++)
-		{
-			gamma[t][i] = alpha[t][i] * beta[t][i];
-			tot += gamma[t][i];
-		}
-		for (i = 0; i < N; i++)
-		{
-			gamma[t][i] /= tot;				
-			
-			//printf("gamma[%d][%d] = %f\n", t, i, gamma[t][i]);				
-		}
-	}		
-	
-	for (t = 0; t < T-1; t++)
-	{
-		tot = 0;
-		for (i = 0; i < N; i++)
-		{
-			for (j = 0; j < N; j++)
-			{
-				xi[t][i][j] = alpha[t][i] * a[i][j] * b[t+1][j] * beta[t+1][j];
-				tot += xi[t][i][j];
-			}
-		}
-		for (i = 0; i < N; i++)
-			for (j = 0; j < N; j++)
-				xi[t][i][j] /= tot;
-	}
-	
-	/*
-	 // CHECK - fine
-	 // gamma[t][i] = \sum_j{xi[t][i][j]}
-	 tot = 0;
-	 for (j = 0; j < N; j++)
-	 tot += xi[3][1][j];
-	 printf("gamma[3][1] = %f, sum_j(xi[3][1][j]) = %f\n", gamma[3][1], tot);
-	 */	
-	
-	for (t = 0; t < T; t++)
-	{
-		free(alpha[t]);
-		free(beta[t]);
-	}
-	free(alpha);
-	free(beta);
-	free(c);
+    for (t = 0; t < T; t++) {
+        free(alpha[t]);
+        free(beta[t]);
+    }
+    free(alpha);
+    free(beta);
+    free(c);
 }
 
 void viterbi_decode(double** x, int T, model_t* model, int* q)
 {
-	int i, j, t;
-	double max;
-	int havemax;
+    int i, j, t;
+    double max;
+    int havemax;
 	
-	int N = model->N;
-	int L = model->L;
-	double* p0 = model->p0;
-	double** a = model->a;
-	double** mu = model->mu;
-	double** cov = model->cov;
+    int N = model->N;
+    int L = model->L;
+    double* p0 = model->p0;
+    double** a = model->a;
+    double** mu = model->mu;
+    double** cov = model->cov;
 	
-	/* inverse covariance and its determinant */
-	double** icov = (double**) malloc(L*sizeof(double*));
-	for (i = 0; i < L; i++)
-		icov[i] = (double*) malloc(L*sizeof(double));
-	double detcov;
+    /* inverse covariance and its determinant */
+    double** icov = (double**) malloc(L*sizeof(double*));
+    for (i = 0; i < L; i++) {
+        icov[i] = (double*) malloc(L*sizeof(double));
+    }
+    double detcov;
 	
-	double** logb = (double**) malloc(T*sizeof(double*));
-	double** phi = (double**) malloc(T*sizeof(double*));
-	int** psi = (int**) malloc(T*sizeof(int*));
-	for (t = 0; t < T; t++)
-	{
-		logb[t] = (double*) malloc(N*sizeof(double));
-		phi[t] = (double*) malloc(N*sizeof(double));
-		psi[t] = (int*) malloc(N*sizeof(int));
-	}
+    double** logb = (double**) malloc(T*sizeof(double*));
+    double** phi = (double**) malloc(T*sizeof(double*));
+    int** psi = (int**) malloc(T*sizeof(int*));
+
+    for (t = 0; t < T; t++) {
+        logb[t] = (double*) malloc(N*sizeof(double));
+        phi[t] = (double*) malloc(N*sizeof(double));
+        psi[t] = (int*) malloc(N*sizeof(int));
+    }
 	
-	/* temporary memory */
-	double* gauss_y = (double*) malloc(L*sizeof(double));
-	double* gauss_z = (double*) malloc(L*sizeof(double));	
+    /* temporary memory */
+    double* gauss_y = (double*) malloc(L*sizeof(double));
+    double* gauss_z = (double*) malloc(L*sizeof(double));	
 	
-	/* calculate observation logprobs */
-	invert(cov, L, icov, &detcov);
-	for (t = 0; t < T; t++)
-		for (i = 0; i < N; i++)
-			logb[t][i] = loggauss(x[t], L, mu[i], icov, detcov, gauss_y, gauss_z);
+    /* calculate observation logprobs */
+    invert(cov, L, icov, &detcov);
+    for (t = 0; t < T; t++) {
+        for (i = 0; i < N; i++) {
+            logb[t][i] = loggauss
+                (x[t], L, mu[i], icov, detcov, gauss_y, gauss_z);
+        }
+    }
 	
-	/* initialise */
-	for (i = 0; i < N; i++)
-	{
-		phi[0][i] = log(p0[i]) + logb[0][i];
-		psi[0][i] = 0;
-	}
+    /* initialise */
+    for (i = 0; i < N; i++) {
+        phi[0][i] = log(p0[i]) + logb[0][i];
+        psi[0][i] = 0;
+    }
 	
-	for (t = 1; t < T; t++)
-	{
-		for (j = 0; j < N; j++)
-		{
-			max = -1000000;
-			havemax = 0;
+    for (t = 1; t < T; t++) {
+        for (j = 0; j < N; j++) {
+            max = -1000000;
+            havemax = 0;
 
-			psi[t][j] = 0;
-			for (i = 0; i < N; i++)
-			{
-				if (phi[t-1][i] + log(a[i][j]) > max || !havemax)
-				{
-					max = phi[t-1][i] + log(a[i][j]);
-					phi[t][j] = max + logb[t][j];
-					psi[t][j] = i;
-					havemax = 1;
-				}
-			}
-		}
-	}
+            psi[t][j] = 0;
+            for (i = 0; i < N; i++) {
+                if (phi[t-1][i] + log(a[i][j]) > max || !havemax) {
+                    max = phi[t-1][i] + log(a[i][j]);
+                    phi[t][j] = max + logb[t][j];
+                    psi[t][j] = i;
+                    havemax = 1;
+                }
+            }
+        }
+    }
 	
-	/* find maximising state at time T-1 */
-	max = phi[T-1][0];
-	q[T-1] = 0;
-	for (i = 1; i < N; i++)
-	{
-		if (phi[T-1][i] > max)
-		{
-			max = phi[T-1][i];
-			q[T-1] = i;
-		}
-	}
-
+    /* find maximising state at time T-1 */
+    max = phi[T-1][0];
+    q[T-1] = 0;
+    for (i = 1; i < N; i++) {
+        if (phi[T-1][i] > max) {
+            max = phi[T-1][i];
+            q[T-1] = i;
+        }
+    }
 	
-	/* track back */
-	t = T - 2;
-	while (t >= 0)
-	{
-		q[t] = psi[t+1][q[t+1]];
-		t--;
-	}
+    /* track back */
+    t = T - 2;
+    while (t >= 0) {
+        q[t] = psi[t+1][q[t+1]];
+        t--;
+    }
 	
-	/* de-allocate memory */
-	for (i = 0; i < L; i++)
-		free(icov[i]);
-	free(icov);
-	for (t = 0; t < T; t++)
-	{
-		free(logb[t]);
-		free(phi[t]);
-		free(psi[t]);
-	}
-	free(logb);
-	free(phi);
-	free(psi);
+    /* de-allocate memory */
+    for (i = 0; i < L; i++) {
+        free(icov[i]);
+    }
+    free(icov);
+    for (t = 0; t < T; t++) {
+        free(logb[t]);
+        free(phi[t]);
+        free(psi[t]);
+    }
+    free(logb);
+    free(phi);
+    free(psi);
 	
-	free(gauss_y);
-	free(gauss_z);
+    free(gauss_y);
+    free(gauss_z);
 }
 
 /* invert matrix and calculate determinant using LAPACK */
 void invert(double** cov, int L, double** icov, double* detcov)
 {
-	/* copy square matrix into a vector in column-major order */
-	double* a = (double*) malloc(L*L*sizeof(double));
-	int i, j;
-	for(j=0; j < L; j++)
-		for (i=0; i < L; i++) 
-			a[j*L+i] = cov[i][j];
+    /* copy square matrix into a vector in column-major order */
+    double* a = (double*) malloc(L*L*sizeof(double));
+    int i, j;
+    for (j=0; j < L; j++) {
+        for (i=0; i < L; i++) {
+            a[j*L+i] = cov[i][j];
+        }
+    }
 	
-	int M = (int) L;	
-	int* ipiv = (int *) malloc(L*L*sizeof(int));
-	int ret;
+    int M = (int) L;	
+    int* ipiv = (int *) malloc(L*L*sizeof(int));
+    int ret;
 	
-	/* LU decomposition */
-	ret = dgetrf_(&M, &M, a, &M, ipiv, &ret);	/* ret should be zero, negative if cov is singular */
-	if (ret < 0)
-	{
-		fprintf(stderr, "Covariance matrix was singular, couldn't invert\n");
-		exit(-1);
-	}
+    /* LU decomposition */
+    ret = dgetrf_(&M, &M, a, &M, ipiv, &ret);	/* ret should be zero, negative if cov is singular */
+    if (ret < 0) {
+        fprintf(stderr, "Covariance matrix was singular, couldn't invert\n");
+        exit(-1);
+    }
 	
-	/* find determinant */
-	double det = 1;
-	for(i = 0; i < L; i++)
-		det *= a[i*L+i];
-	// TODO: get this to work!!! If detcov < 0 then cov is bad anyway...
-	/*
-	int sign = 1;
-	for (i = 0; i < L; i++)
-		if (ipiv[i] != i)
-			sign = -sign;
-	det *= sign;
-	 */
-	if (det < 0)
-		det = -det;
-	*detcov = det;
+    /* find determinant */
+    double det = 1;
+    for (i = 0; i < L; i++) {
+        det *= a[i*L+i];
+    }
+
+    // TODO: get this to work!!! If detcov < 0 then cov is bad anyway...
+    if (det < 0) {
+        det = -det;
+    }
+    *detcov = det;
 	
-	/* allocate required working storage */
+    /* allocate required working storage */
 #ifndef HAVE_ATLAS
-	int lwork = -1;
-	double lwbest = 0.0;
-	dgetri_(&M, a, &M, ipiv, &lwbest, &lwork, &ret);
-	lwork = (int) lwbest;	
-	double* work  = (double*) malloc(lwork*sizeof(double));
+    int lwork = -1;
+    double lwbest = 0.0;
+    dgetri_(&M, a, &M, ipiv, &lwbest, &lwork, &ret);
+    lwork = (int) lwbest;	
+    double* work  = (double*) malloc(lwork*sizeof(double));
 #endif
 	
-	/* find inverse */
-	dgetri_(&M, a, &M, ipiv, work, &lwork, &ret);
+    /* find inverse */
+    dgetri_(&M, a, &M, ipiv, work, &lwork, &ret);
 
-	for(j=0; j < L; j++)
-		for (i=0; i < L; i++) 
-			icov[i][j] = a[j*L+i];	
+    for (j=0; j < L; j++) {
+        for (i=0; i < L; i++) {
+            icov[i][j] = a[j*L+i];
+        }
+    }
 	
 #ifndef HAVE_ATLAS	
-	free(work);
+    free(work);
 #endif
-	free(a);	
+    free(a);	
 }
 
 /* probability of multivariate Gaussian given mean, inverse and determinant of covariance */
 double gauss(double* x, int L, double* mu, double** icov, double detcov, double* y, double* z)
 {
-	int i;
-	double s = 0;
-	for (i = 0; i < L; i++)
-		y[i] = x[i] - mu[i];
-	for (i = 0; i < L; i++)
-	{
-		//z[i] = 0;
-		//for (j = 0; j < L; j++)
-		//	z[i] += icov[i][j] *  y[j];
-		z[i] = cblas_ddot(L, &icov[i][0], 1, y, 1);
-	}
-	s = cblas_ddot(L, z, 1, y, 1);
-	//for (i = 0; i < L; i++)
-	//	s += z[i] * y[i];	
+    int i;
+    double s = 0;
+
+    for (i = 0; i < L; i++) {
+        y[i] = x[i] - mu[i];
+    }
+
+    for (i = 0; i < L; i++) {
+        z[i] = cblas_ddot(L, &icov[i][0], 1, y, 1);
+    }
+
+    s = cblas_ddot(L, z, 1, y, 1);
 	
-	return exp(-s/2.0) / (pow(2*PI, L/2.0) * sqrt(detcov));
+    return exp(-s/2.0) / (pow(2*PI, L/2.0) * sqrt(detcov));
 }
 
 /* log probability of multivariate Gaussian given mean, inverse and determinant of covariance */
 double loggauss(double* x, int L, double* mu, double** icov, double detcov, double* y, double* z)
 {
-	int i;
-	double s = 0;
-	double ret;
-	for (i = 0; i < L; i++)
-		y[i] = x[i] - mu[i];
-	for (i = 0; i < L; i++)
-	{
-		//z[i] = 0;
-		//for (j = 0; j < L; j++)
-		//	z[i] += icov[i][j] *  y[j];
-		z[i] = cblas_ddot(L, &icov[i][0], 1, y, 1);
-	}
-	s = cblas_ddot(L, z, 1, y, 1);
-	//for (i = 0; i < L; i++)
-	//	s += z[i] * y[i];	
+    int i;
+    double s = 0;
+    double ret;
+
+    for (i = 0; i < L; i++) {
+        y[i] = x[i] - mu[i];
+    }
+
+    for (i = 0; i < L; i++) {
+        z[i] = cblas_ddot(L, &icov[i][0], 1, y, 1);
+    }
+
+    s = cblas_ddot(L, z, 1, y, 1);
 	
-	ret = -0.5 * (s + L * log(2*PI) + log(detcov));
+    ret = -0.5 * (s + L * log(2*PI) + log(detcov));
 	
-	/*
-	// TEST
-	if (ISINF(ret) > 0)
-		printf("loggauss returning infinity\n");
-	if (ISINF(ret) < 0)
-		printf("loggauss returning -infinity\n");
-	if (ISNAN(ret))
-		printf("loggauss returning nan\n");	
-	*/
-	
-	return ret;
+    return ret;
 }
 
 void hmm_print(model_t* model)
 {
-	int i, j;
-	printf("p0:\n");
-	for (i = 0; i < model->N; i++)
-		printf("%f ", model->p0[i]);
-	printf("\n\n");
-	printf("a:\n");
-	for (i = 0; i < model->N; i++)
-	{
-		for (j = 0; j < model->N; j++)
-			printf("%f ", model->a[i][j]);
-		printf("\n");
-	}
-	printf("\n\n");
-	printf("mu:\n");
-	for (i = 0; i < model->N; i++)
-	{
-		for (j = 0; j < model->L; j++)
-			printf("%f ", model->mu[i][j]);
-		printf("\n");
-	}
-	printf("\n\n");
-	printf("cov:\n");
-	for (i = 0; i < model->L; i++)
-	{
-		for (j = 0; j < model->L; j++)
-			printf("%f ", model->cov[i][j]);
-		printf("\n");
-	}
-	printf("\n\n");
+    int i, j;
+    printf("p0:\n");
+    for (i = 0; i < model->N; i++) {
+        printf("%f ", model->p0[i]);
+    }
+    printf("\n\n");
+    printf("a:\n");
+    for (i = 0; i < model->N; i++) {
+        for (j = 0; j < model->N; j++) {
+            printf("%f ", model->a[i][j]);
+        }
+        printf("\n");
+    }
+    printf("\n\n");
+    printf("mu:\n");
+    for (i = 0; i < model->N; i++) {
+        for (j = 0; j < model->L; j++) {
+            printf("%f ", model->mu[i][j]);
+        }
+        printf("\n");
+    }
+    printf("\n\n");
+    printf("cov:\n");
+    for (i = 0; i < model->L; i++) {
+        for (j = 0; j < model->L; j++) {
+            printf("%f ", model->cov[i][j]);
+        }
+        printf("\n");
+    }
+    printf("\n\n");
 }
 
 
--- a/hmm/hmm.h	Fri May 31 11:02:28 2019 +0100
+++ b/hmm/hmm.h	Fri May 31 11:54:32 2019 +0100
@@ -1,10 +1,3 @@
-#ifndef _HMM_H
-#define _HMM_H
-
-#ifdef __cplusplus
-extern "C" {
-#endif
-
 /*
  *  hmm.h
  *
@@ -19,29 +12,49 @@
  *
  */
 
+#ifndef _HMM_H
+#define _HMM_H
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 #ifndef PI    
 #define PI 3.14159265358979323846264338327950288
 #endif 
 
 typedef struct _model_t {
-	int N;			/* number of states */
-	double* p0;		/* initial probs */
-	double** a;		/* transition probs */
-	int L;			/* dimensionality of data */
-	double** mu;	/* state means */
-	double** cov;	/* covariance, tied between all states */
+    int N;          /* number of states */
+    double* p0;     /* initial probs */
+    double** a;     /* transition probs */
+    int L;          /* dimensionality of data */
+    double** mu;    /* state means */
+    double** cov;   /* covariance, tied between all states */
 } model_t;
 
-void hmm_train(double** x, int T, model_t* model);							/* with scaling */
-void forward_backwards(double*** xi, double** gamma, double* loglik, double* loglik1, double* loglik2, int iter, 
-					   int N, int T, double* p0, double** a, double** b);
-void baum_welch(double* p0, double** a, double** mu, double** cov, int N, int T, int L, double** x, double*** xi, double** gamma);
-void viterbi_decode(double** x, int T, model_t* model, int* q);				/* using logs */
+void hmm_train(double** x, int T, model_t* model); /* with scaling */
+
+void forward_backwards(double*** xi, double** gamma,
+                       double* loglik, double* loglik1, double* loglik2,
+                       int iter, int N, int T,
+                       double* p0, double** a, double** b);
+    
+void baum_welch(double* p0, double** a, double** mu, double** cov,
+                int N, int T, int L, double** x, double*** xi, double** gamma);
+
+void viterbi_decode(double** x, int T, model_t* model, int* q); /* using logs */
+
 model_t* hmm_init(double** x, int T, int L, int N);
 void hmm_close(model_t* model);
-void invert(double** cov, int L, double** icov, double* detcov);			/* uses LAPACK (included with Mac OSX) */
-double gauss(double* x, int L, double* mu, double** icov, double detcov, double* y, double* z);
-double loggauss(double* x, int L, double* mu, double** icov, double detcov, double* y, double* z);
+    
+void invert(double** cov, int L, double** icov, double* detcov); /* uses LAPACK */
+    
+double gauss(double* x, int L, double* mu, double** icov,
+             double detcov, double* y, double* z);
+    
+double loggauss(double* x, int L, double* mu, double** icov,
+                double detcov, double* y, double* z);
+    
 void hmm_print(model_t* model);
 
 #ifdef __cplusplus
--- a/maths/pca/pca.c	Fri May 31 11:02:28 2019 +0100
+++ b/maths/pca/pca.c	Fri May 31 11:54:32 2019 +0100
@@ -32,58 +32,58 @@
 /* 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;
+    double *mean;
+    int i, j, j1, j2;
 
 /* Allocate storage for mean vector */
 
-mean = (double*) malloc(m*sizeof(double));
+    mean = (double*) malloc(m*sizeof(double));
 
 /* Determine mean of column vectors of input data matrix */
 
-for (j = 0; j < m; j++)
+    for (j = 0; j < m; j++)
     {
-    mean[j] = 0.0;
-    for (i = 0; i < n; i++)
+        mean[j] = 0.0;
+        for (i = 0; i < n; i++)
         {
-        mean[j] += data[i][j];
+            mean[j] += data[i][j];
         }
-    mean[j] /= (double)n;
+        mean[j] /= (double)n;
     }
 
 /*
-printf("\nMeans of column vectors:\n");
-for (j = 0; j < m; j++)  {
-    printf("%12.1f",mean[j]);  }   printf("\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 (i = 0; i < n; i++)
     {
-    for (j = 0; j < m; j++)
+        for (j = 0; j < m; j++)
         {
-        data[i][j] -= mean[j];
+            data[i][j] -= mean[j];
         }
     }
 
 /* Calculate the m * m covariance matrix. */
-for (j1 = 0; j1 < m; j1++)
+    for (j1 = 0; j1 < m; j1++)
     {
-    for (j2 = j1; j2 < m; j2++)
+        for (j2 = j1; j2 < m; j2++)
         {
-        symmat[j1][j2] = 0.0;
-        for (i = 0; i < n; i++)
+            symmat[j1][j2] = 0.0;
+            for (i = 0; i < n; i++)
             {
-            symmat[j1][j2] += data[i][j1] * data[i][j2];
+                symmat[j1][j2] += data[i][j1] * data[i][j2];
             }
-        symmat[j2][j1] = symmat[j1][j2];
+            symmat[j2][j1] = symmat[j1][j2];
         }
     }
 
-free(mean);
+    free(mean);
 
-return;
+    return;
 
 }
 
@@ -101,84 +101,84 @@
 /**  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.  */
+   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--)
+    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;
+        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++)
+    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;
+        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;
     }
 }
 
@@ -186,171 +186,168 @@
 
 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++)
+    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))
+        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);
+                        c = g / f;
+                        r = sqrt((c * c) + 1.0);
+                        e[i+1] = f * r;
+                        c *= (s = 1.0/r);
                     }
-					else
+                    else
                     {
-						s = f / g;
-						r = sqrt((s * s) + 1.0);
-						e[i+1] = g * r;
-						s *= (c = 1.0/r);
+                        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);
-	}
+                    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
-		**********************************************************************/
-	
+    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 */
-	//	}
-	//}
+    //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 */
+    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");
-	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");  }
- */
+  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 (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]; }
+            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");  }
- */
+/*      
+        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 (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 */
+//      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");  }
-	*/
+  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);
+    for (i = 0; i < m; i++)
+        free(symmat[i]);
+    free(symmat);
 //FREE_ARRAY(symmat2,m);
-free(evals);
-free(interm);
+    free(evals);
+    free(interm);
 
 }
-
-
-