cannam@0: /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
cannam@0: 
cannam@0: /*
cannam@0:     QM DSP Library
cannam@0: 
cannam@0:     Centre for Digital Music, Queen Mary, University of London.
Chris@84:     This file 2005-2006 Christian Landone.
Chris@84: 
Chris@84:     This program is free software; you can redistribute it and/or
Chris@84:     modify it under the terms of the GNU General Public License as
Chris@84:     published by the Free Software Foundation; either version 2 of the
Chris@84:     License, or (at your option) any later version.  See the file
Chris@84:     COPYING included with this distribution for more information.
cannam@0: */
cannam@0: 
cannam@0: #include "MathUtilities.h"
cannam@0: 
cannam@0: #include <iostream>
cannam@0: #include <cmath>
cannam@0: 
cannam@0: 
cannam@0: double MathUtilities::mod(double x, double y)
cannam@0: {
cannam@0:     double a = floor( x / y );
cannam@0: 
cannam@0:     double b = x - ( y * a );
cannam@0:     return b;
cannam@0: }
cannam@0: 
cannam@0: double MathUtilities::princarg(double ang)
cannam@0: {
cannam@0:     double ValOut;
cannam@0: 
cannam@0:     ValOut = mod( ang + M_PI, -2 * M_PI ) + M_PI;
cannam@0: 
cannam@0:     return ValOut;
cannam@0: }
cannam@0: 
cannam@0: void MathUtilities::getAlphaNorm(const double *data, unsigned int len, unsigned int alpha, double* ANorm)
cannam@0: {
cannam@0:     unsigned int i;
cannam@0:     double temp = 0.0;
cannam@0:     double a=0.0;
cannam@0: 	
cannam@0:     for( i = 0; i < len; i++)
cannam@0:     {
cannam@0: 	temp = data[ i ];
cannam@0: 		
cannam@79: 	a  += ::pow( fabs(temp), double(alpha) );
cannam@0:     }
cannam@0:     a /= ( double )len;
cannam@0:     a = ::pow( a, ( 1.0 / (double) alpha ) );
cannam@0: 
cannam@0:     *ANorm = a;
cannam@0: }
cannam@0: 
cannam@0: double MathUtilities::getAlphaNorm( const std::vector <double> &data, unsigned int alpha )
cannam@0: {
cannam@0:     unsigned int i;
cannam@0:     unsigned int len = data.size();
cannam@0:     double temp = 0.0;
cannam@0:     double a=0.0;
cannam@0: 	
cannam@0:     for( i = 0; i < len; i++)
cannam@0:     {
cannam@0: 	temp = data[ i ];
cannam@0: 		
cannam@79: 	a  += ::pow( fabs(temp), double(alpha) );
cannam@0:     }
cannam@0:     a /= ( double )len;
cannam@0:     a = ::pow( a, ( 1.0 / (double) alpha ) );
cannam@0: 
cannam@0:     return a;
cannam@0: }
cannam@0: 
cannam@0: double MathUtilities::round(double x)
cannam@0: {
cannam@0:     double val = (double)floor(x + 0.5);
cannam@0:   
cannam@0:     return val;
cannam@0: }
cannam@0: 
cannam@0: double MathUtilities::median(const double *src, unsigned int len)
cannam@0: {
cannam@0:     unsigned int i, j;
cannam@0:     double tmp = 0.0;
cannam@0:     double tempMedian;
cannam@0:     double medianVal;
cannam@0:  
cannam@0:     double* scratch = new double[ len ];//Vector < double > sortedX = Vector < double > ( size );
cannam@0: 
cannam@0:     for ( i = 0; i < len; i++ )
cannam@0:     {
cannam@0: 	scratch[i] = src[i];
cannam@0:     }
cannam@0: 
cannam@0:     for ( i = 0; i < len - 1; i++ )
cannam@0:     {
cannam@0: 	for ( j = 0; j < len - 1 - i; j++ )
cannam@0: 	{
cannam@0: 	    if ( scratch[j + 1] < scratch[j] )
cannam@0: 	    {
cannam@0: 		// compare the two neighbors
cannam@0: 		tmp = scratch[j]; // swap a[j] and a[j+1]
cannam@0: 		scratch[j] = scratch[j + 1];
cannam@0: 		scratch[j + 1] = tmp;
cannam@0: 	    }
cannam@0: 	}
cannam@0:     }
cannam@0:     int middle;
cannam@0:     if ( len % 2 == 0 )
cannam@0:     {
cannam@0: 	middle = len / 2;
cannam@0: 	tempMedian = ( scratch[middle] + scratch[middle - 1] ) / 2;
cannam@0:     }
cannam@0:     else
cannam@0:     {
cannam@0: 	middle = ( int )floor( len / 2.0 );
cannam@0: 	tempMedian = scratch[middle];
cannam@0:     }
cannam@0: 
cannam@0:     medianVal = tempMedian;
cannam@0: 
cannam@0:     delete [] scratch;
cannam@0:     return medianVal;
cannam@0: }
cannam@0: 
cannam@0: double MathUtilities::sum(const double *src, unsigned int len)
cannam@0: {
cannam@0:     unsigned int i ;
cannam@0:     double retVal =0.0;
cannam@0: 
cannam@0:     for(  i = 0; i < len; i++)
cannam@0:     {
cannam@0: 	retVal += src[ i ];
cannam@0:     }
cannam@0: 
cannam@0:     return retVal;
cannam@0: }
cannam@0: 
cannam@0: double MathUtilities::mean(const double *src, unsigned int len)
cannam@0: {
cannam@0:     double retVal =0.0;
cannam@0: 
cannam@0:     double s = sum( src, len );
cannam@0: 	
cannam@0:     retVal =  s  / (double)len;
cannam@0: 
cannam@0:     return retVal;
cannam@0: }
cannam@0: 
cannam@54: double MathUtilities::mean(const std::vector<double> &src,
cannam@54:                            unsigned int start,
cannam@54:                            unsigned int count)
cannam@54: {
cannam@54:     double sum = 0.;
cannam@54: 	
Chris@102:     for (int i = 0; i < (int)count; ++i)
cannam@54:     {
cannam@54:         sum += src[start + i];
cannam@54:     }
cannam@54: 
cannam@54:     return sum / count;
cannam@54: }
cannam@54: 
cannam@0: void MathUtilities::getFrameMinMax(const double *data, unsigned int len, double *min, double *max)
cannam@0: {
cannam@0:     unsigned int i;
cannam@0:     double temp = 0.0;
cannam@58: 
cannam@58:     if (len == 0) {
cannam@58:         *min = *max = 0;
cannam@58:         return;
cannam@58:     }
cannam@0: 	
cannam@0:     *min = data[0];
cannam@0:     *max = data[0];
cannam@0: 
cannam@0:     for( i = 0; i < len; i++)
cannam@0:     {
cannam@0: 	temp = data[ i ];
cannam@0: 
cannam@0: 	if( temp < *min )
cannam@0: 	{
cannam@0: 	    *min =  temp ;
cannam@0: 	}
cannam@0: 	if( temp > *max )
cannam@0: 	{
cannam@0: 	    *max =  temp ;
cannam@0: 	}
cannam@0: 		
cannam@0:     }
cannam@0: }
cannam@7: 
cannam@7: int MathUtilities::getMax( double* pData, unsigned int Length, double* pMax )
cannam@7: {
cannam@7: 	unsigned int index = 0;
cannam@7: 	unsigned int i;
cannam@7: 	double temp = 0.0;
cannam@7: 	
cannam@7: 	double max = pData[0];
cannam@7: 
cannam@7: 	for( i = 0; i < Length; i++)
cannam@7: 	{
cannam@7: 		temp = pData[ i ];
cannam@7: 
cannam@7: 		if( temp > max )
cannam@7: 		{
cannam@7: 			max =  temp ;
cannam@7: 			index = i;
cannam@7: 		}
cannam@7: 		
cannam@7:    	}
cannam@7: 
cannam@54: 	if (pMax) *pMax = max;
cannam@54: 
cannam@54: 
cannam@54: 	return index;
cannam@54: }
cannam@54: 
cannam@54: int MathUtilities::getMax( const std::vector<double> & data, double* pMax )
cannam@54: {
cannam@54: 	unsigned int index = 0;
cannam@54: 	unsigned int i;
cannam@54: 	double temp = 0.0;
cannam@54: 	
cannam@54: 	double max = data[0];
cannam@54: 
cannam@54: 	for( i = 0; i < data.size(); i++)
cannam@54: 	{
cannam@54: 		temp = data[ i ];
cannam@54: 
cannam@54: 		if( temp > max )
cannam@54: 		{
cannam@54: 			max =  temp ;
cannam@54: 			index = i;
cannam@54: 		}
cannam@54: 		
cannam@54:    	}
cannam@54: 
cannam@54: 	if (pMax) *pMax = max;
cannam@7: 
cannam@7: 
cannam@7: 	return index;
cannam@7: }
cannam@7: 
cannam@7: void MathUtilities::circShift( double* pData, int length, int shift)
cannam@7: {
cannam@7: 	shift = shift % length;
cannam@7: 	double temp;
cannam@7: 	int i,n;
cannam@7: 
cannam@7: 	for( i = 0; i < shift; i++)
cannam@7: 	{
cannam@7: 		temp=*(pData + length - 1);
cannam@7: 
cannam@7: 		for( n = length-2; n >= 0; n--)
cannam@7: 		{
cannam@7: 			*(pData+n+1)=*(pData+n);
cannam@7: 		}
cannam@7: 
cannam@7:         *pData = temp;
cannam@7:     }
cannam@7: }
cannam@7: 
cannam@7: int MathUtilities::compareInt (const void * a, const void * b)
cannam@7: {
cannam@7:   return ( *(int*)a - *(int*)b );
cannam@7: }
cannam@7: 
cannam@34: void MathUtilities::normalise(double *data, int length, NormaliseType type)
cannam@34: {
cannam@34:     switch (type) {
cannam@34: 
cannam@34:     case NormaliseNone: return;
cannam@34: 
cannam@34:     case NormaliseUnitSum:
cannam@34:     {
cannam@34:         double sum = 0.0;
cannam@34:         for (int i = 0; i < length; ++i) {
cannam@34:             sum += data[i];
cannam@34:         }
cannam@34:         if (sum != 0.0) {
cannam@34:             for (int i = 0; i < length; ++i) {
cannam@34:                 data[i] /= sum;
cannam@34:             }
cannam@34:         }
cannam@34:     }
cannam@34:     break;
cannam@34: 
cannam@34:     case NormaliseUnitMax:
cannam@34:     {
cannam@34:         double max = 0.0;
cannam@34:         for (int i = 0; i < length; ++i) {
cannam@34:             if (fabs(data[i]) > max) {
cannam@34:                 max = fabs(data[i]);
cannam@34:             }
cannam@34:         }
cannam@34:         if (max != 0.0) {
cannam@34:             for (int i = 0; i < length; ++i) {
cannam@34:                 data[i] /= max;
cannam@34:             }
cannam@34:         }
cannam@34:     }
cannam@34:     break;
cannam@34: 
cannam@34:     }
cannam@34: }
cannam@34: 
cannam@34: void MathUtilities::normalise(std::vector<double> &data, NormaliseType type)
cannam@34: {
cannam@34:     switch (type) {
cannam@34: 
cannam@34:     case NormaliseNone: return;
cannam@34: 
cannam@34:     case NormaliseUnitSum:
cannam@34:     {
cannam@34:         double sum = 0.0;
Chris@102:         for (int i = 0; i < (int)data.size(); ++i) sum += data[i];
cannam@34:         if (sum != 0.0) {
Chris@102:             for (int i = 0; i < (int)data.size(); ++i) data[i] /= sum;
cannam@34:         }
cannam@34:     }
cannam@34:     break;
cannam@34: 
cannam@34:     case NormaliseUnitMax:
cannam@34:     {
cannam@34:         double max = 0.0;
Chris@102:         for (int i = 0; i < (int)data.size(); ++i) {
cannam@34:             if (fabs(data[i]) > max) max = fabs(data[i]);
cannam@34:         }
cannam@34:         if (max != 0.0) {
Chris@102:             for (int i = 0; i < (int)data.size(); ++i) data[i] /= max;
cannam@34:         }
cannam@34:     }
cannam@34:     break;
cannam@34: 
cannam@34:     }
cannam@34: }
cannam@34: 
cannam@54: void MathUtilities::adaptiveThreshold(std::vector<double> &data)
cannam@54: {
cannam@54:     int sz = int(data.size());
cannam@54:     if (sz == 0) return;
cannam@54: 
cannam@54:     std::vector<double> smoothed(sz);
cannam@54: 	
cannam@54:     int p_pre = 8;
cannam@54:     int p_post = 7;
cannam@54: 
cannam@54:     for (int i = 0; i < sz; ++i) {
cannam@54: 
cannam@54:         int first = std::max(0,      i - p_pre);
cannam@54:         int last  = std::min(sz - 1, i + p_post);
cannam@54: 
cannam@54:         smoothed[i] = mean(data, first, last - first + 1);
cannam@54:     }
cannam@54: 
cannam@54:     for (int i = 0; i < sz; i++) {
cannam@54:         data[i] -= smoothed[i];
cannam@54:         if (data[i] < 0.0) data[i] = 0.0;
cannam@54:     }
cannam@54: }
cannam@34: 
cannam@55: bool
cannam@55: MathUtilities::isPowerOfTwo(int x)
cannam@55: {
cannam@55:     if (x < 2) return false;
cannam@55:     if (x & (x-1)) return false;
cannam@55:     return true;
cannam@55: }
cannam@55: 
cannam@55: int
cannam@55: MathUtilities::nextPowerOfTwo(int x)
cannam@55: {
cannam@55:     if (isPowerOfTwo(x)) return x;
cannam@55:     int n = 1;
cannam@55:     while (x) { x >>= 1; n <<= 1; }
cannam@55:     return n;
cannam@55: }
cannam@55: 
cannam@55: int
cannam@55: MathUtilities::previousPowerOfTwo(int x)
cannam@55: {
cannam@55:     if (isPowerOfTwo(x)) return x;
cannam@55:     int n = 1;
cannam@55:     x >>= 1;
cannam@55:     while (x) { x >>= 1; n <<= 1; }
cannam@55:     return n;
cannam@55: }
cannam@55: 
cannam@55: int
cannam@55: MathUtilities::nearestPowerOfTwo(int x)
cannam@55: {
cannam@55:     if (isPowerOfTwo(x)) return x;
cannam@55:     int n0 = previousPowerOfTwo(x), n1 = nearestPowerOfTwo(x);
cannam@55:     if (x - n0 < n1 - x) return n0;
cannam@55:     else return n1;
cannam@55: }
cannam@55: