c@241: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
c@241: 
c@241: /*
c@241:     QM DSP Library
c@241: 
c@241:     Centre for Digital Music, Queen Mary, University of London.
c@309:     This file 2005-2006 Christian Landone.
c@309: 
c@309:     This program is free software; you can redistribute it and/or
c@309:     modify it under the terms of the GNU General Public License as
c@309:     published by the Free Software Foundation; either version 2 of the
c@309:     License, or (at your option) any later version.  See the file
c@309:     COPYING included with this distribution for more information.
c@241: */
c@241: 
c@241: #include "MathUtilities.h"
c@241: 
c@241: #include <iostream>
c@348: #include <algorithm>
c@348: #include <vector>
c@241: #include <cmath>
c@241: 
c@418: using namespace std;
c@241: 
c@241: double MathUtilities::mod(double x, double y)
c@241: {
c@241:     double a = floor( x / y );
c@241: 
c@241:     double b = x - ( y * a );
c@241:     return b;
c@241: }
c@241: 
c@241: double MathUtilities::princarg(double ang)
c@241: {
c@241:     double ValOut;
c@241: 
c@241:     ValOut = mod( ang + M_PI, -2 * M_PI ) + M_PI;
c@241: 
c@241:     return ValOut;
c@241: }
c@241: 
c@414: void MathUtilities::getAlphaNorm(const double *data, int len, int alpha, double* ANorm)
c@241: {
c@414:     int i;
c@241:     double temp = 0.0;
c@241:     double a=0.0;
cannam@477:         
cannam@477:     for( i = 0; i < len; i++) {
cannam@477:         temp = data[ i ];
cannam@477:         a  += ::pow( fabs(temp), double(alpha) );
c@241:     }
c@241:     a /= ( double )len;
c@241:     a = ::pow( a, ( 1.0 / (double) alpha ) );
c@241: 
c@241:     *ANorm = a;
c@241: }
c@241: 
c@418: double MathUtilities::getAlphaNorm( const vector <double> &data, int alpha )
c@241: {
c@414:     int i;
c@414:     int len = data.size();
c@241:     double temp = 0.0;
c@241:     double a=0.0;
cannam@477:         
cannam@477:     for( i = 0; i < len; i++) {
cannam@477:         temp = data[ i ];
cannam@477:         a  += ::pow( fabs(temp), double(alpha) );
c@241:     }
c@241:     a /= ( double )len;
c@241:     a = ::pow( a, ( 1.0 / (double) alpha ) );
c@241: 
c@241:     return a;
c@241: }
c@241: 
c@241: double MathUtilities::round(double x)
c@241: {
c@348:     if (x < 0) {
c@348:         return -floor(-x + 0.5);
c@348:     } else {
c@348:         return floor(x + 0.5);
c@348:     }
c@241: }
c@241: 
c@414: double MathUtilities::median(const double *src, int len)
c@241: {
c@348:     if (len == 0) return 0;
c@348:     
c@418:     vector<double> scratch;
c@348:     for (int i = 0; i < len; ++i) scratch.push_back(src[i]);
c@418:     sort(scratch.begin(), scratch.end());
c@241: 
c@348:     int middle = len/2;
c@348:     if (len % 2 == 0) {
c@348:         return (scratch[middle] + scratch[middle - 1]) / 2;
c@348:     } else {
c@348:         return scratch[middle];
c@241:     }
c@241: }
c@241: 
c@414: double MathUtilities::sum(const double *src, int len)
c@241: {
c@414:     int i ;
c@241:     double retVal =0.0;
c@241: 
cannam@477:     for(  i = 0; i < len; i++) {
cannam@477:         retVal += src[ i ];
c@241:     }
c@241: 
c@241:     return retVal;
c@241: }
c@241: 
c@414: double MathUtilities::mean(const double *src, int len)
c@241: {
c@241:     double retVal =0.0;
c@241: 
c@348:     if (len == 0) return 0;
c@348: 
c@241:     double s = sum( src, len );
cannam@477:         
c@241:     retVal =  s  / (double)len;
c@241: 
c@241:     return retVal;
c@241: }
c@241: 
c@418: double MathUtilities::mean(const vector<double> &src,
c@414:                            int start,
c@414:                            int count)
c@279: {
c@279:     double sum = 0.;
cannam@477:         
c@348:     if (count == 0) return 0;
c@348:     
cannam@477:     for (int i = 0; i < (int)count; ++i) {
c@279:         sum += src[start + i];
c@279:     }
c@279: 
c@279:     return sum / count;
c@279: }
c@279: 
c@414: void MathUtilities::getFrameMinMax(const double *data, int len, double *min, double *max)
c@241: {
c@414:     int i;
c@241:     double temp = 0.0;
c@283: 
c@283:     if (len == 0) {
c@283:         *min = *max = 0;
c@283:         return;
c@283:     }
cannam@477:         
c@241:     *min = data[0];
c@241:     *max = data[0];
c@241: 
cannam@477:     for( i = 0; i < len; i++) {
cannam@477:         temp = data[ i ];
c@241: 
cannam@477:         if( temp < *min ) {
cannam@477:             *min =  temp ;
cannam@477:         }
cannam@477:         if( temp > *max ) {
cannam@477:             *max =  temp ;
cannam@477:         }
c@241:     }
c@241: }
c@241: 
c@414: int MathUtilities::getMax( double* pData, int Length, double* pMax )
c@241: {
cannam@477:     int index = 0;
cannam@477:     int i;
cannam@477:     double temp = 0.0;
cannam@477:         
cannam@477:     double max = pData[0];
c@241: 
cannam@477:     for( i = 0; i < Length; i++) {
cannam@477:         temp = pData[ i ];
c@241: 
cannam@477:         if( temp > max ) {
cannam@477:             max =  temp ;
cannam@477:             index = i;
cannam@477:         }
cannam@477:     }
c@241: 
cannam@477:     if (pMax) *pMax = max;
c@279: 
c@279: 
cannam@477:     return index;
c@279: }
c@279: 
c@418: int MathUtilities::getMax( const vector<double> & data, double* pMax )
c@279: {
cannam@477:     int index = 0;
cannam@477:     int i;
cannam@477:     double temp = 0.0;
cannam@477:         
cannam@477:     double max = data[0];
c@279: 
cannam@477:     for( i = 0; i < int(data.size()); i++) {
c@279: 
cannam@477:         temp = data[ i ];
c@279: 
cannam@477:         if( temp > max ) {
cannam@477:             max =  temp ;
cannam@477:             index = i;
cannam@477:         }
cannam@477:     }
c@241: 
cannam@477:     if (pMax) *pMax = max;
c@241: 
cannam@477: 
cannam@477:     return index;
c@241: }
c@241: 
c@241: void MathUtilities::circShift( double* pData, int length, int shift)
c@241: {
cannam@477:     shift = shift % length;
cannam@477:     double temp;
cannam@477:     int i,n;
c@241: 
cannam@477:     for( i = 0; i < shift; i++) {
cannam@477:         
cannam@477:         temp=*(pData + length - 1);
c@241: 
cannam@477:         for( n = length-2; n >= 0; n--) {
cannam@477:             *(pData+n+1)=*(pData+n);
cannam@477:         }
c@241: 
c@241:         *pData = temp;
c@241:     }
c@241: }
c@241: 
c@241: int MathUtilities::compareInt (const void * a, const void * b)
c@241: {
cannam@477:     return ( *(int*)a - *(int*)b );
c@241: }
c@241: 
c@259: void MathUtilities::normalise(double *data, int length, NormaliseType type)
c@259: {
c@259:     switch (type) {
c@259: 
c@259:     case NormaliseNone: return;
c@259: 
c@259:     case NormaliseUnitSum:
c@259:     {
c@259:         double sum = 0.0;
c@259:         for (int i = 0; i < length; ++i) {
c@259:             sum += data[i];
c@259:         }
c@259:         if (sum != 0.0) {
c@259:             for (int i = 0; i < length; ++i) {
c@259:                 data[i] /= sum;
c@259:             }
c@259:         }
c@259:     }
c@259:     break;
c@259: 
c@259:     case NormaliseUnitMax:
c@259:     {
c@259:         double max = 0.0;
c@259:         for (int i = 0; i < length; ++i) {
c@259:             if (fabs(data[i]) > max) {
c@259:                 max = fabs(data[i]);
c@259:             }
c@259:         }
c@259:         if (max != 0.0) {
c@259:             for (int i = 0; i < length; ++i) {
c@259:                 data[i] /= max;
c@259:             }
c@259:         }
c@259:     }
c@259:     break;
c@259: 
c@259:     }
c@259: }
c@259: 
c@418: void MathUtilities::normalise(vector<double> &data, NormaliseType type)
c@259: {
c@259:     switch (type) {
c@259: 
c@259:     case NormaliseNone: return;
c@259: 
c@259:     case NormaliseUnitSum:
c@259:     {
c@259:         double sum = 0.0;
c@325:         for (int i = 0; i < (int)data.size(); ++i) sum += data[i];
c@259:         if (sum != 0.0) {
c@325:             for (int i = 0; i < (int)data.size(); ++i) data[i] /= sum;
c@259:         }
c@259:     }
c@259:     break;
c@259: 
c@259:     case NormaliseUnitMax:
c@259:     {
c@259:         double max = 0.0;
c@325:         for (int i = 0; i < (int)data.size(); ++i) {
c@259:             if (fabs(data[i]) > max) max = fabs(data[i]);
c@259:         }
c@259:         if (max != 0.0) {
c@325:             for (int i = 0; i < (int)data.size(); ++i) data[i] /= max;
c@259:         }
c@259:     }
c@259:     break;
c@259: 
c@259:     }
c@259: }
c@259: 
c@418: double MathUtilities::getLpNorm(const vector<double> &data, int p)
c@418: {
c@418:     double tot = 0.0;
c@418:     for (int i = 0; i < int(data.size()); ++i) {
c@418:         tot += abs(pow(data[i], p));
c@418:     }
c@418:     return pow(tot, 1.0 / p);
c@418: }
c@418: 
c@418: vector<double> MathUtilities::normaliseLp(const vector<double> &data,
cannam@477:                                           int p,
cannam@477:                                           double threshold)
c@418: {
c@418:     int n = int(data.size());
c@418:     if (n == 0 || p == 0) return data;
c@418:     double norm = getLpNorm(data, p);
c@418:     if (norm < threshold) {
c@418:         return vector<double>(n, 1.0 / pow(n, 1.0 / p)); // unit vector
c@418:     }
c@418:     vector<double> out(n);
c@418:     for (int i = 0; i < n; ++i) {
c@418:         out[i] = data[i] / norm;
c@418:     }
c@418:     return out;
c@418: }
c@418:     
c@418: void MathUtilities::adaptiveThreshold(vector<double> &data)
c@279: {
c@279:     int sz = int(data.size());
c@279:     if (sz == 0) return;
c@279: 
c@418:     vector<double> smoothed(sz);
cannam@477:         
c@279:     int p_pre = 8;
c@279:     int p_post = 7;
c@279: 
c@279:     for (int i = 0; i < sz; ++i) {
c@279: 
c@418:         int first = max(0,      i - p_pre);
c@418:         int last  = min(sz - 1, i + p_post);
c@279: 
c@279:         smoothed[i] = mean(data, first, last - first + 1);
c@279:     }
c@279: 
c@279:     for (int i = 0; i < sz; i++) {
c@279:         data[i] -= smoothed[i];
c@279:         if (data[i] < 0.0) data[i] = 0.0;
c@279:     }
c@279: }
c@259: 
c@280: bool
c@280: MathUtilities::isPowerOfTwo(int x)
c@280: {
c@348:     if (x < 1) return false;
c@280:     if (x & (x-1)) return false;
c@280:     return true;
c@280: }
c@280: 
c@280: int
c@280: MathUtilities::nextPowerOfTwo(int x)
c@280: {
c@280:     if (isPowerOfTwo(x)) return x;
c@348:     if (x < 1) return 1;
c@280:     int n = 1;
c@280:     while (x) { x >>= 1; n <<= 1; }
c@280:     return n;
c@280: }
c@280: 
c@280: int
c@280: MathUtilities::previousPowerOfTwo(int x)
c@280: {
c@280:     if (isPowerOfTwo(x)) return x;
c@348:     if (x < 1) return 1;
c@280:     int n = 1;
c@280:     x >>= 1;
c@280:     while (x) { x >>= 1; n <<= 1; }
c@280:     return n;
c@280: }
c@280: 
c@280: int
c@280: MathUtilities::nearestPowerOfTwo(int x)
c@280: {
c@280:     if (isPowerOfTwo(x)) return x;
cannam@477:     int n0 = previousPowerOfTwo(x);
cannam@477:     int n1 = nextPowerOfTwo(x);
c@280:     if (x - n0 < n1 - x) return n0;
c@280:     else return n1;
c@280: }
c@280: 
c@360: double
c@348: MathUtilities::factorial(int x)
c@348: {
c@348:     if (x < 0) return 0;
c@360:     double f = 1;
c@348:     for (int i = 1; i <= x; ++i) {
cannam@477:         f = f * i;
c@348:     }
c@348:     return f;
c@348: }
c@348: 
c@350: int
c@350: MathUtilities::gcd(int a, int b)
c@350: {
c@350:     int c = a % b;
c@350:     if (c == 0) {
c@350:         return b;
c@350:     } else {
c@350:         return gcd(b, c);
c@350:     }
c@350: }
c@350: