libs/qm-dsp/maths/MathUtilities.cpp

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