1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */
6 Centre for Digital Music, Queen Mary, University of London.
7 This file 2005-2006 Christian Landone.
9 This program is free software; you can redistribute it and/or
10 modify it under the terms of the GNU General Public License as
11 published by the Free Software Foundation; either version 2 of the
12 License, or (at your option) any later version. See the file
13 COPYING included with this distribution for more information.
16 #include "MathUtilities.h"
25 double MathUtilities::mod(double x, double y)
27 double a = floor( x / y );
29 double b = x - ( y * a );
33 double MathUtilities::princarg(double ang)
37 ValOut = mod( ang + M_PI, -2 * M_PI ) + M_PI;
42 void MathUtilities::getAlphaNorm(const double *data, int len, int alpha, double* ANorm)
48 for( i = 0; i < len; i++)
52 a += ::pow( fabs(temp), double(alpha) );
55 a = ::pow( a, ( 1.0 / (double) alpha ) );
60 double MathUtilities::getAlphaNorm( const vector <double> &data, int alpha )
63 int len = data.size();
67 for( i = 0; i < len; i++)
70 a += ::pow( fabs(temp), double(alpha) );
73 a = ::pow( a, ( 1.0 / (double) alpha ) );
78 double MathUtilities::round(double x)
81 return -floor(-x + 0.5);
83 return floor(x + 0.5);
87 double MathUtilities::median(const double *src, int len)
89 if (len == 0) return 0;
91 vector<double> scratch;
92 for (int i = 0; i < len; ++i) scratch.push_back(src[i]);
93 sort(scratch.begin(), scratch.end());
97 return (scratch[middle] + scratch[middle - 1]) / 2;
99 return scratch[middle];
103 double MathUtilities::sum(const double *src, int len)
108 for( i = 0; i < len; i++)
116 double MathUtilities::mean(const double *src, int len)
120 if (len == 0) return 0;
122 double s = sum( src, len );
124 retVal = s / (double)len;
129 double MathUtilities::mean(const vector<double> &src,
135 if (count == 0) return 0;
137 for (int i = 0; i < (int)count; ++i)
139 sum += src[start + i];
145 void MathUtilities::getFrameMinMax(const double *data, int len, double *min, double *max)
158 for( i = 0; i < len; i++)
174 int MathUtilities::getMax( double* pData, int Length, double* pMax )
180 double max = pData[0];
182 for( i = 0; i < Length; i++)
194 if (pMax) *pMax = max;
200 int MathUtilities::getMax( const vector<double> & data, double* pMax )
206 double max = data[0];
208 for( i = 0; i < int(data.size()); i++)
220 if (pMax) *pMax = max;
226 void MathUtilities::circShift( double* pData, int length, int shift)
228 shift = shift % length;
232 for( i = 0; i < shift; i++)
234 temp=*(pData + length - 1);
236 for( n = length-2; n >= 0; n--)
238 *(pData+n+1)=*(pData+n);
245 int MathUtilities::compareInt (const void * a, const void * b)
247 return ( *(int*)a - *(int*)b );
250 void MathUtilities::normalise(double *data, int length, NormaliseType type)
254 case NormaliseNone: return;
256 case NormaliseUnitSum:
259 for (int i = 0; i < length; ++i) {
263 for (int i = 0; i < length; ++i) {
270 case NormaliseUnitMax:
273 for (int i = 0; i < length; ++i) {
274 if (fabs(data[i]) > max) {
279 for (int i = 0; i < length; ++i) {
289 void MathUtilities::normalise(vector<double> &data, NormaliseType type)
293 case NormaliseNone: return;
295 case NormaliseUnitSum:
298 for (int i = 0; i < (int)data.size(); ++i) sum += data[i];
300 for (int i = 0; i < (int)data.size(); ++i) data[i] /= sum;
305 case NormaliseUnitMax:
308 for (int i = 0; i < (int)data.size(); ++i) {
309 if (fabs(data[i]) > max) max = fabs(data[i]);
312 for (int i = 0; i < (int)data.size(); ++i) data[i] /= max;
320 double MathUtilities::getLpNorm(const vector<double> &data, int p)
323 for (int i = 0; i < int(data.size()); ++i) {
324 tot += abs(pow(data[i], p));
326 return pow(tot, 1.0 / p);
329 vector<double> MathUtilities::normaliseLp(const vector<double> &data,
333 int n = int(data.size());
334 if (n == 0 || p == 0) return data;
335 double norm = getLpNorm(data, p);
336 if (norm < threshold) {
337 return vector<double>(n, 1.0 / pow(n, 1.0 / p)); // unit vector
339 vector<double> out(n);
340 for (int i = 0; i < n; ++i) {
341 out[i] = data[i] / norm;
346 void MathUtilities::adaptiveThreshold(vector<double> &data)
348 int sz = int(data.size());
351 vector<double> smoothed(sz);
356 for (int i = 0; i < sz; ++i) {
358 int first = max(0, i - p_pre);
359 int last = min(sz - 1, i + p_post);
361 smoothed[i] = mean(data, first, last - first + 1);
364 for (int i = 0; i < sz; i++) {
365 data[i] -= smoothed[i];
366 if (data[i] < 0.0) data[i] = 0.0;
371 MathUtilities::isPowerOfTwo(int x)
373 if (x < 1) return false;
374 if (x & (x-1)) return false;
379 MathUtilities::nextPowerOfTwo(int x)
381 if (isPowerOfTwo(x)) return x;
384 while (x) { x >>= 1; n <<= 1; }
389 MathUtilities::previousPowerOfTwo(int x)
391 if (isPowerOfTwo(x)) return x;
395 while (x) { x >>= 1; n <<= 1; }
400 MathUtilities::nearestPowerOfTwo(int x)
402 if (isPowerOfTwo(x)) return x;
403 int n0 = previousPowerOfTwo(x), n1 = nextPowerOfTwo(x);
404 if (x - n0 < n1 - x) return n0;
409 MathUtilities::factorial(int x)
413 for (int i = 1; i <= x; ++i) {
420 MathUtilities::gcd(int a, int b)