qm-dsp
1.8
|
00001 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */ 00002 00003 /* 00004 QM DSP Library 00005 00006 Centre for Digital Music, Queen Mary, University of London. 00007 This file 2005-2006 Christian Landone. 00008 00009 This program is free software; you can redistribute it and/or 00010 modify it under the terms of the GNU General Public License as 00011 published by the Free Software Foundation; either version 2 of the 00012 License, or (at your option) any later version. See the file 00013 COPYING included with this distribution for more information. 00014 */ 00015 00016 #include "MathUtilities.h" 00017 00018 #include <iostream> 00019 #include <algorithm> 00020 #include <vector> 00021 #include <cmath> 00022 00023 00024 double MathUtilities::mod(double x, double y) 00025 { 00026 double a = floor( x / y ); 00027 00028 double b = x - ( y * a ); 00029 return b; 00030 } 00031 00032 double MathUtilities::princarg(double ang) 00033 { 00034 double ValOut; 00035 00036 ValOut = mod( ang + M_PI, -2 * M_PI ) + M_PI; 00037 00038 return ValOut; 00039 } 00040 00041 void MathUtilities::getAlphaNorm(const double *data, unsigned int len, unsigned int alpha, double* ANorm) 00042 { 00043 unsigned int i; 00044 double temp = 0.0; 00045 double a=0.0; 00046 00047 for( i = 0; i < len; i++) 00048 { 00049 temp = data[ i ]; 00050 00051 a += ::pow( fabs(temp), double(alpha) ); 00052 } 00053 a /= ( double )len; 00054 a = ::pow( a, ( 1.0 / (double) alpha ) ); 00055 00056 *ANorm = a; 00057 } 00058 00059 double MathUtilities::getAlphaNorm( const std::vector <double> &data, unsigned int alpha ) 00060 { 00061 unsigned int i; 00062 unsigned int len = data.size(); 00063 double temp = 0.0; 00064 double a=0.0; 00065 00066 for( i = 0; i < len; i++) 00067 { 00068 temp = data[ i ]; 00069 00070 a += ::pow( fabs(temp), double(alpha) ); 00071 } 00072 a /= ( double )len; 00073 a = ::pow( a, ( 1.0 / (double) alpha ) ); 00074 00075 return a; 00076 } 00077 00078 double MathUtilities::round(double x) 00079 { 00080 if (x < 0) { 00081 return -floor(-x + 0.5); 00082 } else { 00083 return floor(x + 0.5); 00084 } 00085 } 00086 00087 double MathUtilities::median(const double *src, unsigned int len) 00088 { 00089 if (len == 0) return 0; 00090 00091 std::vector<double> scratch; 00092 for (int i = 0; i < len; ++i) scratch.push_back(src[i]); 00093 std::sort(scratch.begin(), scratch.end()); 00094 00095 int middle = len/2; 00096 if (len % 2 == 0) { 00097 return (scratch[middle] + scratch[middle - 1]) / 2; 00098 } else { 00099 return scratch[middle]; 00100 } 00101 } 00102 00103 double MathUtilities::sum(const double *src, unsigned int len) 00104 { 00105 unsigned int i ; 00106 double retVal =0.0; 00107 00108 for( i = 0; i < len; i++) 00109 { 00110 retVal += src[ i ]; 00111 } 00112 00113 return retVal; 00114 } 00115 00116 double MathUtilities::mean(const double *src, unsigned int len) 00117 { 00118 double retVal =0.0; 00119 00120 if (len == 0) return 0; 00121 00122 double s = sum( src, len ); 00123 00124 retVal = s / (double)len; 00125 00126 return retVal; 00127 } 00128 00129 double MathUtilities::mean(const std::vector<double> &src, 00130 unsigned int start, 00131 unsigned int count) 00132 { 00133 double sum = 0.; 00134 00135 if (count == 0) return 0; 00136 00137 for (int i = 0; i < (int)count; ++i) 00138 { 00139 sum += src[start + i]; 00140 } 00141 00142 return sum / count; 00143 } 00144 00145 void MathUtilities::getFrameMinMax(const double *data, unsigned int len, double *min, double *max) 00146 { 00147 unsigned int i; 00148 double temp = 0.0; 00149 00150 if (len == 0) { 00151 *min = *max = 0; 00152 return; 00153 } 00154 00155 *min = data[0]; 00156 *max = data[0]; 00157 00158 for( i = 0; i < len; i++) 00159 { 00160 temp = data[ i ]; 00161 00162 if( temp < *min ) 00163 { 00164 *min = temp ; 00165 } 00166 if( temp > *max ) 00167 { 00168 *max = temp ; 00169 } 00170 00171 } 00172 } 00173 00174 int MathUtilities::getMax( double* pData, unsigned int Length, double* pMax ) 00175 { 00176 unsigned int index = 0; 00177 unsigned int i; 00178 double temp = 0.0; 00179 00180 double max = pData[0]; 00181 00182 for( i = 0; i < Length; i++) 00183 { 00184 temp = pData[ i ]; 00185 00186 if( temp > max ) 00187 { 00188 max = temp ; 00189 index = i; 00190 } 00191 00192 } 00193 00194 if (pMax) *pMax = max; 00195 00196 00197 return index; 00198 } 00199 00200 int MathUtilities::getMax( const std::vector<double> & data, double* pMax ) 00201 { 00202 unsigned int index = 0; 00203 unsigned int i; 00204 double temp = 0.0; 00205 00206 double max = data[0]; 00207 00208 for( i = 0; i < data.size(); i++) 00209 { 00210 temp = data[ i ]; 00211 00212 if( temp > max ) 00213 { 00214 max = temp ; 00215 index = i; 00216 } 00217 00218 } 00219 00220 if (pMax) *pMax = max; 00221 00222 00223 return index; 00224 } 00225 00226 void MathUtilities::circShift( double* pData, int length, int shift) 00227 { 00228 shift = shift % length; 00229 double temp; 00230 int i,n; 00231 00232 for( i = 0; i < shift; i++) 00233 { 00234 temp=*(pData + length - 1); 00235 00236 for( n = length-2; n >= 0; n--) 00237 { 00238 *(pData+n+1)=*(pData+n); 00239 } 00240 00241 *pData = temp; 00242 } 00243 } 00244 00245 int MathUtilities::compareInt (const void * a, const void * b) 00246 { 00247 return ( *(int*)a - *(int*)b ); 00248 } 00249 00250 void MathUtilities::normalise(double *data, int length, NormaliseType type) 00251 { 00252 switch (type) { 00253 00254 case NormaliseNone: return; 00255 00256 case NormaliseUnitSum: 00257 { 00258 double sum = 0.0; 00259 for (int i = 0; i < length; ++i) { 00260 sum += data[i]; 00261 } 00262 if (sum != 0.0) { 00263 for (int i = 0; i < length; ++i) { 00264 data[i] /= sum; 00265 } 00266 } 00267 } 00268 break; 00269 00270 case NormaliseUnitMax: 00271 { 00272 double max = 0.0; 00273 for (int i = 0; i < length; ++i) { 00274 if (fabs(data[i]) > max) { 00275 max = fabs(data[i]); 00276 } 00277 } 00278 if (max != 0.0) { 00279 for (int i = 0; i < length; ++i) { 00280 data[i] /= max; 00281 } 00282 } 00283 } 00284 break; 00285 00286 } 00287 } 00288 00289 void MathUtilities::normalise(std::vector<double> &data, NormaliseType type) 00290 { 00291 switch (type) { 00292 00293 case NormaliseNone: return; 00294 00295 case NormaliseUnitSum: 00296 { 00297 double sum = 0.0; 00298 for (int i = 0; i < (int)data.size(); ++i) sum += data[i]; 00299 if (sum != 0.0) { 00300 for (int i = 0; i < (int)data.size(); ++i) data[i] /= sum; 00301 } 00302 } 00303 break; 00304 00305 case NormaliseUnitMax: 00306 { 00307 double max = 0.0; 00308 for (int i = 0; i < (int)data.size(); ++i) { 00309 if (fabs(data[i]) > max) max = fabs(data[i]); 00310 } 00311 if (max != 0.0) { 00312 for (int i = 0; i < (int)data.size(); ++i) data[i] /= max; 00313 } 00314 } 00315 break; 00316 00317 } 00318 } 00319 00320 void MathUtilities::adaptiveThreshold(std::vector<double> &data) 00321 { 00322 int sz = int(data.size()); 00323 if (sz == 0) return; 00324 00325 std::vector<double> smoothed(sz); 00326 00327 int p_pre = 8; 00328 int p_post = 7; 00329 00330 for (int i = 0; i < sz; ++i) { 00331 00332 int first = std::max(0, i - p_pre); 00333 int last = std::min(sz - 1, i + p_post); 00334 00335 smoothed[i] = mean(data, first, last - first + 1); 00336 } 00337 00338 for (int i = 0; i < sz; i++) { 00339 data[i] -= smoothed[i]; 00340 if (data[i] < 0.0) data[i] = 0.0; 00341 } 00342 } 00343 00344 bool 00345 MathUtilities::isPowerOfTwo(int x) 00346 { 00347 if (x < 1) return false; 00348 if (x & (x-1)) return false; 00349 return true; 00350 } 00351 00352 int 00353 MathUtilities::nextPowerOfTwo(int x) 00354 { 00355 if (isPowerOfTwo(x)) return x; 00356 if (x < 1) return 1; 00357 int n = 1; 00358 while (x) { x >>= 1; n <<= 1; } 00359 return n; 00360 } 00361 00362 int 00363 MathUtilities::previousPowerOfTwo(int x) 00364 { 00365 if (isPowerOfTwo(x)) return x; 00366 if (x < 1) return 1; 00367 int n = 1; 00368 x >>= 1; 00369 while (x) { x >>= 1; n <<= 1; } 00370 return n; 00371 } 00372 00373 int 00374 MathUtilities::nearestPowerOfTwo(int x) 00375 { 00376 if (isPowerOfTwo(x)) return x; 00377 int n0 = previousPowerOfTwo(x), n1 = nextPowerOfTwo(x); 00378 if (x - n0 < n1 - x) return n0; 00379 else return n1; 00380 } 00381 00382 double 00383 MathUtilities::factorial(int x) 00384 { 00385 if (x < 0) return 0; 00386 double f = 1; 00387 for (int i = 1; i <= x; ++i) { 00388 f = f * i; 00389 } 00390 return f; 00391 } 00392 00393 int 00394 MathUtilities::gcd(int a, int b) 00395 { 00396 int c = a % b; 00397 if (c == 0) { 00398 return b; 00399 } else { 00400 return gcd(b, c); 00401 } 00402 } 00403