Mercurial > hg > qm-dsp
changeset 483:fdaa63607c15
Untabify, indent, tidy
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); } - - -