# HG changeset patch # User Chris Cannam # Date 1559300072 -3600 # Node ID fdaa63607c15e441e8f351c2b38838c9bf289b28 # Parent cbe668c7d7249604456b1114b35a5069a967bd8e Untabify, indent, tidy diff -r cbe668c7d724 -r fdaa63607c15 base/KaiserWindow.cpp --- 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]); diff -r cbe668c7d724 -r fdaa63607c15 base/KaiserWindow.h --- 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: diff -r cbe668c7d724 -r fdaa63607c15 base/Pitch.cpp --- 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; diff -r cbe668c7d724 -r fdaa63607c15 base/Pitch.h --- 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); }; diff -r cbe668c7d724 -r fdaa63607c15 base/SincWindow.cpp --- 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); + } } } diff -r cbe668c7d724 -r fdaa63607c15 base/SincWindow.h --- 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: diff -r cbe668c7d724 -r fdaa63607c15 base/Window.h --- 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; } diff -r cbe668c7d724 -r fdaa63607c15 dsp/chromagram/Chromagram.cpp --- 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); diff -r cbe668c7d724 -r fdaa63607c15 dsp/chromagram/Chromagram.h --- 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 *m_window; double *m_windowbuf; - + double* m_chromadata; double m_FMin; double m_FMax; diff -r cbe668c7d724 -r fdaa63607c15 dsp/chromagram/ConstantQ.cpp --- 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; iis.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 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 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 diff -r cbe668c7d724 -r fdaa63607c15 dsp/signalconditioning/DFProcess.cpp --- 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(Config.LPACoeffs, Config.LPACoeffs + Config.LPOrd + 1); - params.b = std::vector(Config.LPBCoeffs, Config.LPBCoeffs + Config.LPOrd + 1); + params.a = std::vector + (Config.LPACoeffs, Config.LPACoeffs + Config.LPOrd + 1); + params.b = std::vector + (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; } } diff -r cbe668c7d724 -r fdaa63607c15 dsp/signalconditioning/DFProcess.h --- 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(); diff -r cbe668c7d724 -r fdaa63607c15 dsp/signalconditioning/FiltFilt.cpp --- 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() diff -r cbe668c7d724 -r fdaa63607c15 dsp/signalconditioning/Filter.cpp --- 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; } } diff -r cbe668c7d724 -r fdaa63607c15 dsp/signalconditioning/Framer.cpp --- 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++; diff -r cbe668c7d724 -r fdaa63607c15 dsp/signalconditioning/Framer.h --- 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 -#include -#include - - 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 diff -r cbe668c7d724 -r fdaa63607c15 dsp/transforms/DCT.cpp --- 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); diff -r cbe668c7d724 -r fdaa63607c15 dsp/wavelet/Wavelet.cpp --- 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); diff -r cbe668c7d724 -r fdaa63607c15 dsp/wavelet/Wavelet.h --- 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 }; diff -r cbe668c7d724 -r fdaa63607c15 hmm/hmm.c --- 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 #include #include -#include /* to seed random number generator */ +#include /* to seed random number generator */ -#include /* LAPACK for matrix inversion */ +#include /* LAPACK for matrix inversion */ #include "maths/nan-inf.h" @@ -37,788 +37,652 @@ #ifdef _MAC_OS_X #include #else -#include /* BLAS for matrix multiplication */ +#include /* 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"); } diff -r cbe668c7d724 -r fdaa63607c15 hmm/hmm.h --- 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 diff -r cbe668c7d724 -r fdaa63607c15 maths/pca/pca.c --- 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); } - - -