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