c@119: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
c@119: /*
c@119:     Constant-Q library
c@119:     Copyright (c) 2013-2014 Queen Mary, University of London
c@119: 
c@119:     Permission is hereby granted, free of charge, to any person
c@119:     obtaining a copy of this software and associated documentation
c@119:     files (the "Software"), to deal in the Software without
c@119:     restriction, including without limitation the rights to use, copy,
c@119:     modify, merge, publish, distribute, sublicense, and/or sell copies
c@119:     of the Software, and to permit persons to whom the Software is
c@119:     furnished to do so, subject to the following conditions:
c@119: 
c@119:     The above copyright notice and this permission notice shall be
c@119:     included in all copies or substantial portions of the Software.
c@119: 
c@119:     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
c@119:     EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
c@119:     MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
c@119:     NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY
c@119:     CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF
c@119:     CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
c@119:     WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
c@119: 
c@119:     Except as contained in this notice, the names of the Centre for
c@119:     Digital Music; Queen Mary, University of London; and Chris Cannam
c@119:     shall not be used in advertising or otherwise to promote the sale,
c@119:     use or other dealings in this Software without prior written
c@119:     authorization.
c@119: */
c@119: 
c@119: #include "MathUtilities.h"
c@119: 
c@119: #include <iostream>
c@119: #include <algorithm>
c@119: #include <vector>
c@119: #include <cmath>
c@119: 
c@119: 
c@119: double MathUtilities::mod(double x, double y)
c@119: {
c@119:     double a = floor( x / y );
c@119: 
c@119:     double b = x - ( y * a );
c@119:     return b;
c@119: }
c@119: 
c@119: double MathUtilities::princarg(double ang)
c@119: {
c@119:     double ValOut;
c@119: 
c@119:     ValOut = mod( ang + M_PI, -2 * M_PI ) + M_PI;
c@119: 
c@119:     return ValOut;
c@119: }
c@119: 
c@119: void MathUtilities::getAlphaNorm(const double *data, unsigned int len, unsigned int alpha, double* ANorm)
c@119: {
c@119:     unsigned int i;
c@119:     double temp = 0.0;
c@119:     double a=0.0;
c@119: 	
c@119:     for( i = 0; i < len; i++)
c@119:     {
c@119: 	temp = data[ i ];
c@119: 		
c@119: 	a  += ::pow( fabs(temp), double(alpha) );
c@119:     }
c@119:     a /= ( double )len;
c@119:     a = ::pow( a, ( 1.0 / (double) alpha ) );
c@119: 
c@119:     *ANorm = a;
c@119: }
c@119: 
c@119: double MathUtilities::getAlphaNorm( const std::vector <double> &data, unsigned int alpha )
c@119: {
c@119:     unsigned int i;
c@119:     unsigned int len = data.size();
c@119:     double temp = 0.0;
c@119:     double a=0.0;
c@119: 	
c@119:     for( i = 0; i < len; i++)
c@119:     {
c@119: 	temp = data[ i ];
c@119: 		
c@119: 	a  += ::pow( fabs(temp), double(alpha) );
c@119:     }
c@119:     a /= ( double )len;
c@119:     a = ::pow( a, ( 1.0 / (double) alpha ) );
c@119: 
c@119:     return a;
c@119: }
c@119: 
c@119: double MathUtilities::round(double x)
c@119: {
c@119:     if (x < 0) {
c@119:         return -floor(-x + 0.5);
c@119:     } else {
c@119:         return floor(x + 0.5);
c@119:     }
c@119: }
c@119: 
c@119: double MathUtilities::median(const double *src, unsigned int len)
c@119: {
c@119:     if (len == 0) return 0;
c@119:     
c@119:     std::vector<double> scratch;
c@126:     for (int i = 0; i < (int)len; ++i) scratch.push_back(src[i]);
c@119:     std::sort(scratch.begin(), scratch.end());
c@119: 
c@119:     int middle = len/2;
c@119:     if (len % 2 == 0) {
c@119:         return (scratch[middle] + scratch[middle - 1]) / 2;
c@119:     } else {
c@119:         return scratch[middle];
c@119:     }
c@119: }
c@119: 
c@119: double MathUtilities::sum(const double *src, unsigned int len)
c@119: {
c@119:     unsigned int i ;
c@119:     double retVal =0.0;
c@119: 
c@119:     for(  i = 0; i < len; i++)
c@119:     {
c@119: 	retVal += src[ i ];
c@119:     }
c@119: 
c@119:     return retVal;
c@119: }
c@119: 
c@119: double MathUtilities::mean(const double *src, unsigned int len)
c@119: {
c@119:     double retVal =0.0;
c@119: 
c@119:     if (len == 0) return 0;
c@119: 
c@119:     double s = sum( src, len );
c@119: 	
c@119:     retVal =  s  / (double)len;
c@119: 
c@119:     return retVal;
c@119: }
c@119: 
c@119: double MathUtilities::mean(const std::vector<double> &src,
c@119:                            unsigned int start,
c@119:                            unsigned int count)
c@119: {
c@119:     double sum = 0.;
c@119: 	
c@119:     if (count == 0) return 0;
c@119:     
c@119:     for (int i = 0; i < (int)count; ++i)
c@119:     {
c@119:         sum += src[start + i];
c@119:     }
c@119: 
c@119:     return sum / count;
c@119: }
c@119: 
c@119: void MathUtilities::getFrameMinMax(const double *data, unsigned int len, double *min, double *max)
c@119: {
c@119:     unsigned int i;
c@119:     double temp = 0.0;
c@119: 
c@119:     if (len == 0) {
c@119:         *min = *max = 0;
c@119:         return;
c@119:     }
c@119: 	
c@119:     *min = data[0];
c@119:     *max = data[0];
c@119: 
c@119:     for( i = 0; i < len; i++)
c@119:     {
c@119: 	temp = data[ i ];
c@119: 
c@119: 	if( temp < *min )
c@119: 	{
c@119: 	    *min =  temp ;
c@119: 	}
c@119: 	if( temp > *max )
c@119: 	{
c@119: 	    *max =  temp ;
c@119: 	}
c@119: 		
c@119:     }
c@119: }
c@119: 
c@119: int MathUtilities::getMax( double* pData, unsigned int Length, double* pMax )
c@119: {
c@119: 	unsigned int index = 0;
c@119: 	unsigned int i;
c@119: 	double temp = 0.0;
c@119: 	
c@119: 	double max = pData[0];
c@119: 
c@119: 	for( i = 0; i < Length; i++)
c@119: 	{
c@119: 		temp = pData[ i ];
c@119: 
c@119: 		if( temp > max )
c@119: 		{
c@119: 			max =  temp ;
c@119: 			index = i;
c@119: 		}
c@119: 		
c@119:    	}
c@119: 
c@119: 	if (pMax) *pMax = max;
c@119: 
c@119: 
c@119: 	return index;
c@119: }
c@119: 
c@119: int MathUtilities::getMax( const std::vector<double> & data, double* pMax )
c@119: {
c@119: 	unsigned int index = 0;
c@119: 	unsigned int i;
c@119: 	double temp = 0.0;
c@119: 	
c@119: 	double max = data[0];
c@119: 
c@119: 	for( i = 0; i < data.size(); i++)
c@119: 	{
c@119: 		temp = data[ i ];
c@119: 
c@119: 		if( temp > max )
c@119: 		{
c@119: 			max =  temp ;
c@119: 			index = i;
c@119: 		}
c@119: 		
c@119:    	}
c@119: 
c@119: 	if (pMax) *pMax = max;
c@119: 
c@119: 
c@119: 	return index;
c@119: }
c@119: 
c@119: void MathUtilities::circShift( double* pData, int length, int shift)
c@119: {
c@119: 	shift = shift % length;
c@119: 	double temp;
c@119: 	int i,n;
c@119: 
c@119: 	for( i = 0; i < shift; i++)
c@119: 	{
c@119: 		temp=*(pData + length - 1);
c@119: 
c@119: 		for( n = length-2; n >= 0; n--)
c@119: 		{
c@119: 			*(pData+n+1)=*(pData+n);
c@119: 		}
c@119: 
c@119:         *pData = temp;
c@119:     }
c@119: }
c@119: 
c@119: int MathUtilities::compareInt (const void * a, const void * b)
c@119: {
c@119:   return ( *(int*)a - *(int*)b );
c@119: }
c@119: 
c@119: void MathUtilities::normalise(double *data, int length, NormaliseType type)
c@119: {
c@119:     switch (type) {
c@119: 
c@119:     case NormaliseNone: return;
c@119: 
c@119:     case NormaliseUnitSum:
c@119:     {
c@119:         double sum = 0.0;
c@119:         for (int i = 0; i < length; ++i) {
c@119:             sum += data[i];
c@119:         }
c@119:         if (sum != 0.0) {
c@119:             for (int i = 0; i < length; ++i) {
c@119:                 data[i] /= sum;
c@119:             }
c@119:         }
c@119:     }
c@119:     break;
c@119: 
c@119:     case NormaliseUnitMax:
c@119:     {
c@119:         double max = 0.0;
c@119:         for (int i = 0; i < length; ++i) {
c@119:             if (fabs(data[i]) > max) {
c@119:                 max = fabs(data[i]);
c@119:             }
c@119:         }
c@119:         if (max != 0.0) {
c@119:             for (int i = 0; i < length; ++i) {
c@119:                 data[i] /= max;
c@119:             }
c@119:         }
c@119:     }
c@119:     break;
c@119: 
c@119:     }
c@119: }
c@119: 
c@119: void MathUtilities::normalise(std::vector<double> &data, NormaliseType type)
c@119: {
c@119:     switch (type) {
c@119: 
c@119:     case NormaliseNone: return;
c@119: 
c@119:     case NormaliseUnitSum:
c@119:     {
c@119:         double sum = 0.0;
c@119:         for (int i = 0; i < (int)data.size(); ++i) sum += data[i];
c@119:         if (sum != 0.0) {
c@119:             for (int i = 0; i < (int)data.size(); ++i) data[i] /= sum;
c@119:         }
c@119:     }
c@119:     break;
c@119: 
c@119:     case NormaliseUnitMax:
c@119:     {
c@119:         double max = 0.0;
c@119:         for (int i = 0; i < (int)data.size(); ++i) {
c@119:             if (fabs(data[i]) > max) max = fabs(data[i]);
c@119:         }
c@119:         if (max != 0.0) {
c@119:             for (int i = 0; i < (int)data.size(); ++i) data[i] /= max;
c@119:         }
c@119:     }
c@119:     break;
c@119: 
c@119:     }
c@119: }
c@119: 
c@119: void MathUtilities::adaptiveThreshold(std::vector<double> &data)
c@119: {
c@119:     int sz = int(data.size());
c@119:     if (sz == 0) return;
c@119: 
c@119:     std::vector<double> smoothed(sz);
c@119: 	
c@119:     int p_pre = 8;
c@119:     int p_post = 7;
c@119: 
c@119:     for (int i = 0; i < sz; ++i) {
c@119: 
c@119:         int first = std::max(0,      i - p_pre);
c@119:         int last  = std::min(sz - 1, i + p_post);
c@119: 
c@119:         smoothed[i] = mean(data, first, last - first + 1);
c@119:     }
c@119: 
c@119:     for (int i = 0; i < sz; i++) {
c@119:         data[i] -= smoothed[i];
c@119:         if (data[i] < 0.0) data[i] = 0.0;
c@119:     }
c@119: }
c@119: 
c@119: bool
c@119: MathUtilities::isPowerOfTwo(int x)
c@119: {
c@119:     if (x < 1) return false;
c@119:     if (x & (x-1)) return false;
c@119:     return true;
c@119: }
c@119: 
c@119: int
c@119: MathUtilities::nextPowerOfTwo(int x)
c@119: {
c@119:     if (isPowerOfTwo(x)) return x;
c@119:     if (x < 1) return 1;
c@119:     int n = 1;
c@119:     while (x) { x >>= 1; n <<= 1; }
c@119:     return n;
c@119: }
c@119: 
c@119: int
c@119: MathUtilities::previousPowerOfTwo(int x)
c@119: {
c@119:     if (isPowerOfTwo(x)) return x;
c@119:     if (x < 1) return 1;
c@119:     int n = 1;
c@119:     x >>= 1;
c@119:     while (x) { x >>= 1; n <<= 1; }
c@119:     return n;
c@119: }
c@119: 
c@119: int
c@119: MathUtilities::nearestPowerOfTwo(int x)
c@119: {
c@119:     if (isPowerOfTwo(x)) return x;
c@119:     int n0 = previousPowerOfTwo(x), n1 = nextPowerOfTwo(x);
c@119:     if (x - n0 < n1 - x) return n0;
c@119:     else return n1;
c@119: }
c@119: 
c@119: double
c@119: MathUtilities::factorial(int x)
c@119: {
c@119:     if (x < 0) return 0;
c@119:     double f = 1;
c@119:     for (int i = 1; i <= x; ++i) {
c@119: 	f = f * i;
c@119:     }
c@119:     return f;
c@119: }
c@119: 
c@119: int
c@119: MathUtilities::gcd(int a, int b)
c@119: {
c@119:     int c = a % b;
c@119:     if (c == 0) {
c@119:         return b;
c@119:     } else {
c@119:         return gcd(b, c);
c@119:     }
c@119: }
c@119: