From 22b07e0233a29d9633ffa825a79503befaf2e16e Mon Sep 17 00:00:00 2001 From: Robin Gareus Date: Mon, 5 Oct 2015 16:17:49 +0200 Subject: NOOP, remove trailing tabs/whitespace. --- libs/qm-dsp/base/Window.h | 16 +- libs/qm-dsp/dsp/chromagram/Chromagram.cpp | 10 +- libs/qm-dsp/dsp/chromagram/Chromagram.h | 8 +- libs/qm-dsp/dsp/chromagram/ConstantQ.cpp | 12 +- libs/qm-dsp/dsp/chromagram/ConstantQ.h | 6 +- libs/qm-dsp/dsp/onsets/DetectionFunction.cpp | 22 +- libs/qm-dsp/dsp/onsets/DetectionFunction.h | 2 +- libs/qm-dsp/dsp/onsets/PeakPicking.cpp | 18 +- libs/qm-dsp/dsp/onsets/PeakPicking.h | 6 +- libs/qm-dsp/dsp/phasevocoder/PhaseVocoder.cpp | 6 +- libs/qm-dsp/dsp/rateconversion/Decimator.h | 2 +- .../dsp/segmentation/ClusterMeltSegmenter.cpp | 26 +-- .../qm-dsp/dsp/segmentation/ClusterMeltSegmenter.h | 6 +- libs/qm-dsp/dsp/segmentation/Segmenter.cpp | 4 +- libs/qm-dsp/dsp/segmentation/Segmenter.h | 2 +- libs/qm-dsp/dsp/segmentation/cluster_melt.c | 68 +++--- libs/qm-dsp/dsp/segmentation/cluster_segmenter.c | 76 +++---- libs/qm-dsp/dsp/signalconditioning/DFProcess.cpp | 16 +- libs/qm-dsp/dsp/signalconditioning/DFProcess.h | 2 +- libs/qm-dsp/dsp/signalconditioning/FiltFilt.cpp | 16 +- libs/qm-dsp/dsp/signalconditioning/Filter.h | 2 +- libs/qm-dsp/dsp/signalconditioning/Framer.cpp | 8 +- libs/qm-dsp/dsp/tempotracking/DownBeat.cpp | 2 +- libs/qm-dsp/dsp/tempotracking/TempoTrackV2.cpp | 8 +- libs/qm-dsp/dsp/tonal/ChangeDetectionFunction.cpp | 34 +-- libs/qm-dsp/dsp/tonal/ChangeDetectionFunction.h | 2 +- libs/qm-dsp/dsp/tonal/TCSgram.cpp | 4 +- libs/qm-dsp/dsp/tonal/TCSgram.h | 2 +- libs/qm-dsp/dsp/tonal/TonalEstimator.cpp | 20 +- libs/qm-dsp/dsp/tonal/TonalEstimator.h | 22 +- libs/qm-dsp/dsp/transforms/FFT.cpp | 2 +- libs/qm-dsp/dsp/wavelet/Wavelet.cpp | 8 +- libs/qm-dsp/hmm/hmm.c | 242 ++++++++++----------- libs/qm-dsp/maths/MathUtilities.cpp | 26 +-- libs/qm-dsp/maths/MathUtilities.h | 2 +- libs/qm-dsp/maths/pca/pca.c | 20 +- 36 files changed, 364 insertions(+), 364 deletions(-) (limited to 'libs/qm-dsp') diff --git a/libs/qm-dsp/base/Window.h b/libs/qm-dsp/base/Window.h index 686d023c6b..02ca426d66 100644 --- a/libs/qm-dsp/base/Window.h +++ b/libs/qm-dsp/base/Window.h @@ -72,53 +72,53 @@ void Window::encache() 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; - + case BartlettWindow: for (i = 0; i < n/2; ++i) { mult[i] = mult[i] * (i / T(n/2)); mult[i + n/2] = mult[i + n/2] * (1.0 - (i / T(n/2))); } break; - + case HammingWindow: for (i = 0; i < n; ++i) { mult[i] = mult[i] * (0.54 - 0.46 * cos(2 * M_PI * i / n)); } break; - + case HanningWindow: for (i = 0; i < n; ++i) { mult[i] = mult[i] * (0.50 - 0.50 * cos(2 * M_PI * i / n)); } break; - + case BlackmanWindow: 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; - + case GaussianWindow: for (i = 0; i < n; ++i) { mult[i] = mult[i] * exp((-1.0 / (n*n)) * ((T(2*i) - n) * (T(2*i) - n))); } break; - + case ParzenWindow: for (i = 0; i < n; ++i) { mult[i] = mult[i] * (1.0 - fabs((T(2*i) - n) / T(n + 1))); } break; } - + m_cache = mult; } diff --git a/libs/qm-dsp/dsp/chromagram/Chromagram.cpp b/libs/qm-dsp/dsp/chromagram/Chromagram.cpp index 83a6661d17..5901082352 100644 --- a/libs/qm-dsp/dsp/chromagram/Chromagram.cpp +++ b/libs/qm-dsp/dsp/chromagram/Chromagram.cpp @@ -27,14 +27,14 @@ Chromagram::Chromagram( ChromaConfig Config ) : } 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_normalise = Config.normalise; // if frame normalisation is required // No. of constant Q bins - m_uK = ( unsigned int ) ceil( m_BPO * log(m_FMax/m_FMin)/log(2.0)); + m_uK = ( unsigned int ) ceil( m_BPO * log(m_FMax/m_FMin)/log(2.0)); // Create array for chroma result m_chromadata = new double[ m_BPO ]; @@ -49,7 +49,7 @@ int Chromagram::initialise( ChromaConfig Config ) ConstantQConfig.max = m_FMax; ConstantQConfig.BPO = m_BPO; ConstantQConfig.CQThresh = Config.CQThresh; - + // Initialise ConstantQ operator m_ConstantQ = new ConstantQ( ConstantQConfig ); @@ -57,7 +57,7 @@ int Chromagram::initialise( ChromaConfig Config ) 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 ]; @@ -161,7 +161,7 @@ double* Chromagram::process( const double *real, const double *imag ) // Calculate ConstantQ frame m_ConstantQ->process( real, imag, m_CQRe, m_CQIm ); - + // add each octave of cq data into Chromagram const unsigned octaves = (int)floor(double( m_uK/m_BPO))-1; for (unsigned octave = 0; octave <= octaves; octave++) diff --git a/libs/qm-dsp/dsp/chromagram/Chromagram.h b/libs/qm-dsp/dsp/chromagram/Chromagram.h index 37af153be3..f04ffda83a 100644 --- a/libs/qm-dsp/dsp/chromagram/Chromagram.h +++ b/libs/qm-dsp/dsp/chromagram/Chromagram.h @@ -32,17 +32,17 @@ struct ChromaConfig{ class Chromagram { -public: +public: Chromagram( ChromaConfig Config ); ~Chromagram(); - + double* process( const double *data ); // time domain double* process( const double *real, const double *imag ); // frequency domain void unityNormalise( double* src ); // Complex arithmetic double kabs( double real, double imag ); - + // Results unsigned int getK() { return m_uK;} unsigned int getFrameSize() { return m_frameSize; } @@ -54,7 +54,7 @@ private: Window *m_window; double *m_windowbuf; - + double* m_chromadata; double m_FMin; double m_FMax; diff --git a/libs/qm-dsp/dsp/chromagram/ConstantQ.cpp b/libs/qm-dsp/dsp/chromagram/ConstantQ.cpp index 222fd80a36..fa6c32c26b 100644 --- a/libs/qm-dsp/dsp/chromagram/ConstantQ.cpp +++ b/libs/qm-dsp/dsp/chromagram/ConstantQ.cpp @@ -109,14 +109,14 @@ void ConstantQ::sparsekernel() 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++) @@ -152,13 +152,13 @@ void ConstantQ::sparsekernel() //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, doubled because they are floats sk->is.push_back(j); sk->js.push_back(k); @@ -266,7 +266,7 @@ double* ConstantQ::process( const double* fftdata ) const double *real = &(sk->real[0]); const double *imag = &(sk->imag[0]); const unsigned int sparseCells = sk->real.size(); - + for (unsigned i = 0; ireal[0]); const double *imag = &(sk->imag[0]); const unsigned int sparseCells = sk->real.size(); - + for (unsigned i = 0; i &onsets ) { if (len < 4) return; - vector m_maxima; + vector 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 ] ); + m_maxima.push_back( m_workBuffer[ u ] ); } - + quadEval( m_maxima, onsets ); for(unsigned int b = 0; b < m_maxima.size(); b++) @@ -92,7 +92,7 @@ int PeakPicking::quadEval( vector &src, vector &idx ) vector m_maxIndex; vector m_onsetPosition; - + vector m_maxFit; vector m_poly; vector m_err; @@ -123,7 +123,7 @@ int PeakPicking::quadEval( vector &src, vector &idx ) for (int k = -2; k <= 2; ++k) { selMax = src[ m_maxIndex[j] + k ] ; - m_maxFit.push_back(selMax); + m_maxFit.push_back(selMax); } double f = m_poly[0]; @@ -133,7 +133,7 @@ int PeakPicking::quadEval( vector &src, vector &idx ) { idx.push_back(m_maxIndex[j]); } - + m_maxFit.clear(); } diff --git a/libs/qm-dsp/dsp/onsets/PeakPicking.h b/libs/qm-dsp/dsp/onsets/PeakPicking.h index be6853ca8f..931c1c16e2 100644 --- a/libs/qm-dsp/dsp/onsets/PeakPicking.h +++ b/libs/qm-dsp/dsp/onsets/PeakPicking.h @@ -56,7 +56,7 @@ class PeakPicking public: PeakPicking( PPickParams Config ); virtual ~PeakPicking(); - + void process( double* src, unsigned int len, vector &onsets ); @@ -64,7 +64,7 @@ private: void initialise( PPickParams Config ); void deInitialise(); int quadEval( vector &src, vector &idx ); - + DFProcConfig m_DFProcessingParams; unsigned int m_DFLength ; @@ -74,7 +74,7 @@ private: double* m_workBuffer; - + DFProcess* m_DFSmoothing; }; diff --git a/libs/qm-dsp/dsp/phasevocoder/PhaseVocoder.cpp b/libs/qm-dsp/dsp/phasevocoder/PhaseVocoder.cpp index 26cfcfbabe..cd52f5240a 100644 --- a/libs/qm-dsp/dsp/phasevocoder/PhaseVocoder.cpp +++ b/libs/qm-dsp/dsp/phasevocoder/PhaseVocoder.cpp @@ -49,7 +49,7 @@ void PhaseVocoder::FFTShift(unsigned int size, double *src) void PhaseVocoder::process(double *src, double *mag, double *theta) { FFTShift( m_n, src); - + m_fft->process(0, src, m_realOut, m_imagOut); getMagnitude( m_n/2, mag, m_realOut, m_imagOut); @@ -57,7 +57,7 @@ void PhaseVocoder::process(double *src, double *mag, double *theta) } void PhaseVocoder::getMagnitude(unsigned int size, double *mag, double *real, double *imag) -{ +{ unsigned int j; for( j = 0; j < size; j++) @@ -75,5 +75,5 @@ void PhaseVocoder::getPhase(unsigned int size, double *theta, double *real, doub for( k = 0; k < size; k++) { theta[ k ] = atan2( -imag[ k ], real[ k ]); - } + } } diff --git a/libs/qm-dsp/dsp/rateconversion/Decimator.h b/libs/qm-dsp/dsp/rateconversion/Decimator.h index 1ed49fd397..f8a113a0db 100644 --- a/libs/qm-dsp/dsp/rateconversion/Decimator.h +++ b/libs/qm-dsp/dsp/rateconversion/Decimator.h @@ -55,7 +55,7 @@ private: double a[ 9 ]; double b[ 9 ]; - + double* decBuffer; }; diff --git a/libs/qm-dsp/dsp/segmentation/ClusterMeltSegmenter.cpp b/libs/qm-dsp/dsp/segmentation/ClusterMeltSegmenter.cpp index 7643691cd3..ce5f370436 100644 --- a/libs/qm-dsp/dsp/segmentation/ClusterMeltSegmenter.cpp +++ b/libs/qm-dsp/dsp/segmentation/ClusterMeltSegmenter.cpp @@ -212,7 +212,7 @@ void ClusterMeltSegmenter::extractFeaturesConstQ(const double* samples, int nsam fft->process(false, frame, real, imag); constq->process(real, imag, cqre, cqim); - + for (int i = 0; i < ncoeff; ++i) { cq[i] += sqrt(cqre[i] * cqre[i] + cqim[i] * cqim[i]); } @@ -287,7 +287,7 @@ void ClusterMeltSegmenter::extractFeaturesMFCC(const double* samples, int nsampl } mfcc->process(frame, ccout); - + for (int i = 0; i < ncoeff; ++i) { cc[i] += ccout[i]; } @@ -335,22 +335,22 @@ void ClusterMeltSegmenter::segment() << " features with " << features[0].size() << " coefficients (ncoeff = " << ncoeff << ", ncomponents = " << ncomponents << ")" << std::endl; */ // copy the features to a native array and use the existing C segmenter... - double** arrFeatures = new double*[features.size()]; + double** arrFeatures = new double*[features.size()]; for (int i = 0; i < features.size(); i++) { if (featureType == FEATURE_TYPE_UNKNOWN) { arrFeatures[i] = new double[features[0].size()]; for (int j = 0; j < features[0].size(); j++) - arrFeatures[i][j] = features[i][j]; + arrFeatures[i][j] = features[i][j]; } else { arrFeatures[i] = new double[ncoeff+1]; // allow space for the normalised envelope for (int j = 0; j < ncoeff; j++) - arrFeatures[i][j] = features[i][j]; + arrFeatures[i][j] = features[i][j]; } } - + q = new int[features.size()]; - + if (featureType == FEATURE_TYPE_UNKNOWN || featureType == FEATURE_TYPE_MFCC) cluster_segment(q, arrFeatures, features.size(), features[0].size(), nHMMStates, histogramLength, @@ -358,16 +358,16 @@ void ClusterMeltSegmenter::segment() else constq_segment(q, arrFeatures, features.size(), nbins, ncoeff, featureType, nHMMStates, histogramLength, nclusters, neighbourhoodLimit); - + // convert the cluster assignment sequence to a segmentation - makeSegmentation(q, features.size()); - + makeSegmentation(q, features.size()); + // de-allocate arrays delete [] q; for (int i = 0; i < features.size(); i++) delete [] arrFeatures[i]; delete [] arrFeatures; - + // clear the features clear(); } @@ -377,11 +377,11 @@ void ClusterMeltSegmenter::makeSegmentation(int* q, int len) segmentation.segments.clear(); segmentation.nsegtypes = nclusters; segmentation.samplerate = samplerate; - + Segment segment; segment.start = 0; segment.type = q[0]; - + for (int i = 1; i < len; i++) { if (q[i] != q[i-1]) diff --git a/libs/qm-dsp/dsp/segmentation/ClusterMeltSegmenter.h b/libs/qm-dsp/dsp/segmentation/ClusterMeltSegmenter.h index 528c09cb55..8f3130871e 100644 --- a/libs/qm-dsp/dsp/segmentation/ClusterMeltSegmenter.h +++ b/libs/qm-dsp/dsp/segmentation/ClusterMeltSegmenter.h @@ -72,7 +72,7 @@ public: protected: void makeSegmentation(int* q, int len); - + void extractFeaturesConstQ(const double *, int); void extractFeaturesMFCC(const double *, int); @@ -82,9 +82,9 @@ protected: MFCC* mfcc; model_t* model; // the HMM int* q; // the decoded HMM state sequence - vector > histograms; + vector > histograms; - feature_types featureType; + feature_types featureType; double hopSize; // in seconds double windowSize; // in seconds diff --git a/libs/qm-dsp/dsp/segmentation/Segmenter.cpp b/libs/qm-dsp/dsp/segmentation/Segmenter.cpp index 120a6617f5..b60fb58162 100644 --- a/libs/qm-dsp/dsp/segmentation/Segmenter.cpp +++ b/libs/qm-dsp/dsp/segmentation/Segmenter.cpp @@ -19,13 +19,13 @@ ostream& operator<<(ostream& os, const Segmentation& s) { os << "structure_name : begin_time end_time\n"; - + for (int i = 0; i < s.segments.size(); i++) { Segment seg = s.segments[i]; os << std::fixed << seg.type << ':' << '\t' << std::setprecision(6) << seg.start / static_cast(s.samplerate) << '\t' << std::setprecision(6) << seg.end / static_cast(s.samplerate) << "\n"; } - + return os; } diff --git a/libs/qm-dsp/dsp/segmentation/Segmenter.h b/libs/qm-dsp/dsp/segmentation/Segmenter.h index 9a77d70372..fd2f39b850 100644 --- a/libs/qm-dsp/dsp/segmentation/Segmenter.h +++ b/libs/qm-dsp/dsp/segmentation/Segmenter.h @@ -35,7 +35,7 @@ class Segmentation public: int nsegtypes; // number of segment types, so possible types are {0,1,...,nsegtypes-1} int samplerate; - vector segments; + vector segments; }; ostream& operator<<(ostream& os, const Segmentation& s); diff --git a/libs/qm-dsp/dsp/segmentation/cluster_melt.c b/libs/qm-dsp/dsp/segmentation/cluster_melt.c index 0509480807..092bc7f078 100644 --- a/libs/qm-dsp/dsp/segmentation/cluster_melt.c +++ b/libs/qm-dsp/dsp/segmentation/cluster_melt.c @@ -25,7 +25,7 @@ double kldist(double* a, double* b, int n) { because a, b represent probability distributions */ double q, d; int i; - + d = 0; for (i = 0; i < n; i++) { @@ -38,8 +38,8 @@ double kldist(double* a, double* b, int n) { d += b[i] * log(b[i] / q); } } - return d; -} + return d; +} void cluster_melt(double *h, int m, int n, double *Bsched, int t, int k, int l, int *c) { double lambda, sum, beta, logsumexp, maxlp; @@ -48,9 +48,9 @@ void cluster_melt(double *h, int m, int n, double *Bsched, int t, int k, int l, int** nc; /* neighbour counts for each histogram */ double** lp; /* soft assignment probs for each histogram */ int* oldc; /* previous hard assignments (to check convergence) */ - + /* NB h is passed as a 1d row major array */ - + /* parameter values */ lambda = DEFAULT_LAMBDA; if (l > 0) @@ -60,22 +60,22 @@ void cluster_melt(double *h, int m, int n, double *Bsched, int t, int k, int l, B = 2 * limit + 1; maxiter0 = 20; /* number of iterations at initial temperature */ maxiter1 = 5; /* number of iterations at subsequent temperatures */ - - /* allocate memory */ + + /* allocate memory */ cl = (double**) malloc(k*sizeof(double*)); for (i= 0; i < k; i++) cl[i] = (double*) malloc(m*sizeof(double)); - + nc = (int**) malloc(n*sizeof(int*)); for (i= 0; i < n; i++) nc[i] = (int*) malloc(k*sizeof(int)); - + lp = (double**) malloc(n*sizeof(double*)); for (i= 0; i < n; i++) lp[i] = (double*) malloc(k*sizeof(double)); - + oldc = (int*) malloc(n * sizeof(int)); - + /* initialise */ for (i = 0; i < k; i++) { @@ -90,40 +90,40 @@ void cluster_melt(double *h, int m, int n, double *Bsched, int t, int k, int l, { cl[i][j] /= sum; /* normalise */ } - } + } //print_array(cl, k, m); - + for (i = 0; i < n; i++) c[i] = 1; /* initially assign all histograms to cluster 1 */ - + for (a = 0; a < t; a++) { beta = Bsched[a]; - + if (a == 0) maxiter = maxiter0; else maxiter = maxiter1; - + for (it = 0; it < maxiter; it++) { //if (it == maxiter - 1) // mexPrintf("hasn't converged after %d iterations\n", maxiter); - + for (i = 0; i < n; i++) { /* save current hard assignments */ oldc[i] = c[i]; - + /* calculate soft assignment logprobs for each cluster */ sum = 0; for (j = 0; j < k; j++) { lp[i][ j] = -beta * kldist(cl[j], &h[i*m], m); - + /* update matching neighbour counts for this histogram, based on current hard assignments */ /* old version: - nc[i][j] = 0; + nc[i][j] = 0; if (i >= limit && i <= n - 1 - limit) { for (b = i - limit; b <= i + limit; b++) @@ -144,14 +144,14 @@ void cluster_melt(double *h, int m, int n, double *Bsched, int t, int k, int l, for (b = b0; b <= b1; b++) if (c[b] == j+1) nc[i][j]--; - + sum += exp(lp[i][j]); } - + /* normalise responsibilities and add duration logprior */ logsumexp = log(sum); for (j = 0; j < k; j++) - lp[i][j] -= logsumexp + lambda * nc[i][j]; + lp[i][j] -= logsumexp + lambda * nc[i][j]; } //print_array(lp, n, k); /* @@ -162,8 +162,8 @@ void cluster_melt(double *h, int m, int n, double *Bsched, int t, int k, int l, mexPrintf("\n"); } */ - - + + /* update the assignments now that we know the duration priors based on the current assignments */ for (i = 0; i < n; i++) @@ -177,14 +177,14 @@ void cluster_melt(double *h, int m, int n, double *Bsched, int t, int k, int l, c[i] = j+1; } } - + /* break if assignments haven't changed */ i = 0; while (i < n && oldc[i] == c[i]) i++; if (i == n) break; - + /* update reference histograms now we know new responsibilities */ for (j = 0; j < k; j++) { @@ -194,21 +194,21 @@ void cluster_melt(double *h, int m, int n, double *Bsched, int t, int k, int l, for (i = 0; i < n; i++) { cl[j][b] += exp(lp[i][j]) * h[i*m+b]; - } + } } - - sum = 0; + + sum = 0; for (i = 0; i < n; i++) sum += exp(lp[i][j]); for (b = 0; b < m; b++) cl[j][b] /= sum; /* normalise */ - } - + } + //print_array(cl, k, m); //mexPrintf("\n\n"); } } - + /* free memory */ for (i = 0; i < k; i++) free(cl[i]); @@ -219,7 +219,7 @@ void cluster_melt(double *h, int m, int n, double *Bsched, int t, int k, int l, for (i = 0; i < n; i++) free(lp[i]); free(lp); - free(oldc); + free(oldc); } diff --git a/libs/qm-dsp/dsp/segmentation/cluster_segmenter.c b/libs/qm-dsp/dsp/segmentation/cluster_segmenter.c index 0d2762ee7f..c9f115c205 100644 --- a/libs/qm-dsp/dsp/segmentation/cluster_segmenter.c +++ b/libs/qm-dsp/dsp/segmentation/cluster_segmenter.c @@ -25,7 +25,7 @@ void cq2chroma(double** cq, int nframes, int ncoeff, int bins, double** chroma) int t, b, oct, ix; //double maxchroma; /* max chroma value at each time, for normalisation */ //double sum; /* for normalisation */ - + for (t = 0; t < nframes; t++) { for (b = 0; b < bins; b++) @@ -50,7 +50,7 @@ void cq2chroma(double** cq, int nframes, int ncoeff, int bins, double** chroma) maxchroma = chroma[t][b]; if (maxchroma > 0) for (b = 0; b < bins; b++) - chroma[t][b] /= maxchroma; + chroma[t][b] /= maxchroma; */ } } @@ -62,13 +62,13 @@ void mpeg7_constq(double** features, int nframes, int ncoeff) double ss; double env; double maxenv = 0; - + /* convert const-Q features to dB scale */ for (i = 0; i < nframes; i++) for (j = 0; j < ncoeff; j++) features[i][j] = 10.0 * log10(features[i][j]+DBL_EPSILON); - - /* normalise each feature vector and add the norm as an extra feature dimension */ + + /* normalise each feature vector and add the norm as an extra feature dimension */ for (i = 0; i < nframes; i++) { ss = 0; @@ -83,7 +83,7 @@ void mpeg7_constq(double** features, int nframes, int ncoeff) } /* normalise the envelopes */ for (i = 0; i < nframes; i++) - features[i][ncoeff] /= maxenv; + features[i][ncoeff] /= maxenv; } /* return histograms h[nx*m] of data x[nx] into m bins using a sliding window of length h_len (MUST BE ODD) */ @@ -109,7 +109,7 @@ void create_histograms(int* x, int nx, int m, int hlen, double* h) for (j = 0; j < m; j++) h[i*m+j] /= norm; } - + /* duplicate histograms at beginning and end to create one histogram for each data value supplied */ for (i = 0; i < hlen/2; i++) for (j = 0; j < m; j++) @@ -124,7 +124,7 @@ void cluster_segment(int* q, double** features, int frames_read, int feature_len int histogram_length, int nclusters, int neighbour_limit) { int i, j; - + /*****************************/ if (0) { /* try just using the predominant bin number as a 'decoded state' */ @@ -141,44 +141,44 @@ void cluster_segment(int* q, double** features, int frames_read, int feature_len { maxval = features[i][j]; maxbin = j; - } + } } if (maxval > chroma_thresh) q[i] = maxbin; else q[i] = feature_length; } - + } if (1) { /*****************************/ - - + + /* scale all the features to 'balance covariances' during HMM training */ double scale = 10; for (i = 0; i < frames_read; i++) for (j = 0; j < feature_length; j++) features[i][j] *= scale; - + /* train an HMM on the features */ - + /* create a model */ model_t* model = hmm_init(features, frames_read, feature_length, nHMM_states); - + /* train the model */ hmm_train(features, frames_read, model); -/* +/* printf("\n\nafter training:\n"); hmm_print(model); -*/ +*/ /* decode the hidden state sequence */ viterbi_decode(features, frames_read, model, q); hmm_close(model); - + /*****************************/ } /*****************************/ - + /* fprintf(stderr, "HMM state sequence:\n"); @@ -186,11 +186,11 @@ void cluster_segment(int* q, double** features, int frames_read, int feature_len fprintf(stderr, "%d ", q[i]); fprintf(stderr, "\n\n"); */ - + /* create histograms of states */ double* h = (double*) malloc(frames_read*nHMM_states*sizeof(double)); /* vector in row major order */ create_histograms(q, frames_read, nHMM_states, histogram_length, h); - + /* cluster the histograms */ int nbsched = 20; /* length of inverse temperature schedule */ double* bsched = (double*) malloc(nbsched*sizeof(double)); /* inverse temperature schedule */ @@ -200,9 +200,9 @@ void cluster_segment(int* q, double** features, int frames_read, int feature_len for (i = 1; i < nbsched; i++) bsched[i] = alpha * bsched[i-1]; cluster_melt(h, nHMM_states, frames_read, bsched, nbsched, nclusters, neighbour_limit, q); - + /* now q holds a sequence of cluster assignments */ - + free(h); free(bsched); } @@ -214,25 +214,25 @@ void constq_segment(int* q, double** features, int frames_read, int bins, int nc int feature_length; double** chroma; int i; - + if (feature_type == FEATURE_TYPE_CONSTQ) { /* fprintf(stderr, "Converting to dB and normalising...\n"); - */ + */ mpeg7_constq(features, frames_read, ncoeff); -/* +/* fprintf(stderr, "Running PCA...\n"); -*/ +*/ /* do PCA on the features (but not the envelope) */ int ncomponents = 20; pca_project(features, frames_read, ncoeff, ncomponents); - + /* copy the envelope so that it immediatly follows the chosen components */ for (i = 0; i < frames_read; i++) - features[i][ncomponents] = features[i][ncoeff]; - + features[i][ncomponents] = features[i][ncoeff]; + feature_length = ncomponents + 1; - + /************************************** //TEST // feature file name @@ -241,7 +241,7 @@ void constq_segment(int* q, double** features, int frames_read, int bins, int nc strcpy(file_name, dir); strcat(file_name, trackname); strcat(file_name, "_features_c20r8h0.2f0.6.mat"); - + // get the features from Matlab from mat-file int frames_in_file; readmatarray_size(file_name, 2, &frames_in_file, &feature_length); @@ -254,27 +254,27 @@ void constq_segment(int* q, double** features, int frames_read, int bins, int nc features[frames_read-missing_frames][i] = features[frames_read-missing_frames-1][i]; --missing_frames; } - + free(file_name); ******************************************/ - + cluster_segment(q, features, frames_read, feature_length, nHMM_states, histogram_length, nclusters, neighbour_limit); } - + if (feature_type == FEATURE_TYPE_CHROMA) { /* fprintf(stderr, "Converting to chroma features...\n"); -*/ +*/ /* convert constant-Q to normalised chroma features */ chroma = (double**) malloc(frames_read*sizeof(double*)); for (i = 0; i < frames_read; i++) chroma[i] = (double*) malloc(bins*sizeof(double)); cq2chroma(features, frames_read, ncoeff, bins, chroma); feature_length = bins; - + cluster_segment(q, chroma, frames_read, feature_length, nHMM_states, histogram_length, nclusters, neighbour_limit); - + for (i = 0; i < frames_read; i++) free(chroma[i]); free(chroma); diff --git a/libs/qm-dsp/dsp/signalconditioning/DFProcess.cpp b/libs/qm-dsp/dsp/signalconditioning/DFProcess.cpp index 971f5fac8e..4b93cb0e12 100644 --- a/libs/qm-dsp/dsp/signalconditioning/DFProcess.cpp +++ b/libs/qm-dsp/dsp/signalconditioning/DFProcess.cpp @@ -25,7 +25,7 @@ DFProcess::DFProcess( DFProcConfig Config ) { filtSrc = NULL; - filtDst = NULL; + filtDst = NULL; m_filtScratchIn = NULL; m_filtScratchOut = NULL; @@ -51,13 +51,13 @@ void DFProcess::initialise( DFProcConfig Config ) filtSrc = new double[ m_length ]; filtDst = new double[ m_length ]; - + //Low Pass Smoothing Filter Config m_FilterConfigParams.ord = Config.LPOrd; m_FilterConfigParams.ACoeffs = Config.LPACoeffs; m_FilterConfigParams.BCoeffs = Config.LPBCoeffs; - - m_FiltFilt = new FiltFilt( m_FilterConfigParams ); + + m_FiltFilt = new FiltFilt( m_FilterConfigParams ); } void DFProcess::deInitialise() @@ -115,7 +115,7 @@ void DFProcess::medianFilter(double *src, double *dst) { if (index >= m_length) break; - + l = 0; for( j = i; j < ( i + m_winPost + m_winPre + 1); j++) { @@ -139,7 +139,7 @@ void DFProcess::medianFilter(double *src, double *dst) l++; } - + scratch[ index++ ] = MathUtilities::median( y, l); } @@ -147,7 +147,7 @@ void DFProcess::medianFilter(double *src, double *dst) for( i = 0; i < m_length; i++ ) { val = src[ i ] - scratch[ i ];// - 0.033; - + if( m_isMedianPositive ) { if( val > 0 ) @@ -164,7 +164,7 @@ void DFProcess::medianFilter(double *src, double *dst) dst[ i ] = val; } } - + delete [] y; delete [] scratch; } diff --git a/libs/qm-dsp/dsp/signalconditioning/DFProcess.h b/libs/qm-dsp/dsp/signalconditioning/DFProcess.h index 3af80743f5..ad40f846b9 100644 --- a/libs/qm-dsp/dsp/signalconditioning/DFProcess.h +++ b/libs/qm-dsp/dsp/signalconditioning/DFProcess.h @@ -38,7 +38,7 @@ public: void process( double* src, double* dst ); - + private: void initialise( DFProcConfig Config ); void deInitialise(); diff --git a/libs/qm-dsp/dsp/signalconditioning/FiltFilt.cpp b/libs/qm-dsp/dsp/signalconditioning/FiltFilt.cpp index 5a9e2c9549..13dbda0e82 100644 --- a/libs/qm-dsp/dsp/signalconditioning/FiltFilt.cpp +++ b/libs/qm-dsp/dsp/signalconditioning/FiltFilt.cpp @@ -24,7 +24,7 @@ FiltFilt::FiltFilt( FiltFiltConfig Config ) m_filtScratchIn = NULL; m_filtScratchOut = NULL; m_ord = 0; - + initialise( Config ); } @@ -39,7 +39,7 @@ void FiltFilt::initialise( FiltFiltConfig Config ) m_filterConfig.ord = Config.ord; m_filterConfig.ACoeffs = Config.ACoeffs; m_filterConfig.BCoeffs = Config.BCoeffs; - + m_filter = new Filter( m_filterConfig ); } @@ -50,7 +50,7 @@ void FiltFilt::deInitialise() void FiltFilt::process(double *src, double *dst, unsigned int length) -{ +{ unsigned int i; if (length == 0) return; @@ -62,7 +62,7 @@ void FiltFilt::process(double *src, double *dst, unsigned int length) m_filtScratchIn = new double[ nExt ]; m_filtScratchOut = new double[ nExt ]; - + for( i = 0; i< nExt; i++ ) { m_filtScratchIn[ i ] = 0.0; @@ -89,11 +89,11 @@ void FiltFilt::process(double *src, double *dst, unsigned int length) { m_filtScratchIn[ i + nFact ] = src[ i ]; } - + //////////////////////////////// // Do 0Ph filtering m_filter->process( m_filtScratchIn, m_filtScratchOut, nExt); - + // reverse the series for FILTFILT for ( i = 0; i < nExt; i++) { @@ -102,7 +102,7 @@ void FiltFilt::process(double *src, double *dst, unsigned int length) // do FILTER again m_filter->process( m_filtScratchIn, m_filtScratchOut, nExt); - + // reverse the series back for ( i = 0; i < nExt; i++) { @@ -117,7 +117,7 @@ void FiltFilt::process(double *src, double *dst, unsigned int length) for( i = 0; i < length; i++ ) { dst[ index++ ] = m_filtScratchOut[ i + nFact ]; - } + } delete [] m_filtScratchIn; delete [] m_filtScratchOut; diff --git a/libs/qm-dsp/dsp/signalconditioning/Filter.h b/libs/qm-dsp/dsp/signalconditioning/Filter.h index 8533ed0502..b1c20d506a 100644 --- a/libs/qm-dsp/dsp/signalconditioning/Filter.h +++ b/libs/qm-dsp/dsp/signalconditioning/Filter.h @@ -35,7 +35,7 @@ public: void reset(); void process( double *src, double *dst, unsigned int length ); - + private: void initialise( FilterConfig Config ); diff --git a/libs/qm-dsp/dsp/signalconditioning/Framer.cpp b/libs/qm-dsp/dsp/signalconditioning/Framer.cpp index 80ea67f72c..7357d2417a 100644 --- a/libs/qm-dsp/dsp/signalconditioning/Framer.cpp +++ b/libs/qm-dsp/dsp/signalconditioning/Framer.cpp @@ -44,14 +44,14 @@ void Framer::configure( unsigned int frameLength, unsigned int hop ) if( m_dataFrame != NULL ) { - delete [] m_dataFrame; + delete [] m_dataFrame; m_dataFrame = NULL; } m_dataFrame = new double[ m_frameLength ]; if( m_strideFrame != NULL ) { - delete [] m_strideFrame; + delete [] m_strideFrame; m_strideFrame = NULL; } m_strideFrame = new double[ m_stepSize ]; @@ -65,7 +65,7 @@ void Framer::getFrame(double *dst) for( unsigned int u = 0; u < m_frameLength; u++) { dst[ u ] = m_srcBuffer[ m_ulSrcIndex++ ]; - } + } m_ulSrcIndex -= ( m_frameLength - m_stepSize ); } else @@ -77,7 +77,7 @@ void Framer::getFrame(double *dst) { dst[ u ] = m_srcBuffer[ m_ulSrcIndex++ ]; } - + for( unsigned int u = 0; u < zero; u++ ) { dst[ rem + u ] = 0; diff --git a/libs/qm-dsp/dsp/tempotracking/DownBeat.cpp b/libs/qm-dsp/dsp/tempotracking/DownBeat.cpp index b2183b51ed..1167bc6ea9 100644 --- a/libs/qm-dsp/dsp/tempotracking/DownBeat.cpp +++ b/libs/qm-dsp/dsp/tempotracking/DownBeat.cpp @@ -293,7 +293,7 @@ DownBeat::measureSpecDiff(d_vec_t oldspec, d_vec_t newspec) } // JENSEN-SHANNON CALCULATION - sd1 = 0.5*oldspec[i] + 0.5*newspec[i]; + sd1 = 0.5*oldspec[i] + 0.5*newspec[i]; SD = SD + (-sd1*log(sd1)) + (0.5*(oldspec[i]*log(oldspec[i]))) + (0.5*(newspec[i]*log(newspec[i]))); } diff --git a/libs/qm-dsp/dsp/tempotracking/TempoTrackV2.cpp b/libs/qm-dsp/dsp/tempotracking/TempoTrackV2.cpp index acaa8ff7cd..a2c2a24b0a 100644 --- a/libs/qm-dsp/dsp/tempotracking/TempoTrackV2.cpp +++ b/libs/qm-dsp/dsp/tempotracking/TempoTrackV2.cpp @@ -177,7 +177,7 @@ TempoTrackV2::get_rcf(const d_vec_t &dfframe_in, const d_vec_t &wv, d_vec_t &rcf // now apply comb filtering int numelem = 4; - + for (unsigned int i = 2;i < rcf.size();i++) // max beat period { for (int a = 1;a <= numelem;a++) // number of comb elements @@ -218,7 +218,7 @@ TempoTrackV2::viterbi_decode(const d_mat_t &rcfmat, const d_vec_t &wv, d_vec_t & { tmat.push_back ( d_vec_t() ); // adds a new column for (unsigned int j=0; j(i); tmat[i][j] = exp( (-1.*pow((j-mu),2.)) / (2.*pow(sigma,2.)) ); } @@ -246,7 +246,7 @@ TempoTrackV2::viterbi_decode(const d_mat_t &rcfmat, const d_vec_t &wv, d_vec_t & delta.push_back( d_vec_t()); psi.push_back( i_vec_t()); for (unsigned int j=0; j filter width = 2*FWHM = 2*2.3548*sigma m_dFilterSigma = double(m_iFilterWidth) / double(2*2.3548); m_vaGaussian.resize(m_iFilterWidth); - + double dScale = 1.0 / (m_dFilterSigma*sqrt(2*PI)); - + for (int x = -(m_iFilterWidth-1)/2; x <= (m_iFilterWidth-1)/2; x++) { double w = dScale * std::exp ( -(x*x)/(2*m_dFilterSigma*m_dFilterSigma) ); m_vaGaussian[x + (m_iFilterWidth-1)/2] = w; } - + #ifdef DEBUG_CHANGE_DETECTION_FUNCTION std::cerr << "Filter sigma: " << m_dFilterSigma << std::endl; std::cerr << "Filter width: " << m_iFilterWidth << std::endl; @@ -59,37 +59,37 @@ ChangeDistance ChangeDetectionFunction::process(const TCSGram& rTCSGram) { ChangeDistance retVal; retVal.resize(rTCSGram.getSize(), 0.0); - + TCSGram smoothedTCSGram; for (int iPosition = 0; iPosition < rTCSGram.getSize(); iPosition++) { int iSkipLower = 0; - + int iLowerPos = iPosition - (m_iFilterWidth-1)/2; int iUpperPos = iPosition + (m_iFilterWidth-1)/2; - + if (iLowerPos < 0) { iSkipLower = -iLowerPos; iLowerPos = 0; } - + if (iUpperPos >= rTCSGram.getSize()) { int iMaxIndex = rTCSGram.getSize() - 1; iUpperPos = iMaxIndex; } - + TCSVector smoothedVector; // for every bin of the vector, calculate the smoothed value for (int iPC = 0; iPC < 6; iPC++) - { + { size_t j = 0; double dSmoothedValue = 0.0; TCSVector rCV; - + for (int i = iLowerPos; i <= iUpperPos; i++) { rTCSGram.getTCSVector(i, rCV); @@ -98,7 +98,7 @@ ChangeDistance ChangeDetectionFunction::process(const TCSGram& rTCSGram) smoothedVector[iPC] = dSmoothedValue; } - + smoothedTCSGram.addTCSVector(smoothedVector); } @@ -109,10 +109,10 @@ ChangeDistance ChangeDetectionFunction::process(const TCSGram& rTCSGram) if the current estimate is not confident enough, look further into the future/the past e.g., High frequency content, zero crossing rate, spectral flatness */ - + TCSVector nextTCS; TCSVector previousTCS; - + int iWindow = 1; // while (previousTCS.magnitude() < 0.1 && (iPosition-iWindow) > 0) @@ -121,9 +121,9 @@ ChangeDistance ChangeDetectionFunction::process(const TCSGram& rTCSGram) // std::cout << previousTCS.magnitude() << std::endl; iWindow++; } - + iWindow = 1; - + // while (nextTCS.magnitude() < 0.1 && (iPosition+iWindow) < (rTCSGram.getSize()-1) ) { smoothedTCSGram.getTCSVector(iPosition+iWindow, nextTCS); @@ -136,7 +136,7 @@ ChangeDistance ChangeDetectionFunction::process(const TCSGram& rTCSGram) { distance += std::pow(nextTCS[j] - previousTCS[j], 2.0); } - + retVal[iPosition] = std::pow(distance, 0.5); } diff --git a/libs/qm-dsp/dsp/tonal/ChangeDetectionFunction.h b/libs/qm-dsp/dsp/tonal/ChangeDetectionFunction.h index 3a84b3096f..5d041ad769 100644 --- a/libs/qm-dsp/dsp/tonal/ChangeDetectionFunction.h +++ b/libs/qm-dsp/dsp/tonal/ChangeDetectionFunction.h @@ -38,7 +38,7 @@ public: ChangeDistance process(const TCSGram& rTCSGram); private: void setFilterWidth(const int iWidth); - + private: valarray m_vaGaussian; double m_dFilterSigma; diff --git a/libs/qm-dsp/dsp/tonal/TCSgram.cpp b/libs/qm-dsp/dsp/tonal/TCSgram.cpp index c226c81402..954ba03e9b 100644 --- a/libs/qm-dsp/dsp/tonal/TCSgram.cpp +++ b/libs/qm-dsp/dsp/tonal/TCSgram.cpp @@ -55,7 +55,7 @@ void TCSGram::addTCSVector(const TCSVector& rTCSVector) std::pair p; p.first = lMilliSeconds; p.second = rTCSVector; - + m_VectorList.push_back(p); } @@ -68,7 +68,7 @@ long TCSGram::getDuration() const void TCSGram::printDebug() { vectorlist_t::iterator vectorIterator = m_VectorList.begin(); - + while (vectorIterator != m_VectorList.end()) { vectorIterator->second.printDebug(); diff --git a/libs/qm-dsp/dsp/tonal/TCSgram.h b/libs/qm-dsp/dsp/tonal/TCSgram.h index 83e8c93f8b..f4825a996a 100644 --- a/libs/qm-dsp/dsp/tonal/TCSgram.h +++ b/libs/qm-dsp/dsp/tonal/TCSgram.h @@ -26,7 +26,7 @@ typedef std::vector > vectorlist_t; class TCSGram { -public: +public: TCSGram(); ~TCSGram(); void getTCSVector(int, TCSVector&) const; diff --git a/libs/qm-dsp/dsp/tonal/TonalEstimator.cpp b/libs/qm-dsp/dsp/tonal/TonalEstimator.cpp index 16d1aa8995..72b6f85c83 100644 --- a/libs/qm-dsp/dsp/tonal/TonalEstimator.cpp +++ b/libs/qm-dsp/dsp/tonal/TonalEstimator.cpp @@ -27,15 +27,15 @@ TonalEstimator::TonalEstimator() m_Basis.resize(6); int i = 0; - - + + // circle of fifths m_Basis[i].resize(12); for (int iP = 0; iP < 12; iP++) { m_Basis[i][iP] = std::sin( (7.0 / 6.0) * iP * PI); } - + i++; m_Basis[i].resize(12); @@ -43,17 +43,17 @@ TonalEstimator::TonalEstimator() { m_Basis[i][iP] = std::cos( (7.0 / 6.0) * iP * PI); } - + i++; - - + + // circle of major thirds m_Basis[i].resize(12); for (int iP = 0; iP < 12; iP++) { m_Basis[i][iP] = 0.6 * std::sin( (2.0 / 3.0) * iP * PI); } - + i++; m_Basis[i].resize(12); @@ -71,7 +71,7 @@ TonalEstimator::TonalEstimator() { m_Basis[i][iP] = 1.1 * std::sin( (3.0 / 2.0) * iP * PI); } - + i++; m_Basis[i].resize(12); @@ -90,7 +90,7 @@ TCSVector TonalEstimator::transform2TCS(const ChromaVector& rVector) { TCSVector vaRetVal; vaRetVal.resize(6, 0.0); - + for (int i = 0; i < 6; i++) { for (int iP = 0; iP < 12; iP++) @@ -98,6 +98,6 @@ TCSVector TonalEstimator::transform2TCS(const ChromaVector& rVector) vaRetVal[i] += m_Basis[i][iP] * rVector[iP]; } } - + return vaRetVal; } diff --git a/libs/qm-dsp/dsp/tonal/TonalEstimator.h b/libs/qm-dsp/dsp/tonal/TonalEstimator.h index 5753dff050..cfb8bba5b6 100644 --- a/libs/qm-dsp/dsp/tonal/TonalEstimator.h +++ b/libs/qm-dsp/dsp/tonal/TonalEstimator.h @@ -27,24 +27,24 @@ class ChromaVector : public std::valarray public: ChromaVector(size_t uSize = 12) : std::valarray() { resize(uSize, 0.0f); } - + virtual ~ChromaVector() {}; - + void printDebug() { for (int i = 0; i < size(); i++) { std::cout << (*this)[i] << ";"; } - + std::cout << std::endl; } - + void normalizeL1() { // normalize the chroma vector (L1 norm) double dSum = 0.0; - + for (size_t i = 0; i < 12; (dSum += std::abs((*this)[i++]))) ; for (size_t i = 0; i < 12; dSum > 0.0000001?((*this)[i] /= dSum):(*this)[i]=0.0, i++) ; @@ -55,7 +55,7 @@ public: for (size_t i = 0; i < 12; ++i) (*this)[i] = 0.0; } - + }; class TCSVector : public std::valarray @@ -63,7 +63,7 @@ class TCSVector : public std::valarray public: TCSVector() : std::valarray() { resize(6, 0.0f); } - + virtual ~TCSVector() {}; void printDebug() @@ -72,19 +72,19 @@ public: { std::cout << (*this)[i] << ";"; } - + std::cout << std::endl; } - + double magnitude() const { double dMag = 0.0; - + for (size_t i = 0; i < 6; i++) { dMag += std::pow((*this)[i], 2.0); } - + return std::sqrt(dMag); } diff --git a/libs/qm-dsp/dsp/transforms/FFT.cpp b/libs/qm-dsp/dsp/transforms/FFT.cpp index 454cfb1422..0ca618f40b 100644 --- a/libs/qm-dsp/dsp/transforms/FFT.cpp +++ b/libs/qm-dsp/dsp/transforms/FFT.cpp @@ -54,7 +54,7 @@ FFTReal::process(bool inverse, } static unsigned int numberOfBitsNeeded(unsigned int p_nSamples) -{ +{ int i; if( p_nSamples < 2 ) diff --git a/libs/qm-dsp/dsp/wavelet/Wavelet.cpp b/libs/qm-dsp/dsp/wavelet/Wavelet.cpp index d3e1fa1baf..57d6337719 100644 --- a/libs/qm-dsp/dsp/wavelet/Wavelet.cpp +++ b/libs/qm-dsp/dsp/wavelet/Wavelet.cpp @@ -82,7 +82,7 @@ Wavelet::createDecompositionFilters(Type wavelet, hpd.clear(); unsigned int flength = 0; - + switch (wavelet) { case Haar: @@ -103,7 +103,7 @@ Wavelet::createDecompositionFilters(Type wavelet, hpd.push_back(-0.22414386804186); hpd.push_back(-0.12940952255092); flength = 4; - break; + break; case Daubechies_3: lpd.push_back(0.03522629188210); @@ -592,7 +592,7 @@ Wavelet::createDecompositionFilters(Type wavelet, hpd.push_back(-0.00000000000000); flength = 80; break; - + case Symlet_2: lpd.push_back(-0.12940952255092); lpd.push_back(0.22414386804186); @@ -692,7 +692,7 @@ Wavelet::createDecompositionFilters(Type wavelet, hpd.push_back(0.01540410932703); flength = 12; break; - + case Symlet_7: lpd.push_back(0.00268181456826); lpd.push_back(-0.00104738488868); diff --git a/libs/qm-dsp/hmm/hmm.c b/libs/qm-dsp/hmm/hmm.c index 69eee02b66..6642e4e1db 100644 --- a/libs/qm-dsp/hmm/hmm.c +++ b/libs/qm-dsp/hmm/hmm.c @@ -46,11 +46,11 @@ 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->L = L; model->p0 = (double*) malloc(N*sizeof(double)); model->a = (double**) malloc(N*sizeof(double*)); model->mu = (double**) malloc(N*sizeof(double*)); @@ -62,10 +62,10 @@ model_t* hmm_init(double** x, int T, int L, int N) 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++) { @@ -74,7 +74,7 @@ model_t* hmm_init(double** x, int T, int L, int N) global_mean[d] += x[t][d]; global_mean[d] /= T; } - + /* calculate global diagonal covariance */ for (d = 0; d < L; d++) { @@ -84,7 +84,7 @@ model_t* hmm_init(double** x, int T, int L, int N) 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++) { @@ -94,8 +94,8 @@ model_t* hmm_init(double** x, int T, int L, int N) /* 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++) @@ -115,16 +115,16 @@ model_t* hmm_init(double** x, int T, int L, int N) } 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]); @@ -134,23 +134,23 @@ void hmm_close(model_t* model) free(model->mu); for (i = 0; i < model->L; i++) free(model->cov[i]); - free(model->cov); - free(model); + 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 */ + + /* allocate memory */ double** gamma = (double**) malloc(T*sizeof(double*)); double*** xi = (double***) malloc(T*sizeof(double**)); for (t = 0; t < T; t++) @@ -160,45 +160,45 @@ void hmm_train(double** x, int T, model_t* model) 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)); - + 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 niter = 50; int iter = 0; double loglik1, loglik2; int foundnan = 0; - while (iter < niter && !foundnan && !(iter > 1 && (loglik - loglik1) < thresh * (loglik1 - loglik2))) + 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; } @@ -213,13 +213,13 @@ void hmm_train(double** x, int T, model_t* model) } */ } -/* +/* 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, "iteration %d: loglik = %f\n", iter, loglik); fprintf(stderr, "re-estimation...\n"); fflush(stderr); */ @@ -227,9 +227,9 @@ void hmm_train(double** x, int T, model_t* model) 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++) @@ -242,7 +242,7 @@ void hmm_train(double** x, int T, model_t* model) */ //hmm_print(model); } - + /* deallocate memory */ for (t = 0; t < T; t++) { @@ -254,12 +254,12 @@ void hmm_train(double** x, int T, model_t* model) } free(gamma); free(xi); - free(b); - + free(b); + for (i = 0; i < L; i++) free(icov[i]); free(icov); - + free(gauss_y); free(gauss_z); } @@ -267,27 +267,27 @@ void hmm_train(double** x, int T, model_t* model) void mlss_reestimate(double* p0, double** a, double** mu, double** cov, int N, int T, int L, int* q, double** x) { /* fit a single Gaussian to observations in each state */ - + /* calculate the mean observation in each state */ - + /* calculate the overall covariance */ - + /* count transitions */ - + /* estimate initial probs from transitions (???) */ } 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)); - + double* yy2 = (double*) malloc(T*L*sizeof(double)); + /* re-estimate transition probs */ for (i = 0; i < N; i++) { @@ -295,7 +295,7 @@ void baum_welch(double* p0, double** a, double** mu, double** cov, int N, int T, for (t = 0; t < T-1; t++) sum_gamma[i] += gamma[t][i]; } - + for (i = 0; i < N; i++) { if (sum_gamma[i] == 0) @@ -310,7 +310,7 @@ void baum_welch(double* p0, double** a, double** mu, double** cov, int N, int T, for (t = 0; t < T-1; t++) a[i][j] += xi[t][i][j]; //s += a[i][j]; - a[i][j] /= sum_gamma[i]; + a[i][j] /= sum_gamma[i]; } /* for (j = 0; j < N; j++) @@ -319,21 +319,21 @@ void baum_welch(double* p0, double** a, double** mu, double** cov, int N, int T, } */ } - + /* 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]; - + sum_sum_gamma += sum_gamma[i]; + /* for (d = 0; d < L; d++) { @@ -343,9 +343,9 @@ void baum_welch(double* p0, double** a, double** mu, double** cov, int N, int T, 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); @@ -365,12 +365,12 @@ void baum_welch(double* p0, double** a, double** mu, double** cov, int N, int T, 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++) @@ -379,20 +379,20 @@ void baum_welch(double* p0, double** a, double** mu, double** cov, int N, int 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; */ - + 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++) { @@ -404,7 +404,7 @@ void baum_welch(double* p0, double** a, double** mu, double** cov, int N, int T, mu[j][d] /= sum_gamma[j]; } } - + /* deallocate memory */ free(sum_gamma); free(yy); @@ -416,7 +416,7 @@ void forward_backwards(double*** xi, double** gamma, double* loglik, double* log { /* 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++) @@ -424,34 +424,34 @@ void forward_backwards(double*** xi, double** gamma, double* loglik, double* log 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]; - + //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]; - + 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++) { @@ -459,10 +459,10 @@ void forward_backwards(double*** xi, double** gamma, double* loglik, double* log 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) { @@ -477,16 +477,16 @@ void forward_backwards(double*** xi, double** gamma, double* loglik, double* log 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; @@ -495,20 +495,20 @@ void forward_backwards(double*** xi, double** gamma, double* loglik, double* log { 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++) @@ -526,7 +526,7 @@ void forward_backwards(double*** xi, double** gamma, double* loglik, double* log } printf("\n\n"); */ - + /* calculate posterior probs */ double tot; for (t = 0; t < T; t++) @@ -539,12 +539,12 @@ void forward_backwards(double*** xi, double** gamma, double* loglik, double* log } for (i = 0; i < N; i++) { - gamma[t][i] /= tot; - - //printf("gamma[%d][%d] = %f\n", t, i, gamma[t][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; @@ -560,7 +560,7 @@ void forward_backwards(double*** xi, double** gamma, double* loglik, double* log for (j = 0; j < N; j++) xi[t][i][j] /= tot; } - + /* // CHECK - fine // gamma[t][i] = \sum_j{xi[t][i][j]} @@ -568,8 +568,8 @@ void forward_backwards(double*** xi, double** gamma, double* loglik, double* log 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]); @@ -585,20 +585,20 @@ void viterbi_decode(double** x, int T, model_t* model, int* q) 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; - + /* 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*)); @@ -608,24 +608,24 @@ void viterbi_decode(double** x, int T, model_t* model, int* q) 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)); - + 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); - + /* 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++) @@ -646,7 +646,7 @@ void viterbi_decode(double** x, int T, model_t* model, int* q) } } } - + /* find maximising state at time T-1 */ max = phi[T-1][0]; q[T-1] = 0; @@ -659,7 +659,7 @@ void viterbi_decode(double** x, int T, model_t* model, int* q) } } - + /* track back */ t = T - 2; while (t >= 0) @@ -667,7 +667,7 @@ void viterbi_decode(double** x, int T, model_t* model, int* q) q[t] = psi[t+1][q[t+1]]; t--; } - + /* de-allocate memory */ for (i = 0; i < L; i++) free(icov[i]); @@ -681,7 +681,7 @@ void viterbi_decode(double** x, int T, model_t* model, int* q) free(logb); free(phi); free(psi); - + free(gauss_y); free(gauss_z); } @@ -695,11 +695,11 @@ void invert(double** cov, int L, double** icov, double* detcov) for(j=0; j < L; j++) for (i=0; i < L; i++) a[j*L+i] = cov[i][j]; - - int M = (int) L; + + 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) @@ -707,7 +707,7 @@ void invert(double** cov, int L, double** icov, double* detcov) fprintf(stderr, "Covariance matrix was singular, couldn't invert\n"); exit(-1); } - + /* find determinant */ double det = 1; for(i = 0; i < L; i++) @@ -723,27 +723,27 @@ void invert(double** cov, int L, double** icov, double* detcov) if (det < 0) det = -det; *detcov = det; - + /* 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; + lwork = (int) lwbest; double* work = (double*) malloc(lwork*sizeof(double)); #endif - + /* 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]; - -#ifndef HAVE_ATLAS + icov[i][j] = a[j*L+i]; + +#ifndef HAVE_ATLAS free(work); #endif - free(a); + free(a); } /* probability of multivariate Gaussian given mean, inverse and determinant of covariance */ @@ -762,8 +762,8 @@ double gauss(double* x, int L, double* mu, double** icov, double detcov, double* } s = cblas_ddot(L, z, 1, y, 1); //for (i = 0; i < L; i++) - // s += z[i] * y[i]; - + // s += z[i] * y[i]; + return exp(-s/2.0) / (pow(2*PI, L/2.0) * sqrt(detcov)); } @@ -784,10 +784,10 @@ double loggauss(double* x, int L, double* mu, double** icov, double detcov, doub } s = cblas_ddot(L, z, 1, y, 1); //for (i = 0; i < L; i++) - // s += z[i] * y[i]; - + // s += z[i] * y[i]; + ret = -0.5 * (s + L * log(2*PI) + log(detcov)); - + /* // TEST if (ISINF(ret) > 0) @@ -795,9 +795,9 @@ double loggauss(double* x, int L, double* mu, double** icov, double detcov, doub if (ISINF(ret) < 0) printf("loggauss returning -infinity\n"); if (ISNAN(ret)) - printf("loggauss returning nan\n"); + printf("loggauss returning nan\n"); */ - + return ret; } diff --git a/libs/qm-dsp/maths/MathUtilities.cpp b/libs/qm-dsp/maths/MathUtilities.cpp index 2249c53a95..a55252b365 100644 --- a/libs/qm-dsp/maths/MathUtilities.cpp +++ b/libs/qm-dsp/maths/MathUtilities.cpp @@ -41,11 +41,11 @@ void MathUtilities::getAlphaNorm(const double *data, unsigned int len, unsigned unsigned int i; double temp = 0.0; double a=0.0; - + for( i = 0; i < len; i++) { temp = data[ i ]; - + a += ::pow( fabs(temp), double(alpha) ); } a /= ( double )len; @@ -60,11 +60,11 @@ double MathUtilities::getAlphaNorm( const std::vector &data, unsigned i unsigned int len = data.size(); double temp = 0.0; double a=0.0; - + for( i = 0; i < len; i++) { temp = data[ i ]; - + a += ::pow( fabs(temp), double(alpha) ); } a /= ( double )len; @@ -143,7 +143,7 @@ double MathUtilities::mean(const double *src, unsigned int len) double retVal =0.0; double s = sum( src, len ); - + retVal = s / (double)len; return retVal; @@ -154,7 +154,7 @@ double MathUtilities::mean(const std::vector &src, unsigned int count) { double sum = 0.; - + for (unsigned int i = 0; i < count; ++i) { sum += src[start + i]; @@ -172,7 +172,7 @@ void MathUtilities::getFrameMinMax(const double *data, unsigned int len, double *min = *max = 0; return; } - + *min = data[0]; *max = data[0]; @@ -188,7 +188,7 @@ void MathUtilities::getFrameMinMax(const double *data, unsigned int len, double { *max = temp ; } - + } } @@ -197,7 +197,7 @@ int MathUtilities::getMax( double* pData, unsigned int Length, double* pMax ) unsigned int index = 0; unsigned int i; double temp = 0.0; - + double max = pData[0]; for( i = 0; i < Length; i++) @@ -209,7 +209,7 @@ int MathUtilities::getMax( double* pData, unsigned int Length, double* pMax ) max = temp ; index = i; } - + } if (pMax) *pMax = max; @@ -223,7 +223,7 @@ int MathUtilities::getMax( const std::vector & data, double* pMax ) unsigned int index = 0; unsigned int i; double temp = 0.0; - + double max = data[0]; for( i = 0; i < data.size(); i++) @@ -235,7 +235,7 @@ int MathUtilities::getMax( const std::vector & data, double* pMax ) max = temp ; index = i; } - + } if (pMax) *pMax = max; @@ -344,7 +344,7 @@ void MathUtilities::adaptiveThreshold(std::vector &data) if (sz == 0) return; std::vector smoothed(sz); - + int p_pre = 8; int p_post = 7; diff --git a/libs/qm-dsp/maths/MathUtilities.h b/libs/qm-dsp/maths/MathUtilities.h index 1d3395e8ad..85774756d4 100644 --- a/libs/qm-dsp/maths/MathUtilities.h +++ b/libs/qm-dsp/maths/MathUtilities.h @@ -22,7 +22,7 @@ class MathUtilities { -public: +public: static double round( double x ); static void getFrameMinMax( const double* data, unsigned int len, double* min, double* max ); diff --git a/libs/qm-dsp/maths/pca/pca.c b/libs/qm-dsp/maths/pca/pca.c index 288ec4d2e4..d2b2f11c47 100644 --- a/libs/qm-dsp/maths/pca/pca.c +++ b/libs/qm-dsp/maths/pca/pca.c @@ -110,7 +110,7 @@ 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--) { l = i - 1; @@ -188,7 +188,7 @@ 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; @@ -253,19 +253,19 @@ 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 **********************************************************************/ - + /* 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 */ @@ -278,7 +278,7 @@ void pca_project(double** data, int n, int m, int ncomponents) tred2(symmat, m, evals, interm); /* Triangular decomposition */ 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"); @@ -289,7 +289,7 @@ columns of symmat now contain the associated eigenvectors. */ 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++) { @@ -310,7 +310,7 @@ for (i = 0; i < n; i++) { } } -/* +/* printf("\nProjections of row-points on first 3 prin. comps.:\n"); for (i = 0; i < n; i++) { for (j = 0; j < 3; j++) { -- cgit v1.2.3