annotate maths/MathUtilities.cpp @ 325:31f22daeba64

Fix compiler warnings
author Chris Cannam <c.cannam@qmul.ac.uk>
date Thu, 13 Jun 2013 10:23:09 +0100
parents d5014ab8b0e5
children 50fae18236ee
rev   line source
c@241 1 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*- vi:set ts=8 sts=4 sw=4: */
c@241 2
c@241 3 /*
c@241 4 QM DSP Library
c@241 5
c@241 6 Centre for Digital Music, Queen Mary, University of London.
c@309 7 This file 2005-2006 Christian Landone.
c@309 8
c@309 9 This program is free software; you can redistribute it and/or
c@309 10 modify it under the terms of the GNU General Public License as
c@309 11 published by the Free Software Foundation; either version 2 of the
c@309 12 License, or (at your option) any later version. See the file
c@309 13 COPYING included with this distribution for more information.
c@241 14 */
c@241 15
c@241 16 #include "MathUtilities.h"
c@241 17
c@241 18 #include <iostream>
c@241 19 #include <cmath>
c@241 20
c@241 21
c@241 22 double MathUtilities::mod(double x, double y)
c@241 23 {
c@241 24 double a = floor( x / y );
c@241 25
c@241 26 double b = x - ( y * a );
c@241 27 return b;
c@241 28 }
c@241 29
c@241 30 double MathUtilities::princarg(double ang)
c@241 31 {
c@241 32 double ValOut;
c@241 33
c@241 34 ValOut = mod( ang + M_PI, -2 * M_PI ) + M_PI;
c@241 35
c@241 36 return ValOut;
c@241 37 }
c@241 38
c@241 39 void MathUtilities::getAlphaNorm(const double *data, unsigned int len, unsigned int alpha, double* ANorm)
c@241 40 {
c@241 41 unsigned int i;
c@241 42 double temp = 0.0;
c@241 43 double a=0.0;
c@241 44
c@241 45 for( i = 0; i < len; i++)
c@241 46 {
c@241 47 temp = data[ i ];
c@241 48
c@304 49 a += ::pow( fabs(temp), double(alpha) );
c@241 50 }
c@241 51 a /= ( double )len;
c@241 52 a = ::pow( a, ( 1.0 / (double) alpha ) );
c@241 53
c@241 54 *ANorm = a;
c@241 55 }
c@241 56
c@241 57 double MathUtilities::getAlphaNorm( const std::vector <double> &data, unsigned int alpha )
c@241 58 {
c@241 59 unsigned int i;
c@241 60 unsigned int len = data.size();
c@241 61 double temp = 0.0;
c@241 62 double a=0.0;
c@241 63
c@241 64 for( i = 0; i < len; i++)
c@241 65 {
c@241 66 temp = data[ i ];
c@241 67
c@304 68 a += ::pow( fabs(temp), double(alpha) );
c@241 69 }
c@241 70 a /= ( double )len;
c@241 71 a = ::pow( a, ( 1.0 / (double) alpha ) );
c@241 72
c@241 73 return a;
c@241 74 }
c@241 75
c@241 76 double MathUtilities::round(double x)
c@241 77 {
c@241 78 double val = (double)floor(x + 0.5);
c@241 79
c@241 80 return val;
c@241 81 }
c@241 82
c@241 83 double MathUtilities::median(const double *src, unsigned int len)
c@241 84 {
c@241 85 unsigned int i, j;
c@241 86 double tmp = 0.0;
c@241 87 double tempMedian;
c@241 88 double medianVal;
c@241 89
c@241 90 double* scratch = new double[ len ];//Vector < double > sortedX = Vector < double > ( size );
c@241 91
c@241 92 for ( i = 0; i < len; i++ )
c@241 93 {
c@241 94 scratch[i] = src[i];
c@241 95 }
c@241 96
c@241 97 for ( i = 0; i < len - 1; i++ )
c@241 98 {
c@241 99 for ( j = 0; j < len - 1 - i; j++ )
c@241 100 {
c@241 101 if ( scratch[j + 1] < scratch[j] )
c@241 102 {
c@241 103 // compare the two neighbors
c@241 104 tmp = scratch[j]; // swap a[j] and a[j+1]
c@241 105 scratch[j] = scratch[j + 1];
c@241 106 scratch[j + 1] = tmp;
c@241 107 }
c@241 108 }
c@241 109 }
c@241 110 int middle;
c@241 111 if ( len % 2 == 0 )
c@241 112 {
c@241 113 middle = len / 2;
c@241 114 tempMedian = ( scratch[middle] + scratch[middle - 1] ) / 2;
c@241 115 }
c@241 116 else
c@241 117 {
c@241 118 middle = ( int )floor( len / 2.0 );
c@241 119 tempMedian = scratch[middle];
c@241 120 }
c@241 121
c@241 122 medianVal = tempMedian;
c@241 123
c@241 124 delete [] scratch;
c@241 125 return medianVal;
c@241 126 }
c@241 127
c@241 128 double MathUtilities::sum(const double *src, unsigned int len)
c@241 129 {
c@241 130 unsigned int i ;
c@241 131 double retVal =0.0;
c@241 132
c@241 133 for( i = 0; i < len; i++)
c@241 134 {
c@241 135 retVal += src[ i ];
c@241 136 }
c@241 137
c@241 138 return retVal;
c@241 139 }
c@241 140
c@241 141 double MathUtilities::mean(const double *src, unsigned int len)
c@241 142 {
c@241 143 double retVal =0.0;
c@241 144
c@241 145 double s = sum( src, len );
c@241 146
c@241 147 retVal = s / (double)len;
c@241 148
c@241 149 return retVal;
c@241 150 }
c@241 151
c@279 152 double MathUtilities::mean(const std::vector<double> &src,
c@279 153 unsigned int start,
c@279 154 unsigned int count)
c@279 155 {
c@279 156 double sum = 0.;
c@279 157
c@325 158 for (int i = 0; i < (int)count; ++i)
c@279 159 {
c@279 160 sum += src[start + i];
c@279 161 }
c@279 162
c@279 163 return sum / count;
c@279 164 }
c@279 165
c@241 166 void MathUtilities::getFrameMinMax(const double *data, unsigned int len, double *min, double *max)
c@241 167 {
c@241 168 unsigned int i;
c@241 169 double temp = 0.0;
c@283 170
c@283 171 if (len == 0) {
c@283 172 *min = *max = 0;
c@283 173 return;
c@283 174 }
c@241 175
c@241 176 *min = data[0];
c@241 177 *max = data[0];
c@241 178
c@241 179 for( i = 0; i < len; i++)
c@241 180 {
c@241 181 temp = data[ i ];
c@241 182
c@241 183 if( temp < *min )
c@241 184 {
c@241 185 *min = temp ;
c@241 186 }
c@241 187 if( temp > *max )
c@241 188 {
c@241 189 *max = temp ;
c@241 190 }
c@241 191
c@241 192 }
c@241 193 }
c@241 194
c@241 195 int MathUtilities::getMax( double* pData, unsigned int Length, double* pMax )
c@241 196 {
c@241 197 unsigned int index = 0;
c@241 198 unsigned int i;
c@241 199 double temp = 0.0;
c@241 200
c@241 201 double max = pData[0];
c@241 202
c@241 203 for( i = 0; i < Length; i++)
c@241 204 {
c@241 205 temp = pData[ i ];
c@241 206
c@241 207 if( temp > max )
c@241 208 {
c@241 209 max = temp ;
c@241 210 index = i;
c@241 211 }
c@241 212
c@241 213 }
c@241 214
c@279 215 if (pMax) *pMax = max;
c@279 216
c@279 217
c@279 218 return index;
c@279 219 }
c@279 220
c@279 221 int MathUtilities::getMax( const std::vector<double> & data, double* pMax )
c@279 222 {
c@279 223 unsigned int index = 0;
c@279 224 unsigned int i;
c@279 225 double temp = 0.0;
c@279 226
c@279 227 double max = data[0];
c@279 228
c@279 229 for( i = 0; i < data.size(); i++)
c@279 230 {
c@279 231 temp = data[ i ];
c@279 232
c@279 233 if( temp > max )
c@279 234 {
c@279 235 max = temp ;
c@279 236 index = i;
c@279 237 }
c@279 238
c@279 239 }
c@279 240
c@279 241 if (pMax) *pMax = max;
c@241 242
c@241 243
c@241 244 return index;
c@241 245 }
c@241 246
c@241 247 void MathUtilities::circShift( double* pData, int length, int shift)
c@241 248 {
c@241 249 shift = shift % length;
c@241 250 double temp;
c@241 251 int i,n;
c@241 252
c@241 253 for( i = 0; i < shift; i++)
c@241 254 {
c@241 255 temp=*(pData + length - 1);
c@241 256
c@241 257 for( n = length-2; n >= 0; n--)
c@241 258 {
c@241 259 *(pData+n+1)=*(pData+n);
c@241 260 }
c@241 261
c@241 262 *pData = temp;
c@241 263 }
c@241 264 }
c@241 265
c@241 266 int MathUtilities::compareInt (const void * a, const void * b)
c@241 267 {
c@241 268 return ( *(int*)a - *(int*)b );
c@241 269 }
c@241 270
c@259 271 void MathUtilities::normalise(double *data, int length, NormaliseType type)
c@259 272 {
c@259 273 switch (type) {
c@259 274
c@259 275 case NormaliseNone: return;
c@259 276
c@259 277 case NormaliseUnitSum:
c@259 278 {
c@259 279 double sum = 0.0;
c@259 280 for (int i = 0; i < length; ++i) {
c@259 281 sum += data[i];
c@259 282 }
c@259 283 if (sum != 0.0) {
c@259 284 for (int i = 0; i < length; ++i) {
c@259 285 data[i] /= sum;
c@259 286 }
c@259 287 }
c@259 288 }
c@259 289 break;
c@259 290
c@259 291 case NormaliseUnitMax:
c@259 292 {
c@259 293 double max = 0.0;
c@259 294 for (int i = 0; i < length; ++i) {
c@259 295 if (fabs(data[i]) > max) {
c@259 296 max = fabs(data[i]);
c@259 297 }
c@259 298 }
c@259 299 if (max != 0.0) {
c@259 300 for (int i = 0; i < length; ++i) {
c@259 301 data[i] /= max;
c@259 302 }
c@259 303 }
c@259 304 }
c@259 305 break;
c@259 306
c@259 307 }
c@259 308 }
c@259 309
c@259 310 void MathUtilities::normalise(std::vector<double> &data, NormaliseType type)
c@259 311 {
c@259 312 switch (type) {
c@259 313
c@259 314 case NormaliseNone: return;
c@259 315
c@259 316 case NormaliseUnitSum:
c@259 317 {
c@259 318 double sum = 0.0;
c@325 319 for (int i = 0; i < (int)data.size(); ++i) sum += data[i];
c@259 320 if (sum != 0.0) {
c@325 321 for (int i = 0; i < (int)data.size(); ++i) data[i] /= sum;
c@259 322 }
c@259 323 }
c@259 324 break;
c@259 325
c@259 326 case NormaliseUnitMax:
c@259 327 {
c@259 328 double max = 0.0;
c@325 329 for (int i = 0; i < (int)data.size(); ++i) {
c@259 330 if (fabs(data[i]) > max) max = fabs(data[i]);
c@259 331 }
c@259 332 if (max != 0.0) {
c@325 333 for (int i = 0; i < (int)data.size(); ++i) data[i] /= max;
c@259 334 }
c@259 335 }
c@259 336 break;
c@259 337
c@259 338 }
c@259 339 }
c@259 340
c@279 341 void MathUtilities::adaptiveThreshold(std::vector<double> &data)
c@279 342 {
c@279 343 int sz = int(data.size());
c@279 344 if (sz == 0) return;
c@279 345
c@279 346 std::vector<double> smoothed(sz);
c@279 347
c@279 348 int p_pre = 8;
c@279 349 int p_post = 7;
c@279 350
c@279 351 for (int i = 0; i < sz; ++i) {
c@279 352
c@279 353 int first = std::max(0, i - p_pre);
c@279 354 int last = std::min(sz - 1, i + p_post);
c@279 355
c@279 356 smoothed[i] = mean(data, first, last - first + 1);
c@279 357 }
c@279 358
c@279 359 for (int i = 0; i < sz; i++) {
c@279 360 data[i] -= smoothed[i];
c@279 361 if (data[i] < 0.0) data[i] = 0.0;
c@279 362 }
c@279 363 }
c@259 364
c@280 365 bool
c@280 366 MathUtilities::isPowerOfTwo(int x)
c@280 367 {
c@280 368 if (x < 2) return false;
c@280 369 if (x & (x-1)) return false;
c@280 370 return true;
c@280 371 }
c@280 372
c@280 373 int
c@280 374 MathUtilities::nextPowerOfTwo(int x)
c@280 375 {
c@280 376 if (isPowerOfTwo(x)) return x;
c@280 377 int n = 1;
c@280 378 while (x) { x >>= 1; n <<= 1; }
c@280 379 return n;
c@280 380 }
c@280 381
c@280 382 int
c@280 383 MathUtilities::previousPowerOfTwo(int x)
c@280 384 {
c@280 385 if (isPowerOfTwo(x)) return x;
c@280 386 int n = 1;
c@280 387 x >>= 1;
c@280 388 while (x) { x >>= 1; n <<= 1; }
c@280 389 return n;
c@280 390 }
c@280 391
c@280 392 int
c@280 393 MathUtilities::nearestPowerOfTwo(int x)
c@280 394 {
c@280 395 if (isPowerOfTwo(x)) return x;
c@280 396 int n0 = previousPowerOfTwo(x), n1 = nearestPowerOfTwo(x);
c@280 397 if (x - n0 < n1 - x) return n0;
c@280 398 else return n1;
c@280 399 }
c@280 400