Mercurial > hg > qm-dsp
changeset 477:fa407c1d9923
Untabify + indent
| author | Chris Cannam <cannam@all-day-breakfast.com> |
|---|---|
| date | Thu, 30 May 2019 18:36:48 +0100 |
| parents | 2de6184b2ce0 |
| children | c92718cc6ef1 |
| files | maths/Correlation.cpp maths/MathAliases.h maths/MathUtilities.cpp maths/MathUtilities.h maths/MedianFilter.h maths/Polyfit.h |
| diffstat | 6 files changed, 260 insertions(+), 275 deletions(-) [+] |
line wrap: on
line diff
--- a/maths/Correlation.cpp Thu May 30 18:30:58 2019 +0100 +++ b/maths/Correlation.cpp Thu May 30 18:36:48 2019 +0100 @@ -36,21 +36,19 @@ unsigned int i,j; - for( i = 0; i < length; i++) - { - for( j = i; j < length; j++) - { - tmp += src[ j-i ] * src[ j ]; - } + for( i = 0; i < length; i++) { + for( j = i; j < length; j++) { + tmp += src[ j-i ] * src[ j ]; + } + outVal = tmp / ( length - i ); - outVal = tmp / ( length - i ); - - if( outVal <= 0 ) - dst[ i ] = EPS; - else - dst[ i ] = outVal; - - tmp = 0.0; + if( outVal <= 0 ) { + dst[ i ] = EPS; + } else { + dst[ i ] = outVal; + } + + tmp = 0.0; } }
--- a/maths/MathAliases.h Thu May 30 18:30:58 2019 +0100 +++ b/maths/MathAliases.h Thu May 30 18:36:48 2019 +0100 @@ -27,34 +27,34 @@ #define PI (3.14159265358979232846) #endif -#define TWO_PI (2. * PI) +#define TWO_PI (2. * PI) #define EPS 2.2204e-016 /* aliases to math.h functions */ -#define EXP exp -#define COS cos -#define SIN sin -#define ABS fabs -#define POW powf -#define SQRT sqrtf -#define LOG10 log10f -#define LOG logf -#define FLOOR floorf -#define TRUNC truncf +#define EXP exp +#define COS cos +#define SIN sin +#define ABS fabs +#define POW powf +#define SQRT sqrtf +#define LOG10 log10f +#define LOG logf +#define FLOOR floorf +#define TRUNC truncf /* aliases to complex.h functions */ /** sample = EXPC(complex) */ -#define EXPC cexpf +#define EXPC cexpf /** complex = CEXPC(complex) */ -#define CEXPC cexp +#define CEXPC cexp /** sample = ARGC(complex) */ -#define ARGC cargf +#define ARGC cargf /** sample = ABSC(complex) norm */ -#define ABSC cabsf +#define ABSC cabsf /** sample = REAL(complex) */ -#define REAL crealf +#define REAL crealf /** sample = IMAG(complex) */ -#define IMAG cimagf +#define IMAG cimagf #endif
--- a/maths/MathUtilities.cpp Thu May 30 18:30:58 2019 +0100 +++ b/maths/MathUtilities.cpp Thu May 30 18:36:48 2019 +0100 @@ -44,12 +44,10 @@ 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) ); + + for( i = 0; i < len; i++) { + temp = data[ i ]; + a += ::pow( fabs(temp), double(alpha) ); } a /= ( double )len; a = ::pow( a, ( 1.0 / (double) alpha ) ); @@ -63,11 +61,10 @@ 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) ); + + for( i = 0; i < len; i++) { + temp = data[ i ]; + a += ::pow( fabs(temp), double(alpha) ); } a /= ( double )len; a = ::pow( a, ( 1.0 / (double) alpha ) ); @@ -105,9 +102,8 @@ int i ; double retVal =0.0; - for( i = 0; i < len; i++) - { - retVal += src[ i ]; + for( i = 0; i < len; i++) { + retVal += src[ i ]; } return retVal; @@ -120,7 +116,7 @@ if (len == 0) return 0; double s = sum( src, len ); - + retVal = s / (double)len; return retVal; @@ -131,11 +127,10 @@ int count) { double sum = 0.; - + if (count == 0) return 0; - for (int i = 0; i < (int)count; ++i) - { + for (int i = 0; i < (int)count; ++i) { sum += src[start + i]; } @@ -151,92 +146,82 @@ *min = *max = 0; return; } - + *min = data[0]; *max = data[0]; - for( i = 0; i < len; i++) - { - temp = data[ i ]; + for( i = 0; i < len; i++) { + temp = data[ i ]; - if( temp < *min ) - { - *min = temp ; - } - if( temp > *max ) - { - *max = temp ; - } - + if( temp < *min ) { + *min = temp ; + } + if( temp > *max ) { + *max = temp ; + } } } int MathUtilities::getMax( double* pData, int Length, double* pMax ) { - int index = 0; - int i; - double temp = 0.0; - - double max = pData[0]; + int index = 0; + int i; + double temp = 0.0; + + double max = pData[0]; - for( i = 0; i < Length; i++) - { - temp = pData[ i ]; + for( i = 0; i < Length; i++) { + temp = pData[ i ]; - if( temp > max ) - { - max = temp ; - index = i; - } - - } + if( temp > max ) { + max = temp ; + index = i; + } + } - if (pMax) *pMax = max; + if (pMax) *pMax = max; - return index; + return index; } int MathUtilities::getMax( const vector<double> & data, double* pMax ) { - int index = 0; - int i; - double temp = 0.0; - - double max = data[0]; + int index = 0; + int i; + double temp = 0.0; + + double max = data[0]; - for( i = 0; i < int(data.size()); i++) - { - temp = data[ i ]; + for( i = 0; i < int(data.size()); i++) { - if( temp > max ) - { - max = temp ; - index = i; - } - - } + temp = data[ i ]; - if (pMax) *pMax = max; + if( temp > max ) { + max = temp ; + index = i; + } + } + if (pMax) *pMax = max; - return index; + + return index; } void MathUtilities::circShift( double* pData, int length, int shift) { - shift = shift % length; - double temp; - int i,n; + shift = shift % length; + double temp; + int i,n; - for( i = 0; i < shift; i++) - { - temp=*(pData + length - 1); + for( i = 0; i < shift; i++) { + + temp=*(pData + length - 1); - for( n = length-2; n >= 0; n--) - { - *(pData+n+1)=*(pData+n); - } + for( n = length-2; n >= 0; n--) { + *(pData+n+1)=*(pData+n); + } *pData = temp; } @@ -244,7 +229,7 @@ int MathUtilities::compareInt (const void * a, const void * b) { - return ( *(int*)a - *(int*)b ); + return ( *(int*)a - *(int*)b ); } void MathUtilities::normalise(double *data, int length, NormaliseType type) @@ -327,8 +312,8 @@ } vector<double> MathUtilities::normaliseLp(const vector<double> &data, - int p, - double threshold) + int p, + double threshold) { int n = int(data.size()); if (n == 0 || p == 0) return data; @@ -349,7 +334,7 @@ if (sz == 0) return; vector<double> smoothed(sz); - + int p_pre = 8; int p_post = 7; @@ -400,7 +385,8 @@ MathUtilities::nearestPowerOfTwo(int x) { if (isPowerOfTwo(x)) return x; - int n0 = previousPowerOfTwo(x), n1 = nextPowerOfTwo(x); + int n0 = previousPowerOfTwo(x); + int n1 = nextPowerOfTwo(x); if (x - n0 < n1 - x) return n0; else return n1; } @@ -411,7 +397,7 @@ if (x < 0) return 0; double f = 1; for (int i = 1; i <= x; ++i) { - f = f * i; + f = f * i; } return f; }
--- a/maths/MathUtilities.h Thu May 30 18:30:58 2019 +0100 +++ b/maths/MathUtilities.h Thu May 30 18:36:48 2019 +0100 @@ -25,7 +25,7 @@ */ class MathUtilities { -public: +public: /** * Round x to the nearest integer. */ @@ -35,7 +35,8 @@ * Return through min and max pointers the highest and lowest * values in the given array of the given length. */ - static void getFrameMinMax( const double* data, int len, double* min, double* max ); + static void getFrameMinMax( const double* data, int len, + double* min, double* max ); /** * Return the mean of the given array of the given length. @@ -79,7 +80,7 @@ * of the input data, and when alpha = 1 this is the mean * magnitude. */ - static void getAlphaNorm(const double *data, int len, int alpha, double* ANorm); + static void getAlphaNorm(const double *data, int len, int alpha, double* ANorm); /** * The alpha norm is the alpha'th root of the mean alpha'th power @@ -124,11 +125,11 @@ */ static void adaptiveThreshold(std::vector<double> &data); - static void circShift( double* data, int length, int shift); + static void circShift( double* data, int length, int shift); - static int getMax( double* data, int length, double* max = 0 ); - static int getMax( const std::vector<double> &data, double* max = 0 ); - static int compareInt(const void * a, const void * b); + static int getMax( double* data, int length, double* max = 0 ); + static int getMax( const std::vector<double> &data, double* max = 0 ); + static int compareInt(const void * a, const void * b); /** * Return true if x is 2^n for some integer n >= 0.
--- a/maths/MedianFilter.h Thu May 30 18:30:58 2019 +0100 +++ b/maths/MedianFilter.h Thu May 30 18:36:48 2019 +0100 @@ -27,17 +27,17 @@ { public: MedianFilter(int size, float percentile = 50.f) : - m_size(size), - m_frame(new T[size]), - m_sorted(new T[size]), - m_sortend(m_sorted + size - 1) { + m_size(size), + m_frame(new T[size]), + m_sorted(new T[size]), + m_sortend(m_sorted + size - 1) { setPercentile(percentile); - reset(); + reset(); } ~MedianFilter() { - delete[] m_frame; - delete[] m_sorted; + delete[] m_frame; + delete[] m_sorted; } void setPercentile(float p) { @@ -52,15 +52,15 @@ // we do need to push something, to maintain the filter length value = T(); } - drop(m_frame[0]); - const int sz1 = m_size-1; - for (int i = 0; i < sz1; ++i) m_frame[i] = m_frame[i+1]; - m_frame[m_size-1] = value; - put(value); + drop(m_frame[0]); + const int sz1 = m_size-1; + for (int i = 0; i < sz1; ++i) m_frame[i] = m_frame[i+1]; + m_frame[m_size-1] = value; + put(value); } T get() const { - return m_sorted[m_index]; + return m_sorted[m_index]; } int getSize() const { @@ -68,15 +68,15 @@ } T getAt(float percentile) { - int ix = int((m_size * percentile) / 100.f); + int ix = int((m_size * percentile) / 100.f); if (ix >= m_size) ix = m_size-1; if (ix < 0) ix = 0; - return m_sorted[ix]; + return m_sorted[ix]; } void reset() { - for (int i = 0; i < m_size; ++i) m_frame[i] = 0; - for (int i = 0; i < m_size; ++i) m_sorted[i] = 0; + for (int i = 0; i < m_size; ++i) m_frame[i] = 0; + for (int i = 0; i < m_size; ++i) m_sorted[i] = 0; } static std::vector<T> filter(int size, const std::vector<T> &in) { @@ -102,25 +102,25 @@ int m_index; void put(T value) { - // precondition: m_sorted contains m_size-1 values, packed at start - // postcondition: m_sorted contains m_size values, one of which is value - T *point = std::lower_bound(m_sorted, m_sortend, value); - const int n = m_sortend - point; - for (int i = n; i > 0; --i) point[i] = point[i-1]; - *point = value; + // precondition: m_sorted contains m_size-1 values, packed at start + // postcondition: m_sorted contains m_size values, one of which is value + T *point = std::lower_bound(m_sorted, m_sortend, value); + const int n = m_sortend - point; + for (int i = n; i > 0; --i) point[i] = point[i-1]; + *point = value; } void drop(T value) { - // precondition: m_sorted contains m_size values, one of which is value - // postcondition: m_sorted contains m_size-1 values, packed at start - T *point = std::lower_bound(m_sorted, m_sortend + 1, value); - if (*point != value) { - std::cerr << "WARNING: MedianFilter::drop: *point is " << *point - << ", expected " << value << std::endl; - } - const int n = m_sortend - point; - for (int i = 0; i < n; ++i) point[i] = point[i+1]; - *m_sortend = T(0); + // precondition: m_sorted contains m_size values, one of which is value + // postcondition: m_sorted contains m_size-1 values, packed at start + T *point = std::lower_bound(m_sorted, m_sortend + 1, value); + if (*point != value) { + std::cerr << "WARNING: MedianFilter::drop: *point is " << *point + << ", expected " << value << std::endl; + } + const int n = m_sortend - point; + for (int i = 0; i < n; ++i) point[i] = point[i+1]; + *m_sortend = T(0); } MedianFilter(const MedianFilter &); // not provided
--- a/maths/Polyfit.h Thu May 30 18:30:58 2019 +0100 +++ b/maths/Polyfit.h Thu May 30 18:36:48 2019 +0100 @@ -50,8 +50,8 @@ public: static double PolyFit2 (const vector<double> &x, // does the work - const vector<double> &y, - vector<double> &coef); + const vector<double> &y, + vector<double> &coef); private: @@ -61,21 +61,21 @@ static void Square (const Matrix &x, // Matrix multiplication routine - const vector<double> &y, - Matrix &a, // A = transpose X times X - vector<double> &g, // G = Y times X - const int nrow, const int ncol); + const vector<double> &y, + Matrix &a, // A = transpose X times X + vector<double> &g, // G = Y times X + const int nrow, const int ncol); // Forms square coefficient matrix static bool GaussJordan (Matrix &b, // square matrix of coefficients - const vector<double> &y, // constant vector - vector<double> &coef); // solution vector + const vector<double> &y, // constant vector + vector<double> &coef); // solution vector // returns false if matrix singular static bool GaussJordan2(Matrix &b, - const vector<double> &y, - Matrix &w, - vector<vector<int> > &index); + const vector<double> &y, + Matrix &w, + vector<vector<int> > &index); }; // some utility functions @@ -114,8 +114,8 @@ // main PolyFit routine double TPolyFit::PolyFit2 (const vector<double> &x, - const vector<double> &y, - vector<double> &coefs) + const vector<double> &y, + vector<double> &coefs) // nterms = coefs.size() // npoints = x.size() { @@ -144,38 +144,38 @@ } for(i = 0; i < npoints; ++i) { - // { setup x matrix } - xi = x[i]; - xmatr[i][0] = 1.0; // { first column } - for(j = 1; j < nterms; ++j) - xmatr[i][j] = xmatr [i][j - 1] * xi; + // { setup x matrix } + xi = x[i]; + xmatr[i][0] = 1.0; // { first column } + for(j = 1; j < nterms; ++j) + xmatr[i][j] = xmatr [i][j - 1] * xi; } Square (xmatr, y, a, g, npoints, nterms); if(!GaussJordan (a, g, coefs)) - return -1; + return -1; sum_y = 0.0; sum_y2 = 0.0; srs = 0.0; for(i = 0; i < npoints; ++i) { - yi = y[i]; - yc = 0.0; - for(j = 0; j < nterms; ++j) - yc += coefs [j] * xmatr [i][j]; - srs += NSUtility::sqr (yc - yi); - sum_y += yi; - sum_y2 += yi * yi; + yi = y[i]; + yc = 0.0; + for(j = 0; j < nterms; ++j) + yc += coefs [j] * xmatr [i][j]; + srs += NSUtility::sqr (yc - yi); + sum_y += yi; + sum_y2 += yi * yi; } // If all Y values are the same, avoid dividing by zero correl_coef = sum_y2 - NSUtility::sqr (sum_y) / npoints; // Either return 0 or the correct value of correlation coefficient if (correl_coef != 0) - correl_coef = srs / correl_coef; + correl_coef = srs / correl_coef; if (correl_coef >= 1) - correl_coef = 0.0; + correl_coef = 0.0; else - correl_coef = sqrt (1.0 - correl_coef); + correl_coef = sqrt (1.0 - correl_coef); return correl_coef; } @@ -189,36 +189,36 @@ // Form square coefficient matrix void TPolyFit::Square (const Matrix &x, - const vector<double> &y, - Matrix &a, - vector<double> &g, - const int nrow, - const int ncol) + const vector<double> &y, + Matrix &a, + vector<double> &g, + const int nrow, + const int ncol) { int i, k, l; for(k = 0; k < ncol; ++k) { - for(l = 0; l < k + 1; ++l) - { - a [k][l] = 0.0; - for(i = 0; i < nrow; ++i) - { - a[k][l] += x[i][l] * x [i][k]; - if(k != l) - a[l][k] = a[k][l]; - } - } - g[k] = 0.0; - for(i = 0; i < nrow; ++i) - g[k] += y[i] * x[i][k]; + for(l = 0; l < k + 1; ++l) + { + a [k][l] = 0.0; + for(i = 0; i < nrow; ++i) + { + a[k][l] += x[i][l] * x [i][k]; + if(k != l) + a[l][k] = a[k][l]; + } + } + g[k] = 0.0; + for(i = 0; i < nrow; ++i) + g[k] += y[i] * x[i][k]; } } //--------------------------------------------------------------------------------- bool TPolyFit::GaussJordan (Matrix &b, - const vector<double> &y, - vector<double> &coef) + const vector<double> &y, + vector<double> &coef) //b square matrix of coefficients //y constant vector //coef solution vector @@ -247,44 +247,44 @@ NSUtility::zeroise(index, ncol, 3); if(!GaussJordan2(b, y, w, index)) - return false; + return false; // Interchange columns int m; for (int i = 0; i < ncol; ++i) { - m = ncol - i - 1; - if(index [m][0] != index [m][1]) - { - irow = index [m][0]; - icol = index [m][1]; - for(int k = 0; k < ncol; ++k) - swap (b[k][irow], b[k][icol]); - } + m = ncol - i - 1; + if(index [m][0] != index [m][1]) + { + irow = index [m][0]; + icol = index [m][1]; + for(int k = 0; k < ncol; ++k) + swap (b[k][irow], b[k][icol]); + } } for(int k = 0; k < ncol; ++k) { - if(index [k][2] != 0) - { + if(index [k][2] != 0) + { std::cerr << "ERROR: Error in PolyFit::GaussJordan: matrix is singular" << std::endl; return false; - } + } } for( int i = 0; i < ncol; ++i) - coef[i] = w[i][0]; + coef[i] = w[i][0]; return true; -} // end; { procedure GaussJordan } +} // end; { procedure GaussJordan } //---------------------------------------------------------------------------------------------- bool TPolyFit::GaussJordan2(Matrix &b, - const vector<double> &y, - Matrix &w, - vector<vector<int> > &index) + const vector<double> &y, + Matrix &w, + vector<vector<int> > &index) { //GaussJordan2; // first half of GaussJordan // actual start of gaussj @@ -297,75 +297,75 @@ int nv = 1; // single constant vector for(int i = 0; i < ncol; ++i) { - w[i][0] = y[i]; // copy constant vector - index[i][2] = -1; + w[i][0] = y[i]; // copy constant vector + index[i][2] = -1; } determ = 1.0; for(int i = 0; i < ncol; ++i) { - // Search for largest element - big = 0.0; - for(int j = 0; j < ncol; ++j) - { - if(index[j][2] != 0) - { - for(int k = 0; k < ncol; ++k) - { - if(index[k][2] > 0) { + // Search for largest element + big = 0.0; + for(int j = 0; j < ncol; ++j) + { + if(index[j][2] != 0) + { + for(int k = 0; k < ncol; ++k) + { + if(index[k][2] > 0) { std::cerr << "ERROR: Error in PolyFit::GaussJordan2: matrix is singular" << std::endl; return false; } - if(index[k][2] < 0 && fabs(b[j][k]) > big) - { - irow = j; - icol = k; - big = fabs(b[j][k]); - } - } // { k-loop } - } - } // { j-loop } - index [icol][2] = index [icol][2] + 1; - index [i][0] = irow; - index [i][1] = icol; + if(index[k][2] < 0 && fabs(b[j][k]) > big) + { + irow = j; + icol = k; + big = fabs(b[j][k]); + } + } // { k-loop } + } + } // { j-loop } + index [icol][2] = index [icol][2] + 1; + index [i][0] = irow; + index [i][1] = icol; - // Interchange rows to put pivot on diagonal - // GJ3 - if(irow != icol) - { - determ = -determ; - for(int m = 0; m < ncol; ++m) - swap (b [irow][m], b[icol][m]); - if (nv > 0) - for (int m = 0; m < nv; ++m) - swap (w[irow][m], w[icol][m]); - } // end GJ3 + // Interchange rows to put pivot on diagonal + // GJ3 + if(irow != icol) + { + determ = -determ; + for(int m = 0; m < ncol; ++m) + swap (b [irow][m], b[icol][m]); + if (nv > 0) + for (int m = 0; m < nv; ++m) + swap (w[irow][m], w[icol][m]); + } // end GJ3 - // divide pivot row by pivot column - pivot = b[icol][icol]; - determ *= pivot; - b[icol][icol] = 1.0; + // divide pivot row by pivot column + pivot = b[icol][icol]; + determ *= pivot; + b[icol][icol] = 1.0; - for(int m = 0; m < ncol; ++m) - b[icol][m] /= pivot; - if(nv > 0) - for(int m = 0; m < nv; ++m) - w[icol][m] /= pivot; + for(int m = 0; m < ncol; ++m) + b[icol][m] /= pivot; + if(nv > 0) + for(int m = 0; m < nv; ++m) + w[icol][m] /= pivot; - // Reduce nonpivot rows - for(int n = 0; n < ncol; ++n) - { - if(n != icol) - { - t = b[n][icol]; - b[n][icol] = 0.0; - for(int m = 0; m < ncol; ++m) - b[n][m] -= b[icol][m] * t; - if(nv > 0) - for(int m = 0; m < nv; ++m) - w[n][m] -= w[icol][m] * t; - } - } + // Reduce nonpivot rows + for(int n = 0; n < ncol; ++n) + { + if(n != icol) + { + t = b[n][icol]; + b[n][icol] = 0.0; + for(int m = 0; m < ncol; ++m) + b[n][m] -= b[icol][m] * t; + if(nv > 0) + for(int m = 0; m < nv; ++m) + w[n][m] -= w[icol][m] * t; + } + } } // { i-loop } return true; }