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 (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     double a=0.0;
 171
 172     if (len == 0) {
 173         *min = *max = 0;
 174         return;
 175     }
 176
 177     *min = data[0];
 178     *max = data[0];
 179
 180     for( i = 0; i < len; i++)
 181     {
 182         temp = data[ i ];
 183
 184         if( temp < *min )
 185         {
 186             *min =  temp ;
 187         }
 188         if( temp > *max )
 189         {
 190             *max =  temp ;
 191         }
 192
 193     }
 194 }
 195
 196 int MathUtilities::getMax( double* pData, unsigned int Length, double* pMax )
 197 {
 198         unsigned int index = 0;
 199         unsigned int i;
 200         double temp = 0.0;
 201
 202         double max = pData[0];
 203
 204         for( i = 0; i < Length; i++)
 205         {
 206                 temp = pData[ i ];
 207
 208                 if( temp > max )
 209                 {
 210                         max =  temp ;
 211                         index = i;
 212                 }
 213
 214         }
 215
 216         if (pMax) *pMax = max;
 217
 218
 219         return index;
 220 }
 221
 222 int MathUtilities::getMax( const std::vector<double> & data, double* pMax )
 223 {
 224         unsigned int index = 0;
 225         unsigned int i;
 226         double temp = 0.0;
 227
 228         double max = data[0];
 229
 230         for( i = 0; i < data.size(); i++)
 231         {
 232                 temp = data[ i ];
 233
 234                 if( temp > max )
 235                 {
 236                         max =  temp ;
 237                         index = i;
 238                 }
 239
 240         }
 241
 242         if (pMax) *pMax = max;
 243
 244
 245         return index;
 246 }
 247
 248 void MathUtilities::circShift( double* pData, int length, int shift)
 249 {
 250         shift = shift % length;
 251         double temp;
 252         int i,n;
 253
 254         for( i = 0; i < shift; i++)
 255         {
 256                 temp=*(pData + length - 1);
 257
 258                 for( n = length-2; n >= 0; n--)
 259                 {
 260                         *(pData+n+1)=*(pData+n);
 261                 }
 262
 263         *pData = temp;
 264     }
 265 }
 266
 267 int MathUtilities::compareInt (const void * a, const void * b)
 268 {
 269   return ( *(int*)a - *(int*)b );
 270 }
 271
 272 void MathUtilities::normalise(double *data, int length, NormaliseType type)
 273 {
 274     switch (type) {
 275
 276     case NormaliseNone: return;
 277
 278     case NormaliseUnitSum:
 279     {
 280         double sum = 0.0;
 281         for (int i = 0; i < length; ++i) {
 282             sum += data[i];
 283         }
 284         if (sum != 0.0) {
 285             for (int i = 0; i < length; ++i) {
 286                 data[i] /= sum;
 287             }
 288         }
 289     }
 290     break;
 291
 292     case NormaliseUnitMax:
 293     {
 294         double max = 0.0;
 295         for (int i = 0; i < length; ++i) {
 296             if (fabs(data[i]) > max) {
 297                 max = fabs(data[i]);
 298             }
 299         }
 300         if (max != 0.0) {
 301             for (int i = 0; i < length; ++i) {
 302                 data[i] /= max;
 303             }
 304         }
 305     }
 306     break;
 307
 308     }
 309 }
 310
 311 void MathUtilities::normalise(std::vector<double> &data, NormaliseType type)
 312 {
 313     switch (type) {
 314
 315     case NormaliseNone: return;
 316
 317     case NormaliseUnitSum:
 318     {
 319         double sum = 0.0;
 320         for (int i = 0; i < data.size(); ++i) sum += data[i];
 321         if (sum != 0.0) {
 322             for (int i = 0; i < data.size(); ++i) data[i] /= sum;
 323         }
 324     }
 325     break;
 326
 327     case NormaliseUnitMax:
 328     {
 329         double max = 0.0;
 330         for (int i = 0; i < data.size(); ++i) {
 331             if (fabs(data[i]) > max) max = fabs(data[i]);
 332         }
 333         if (max != 0.0) {
 334             for (int i = 0; i < data.size(); ++i) data[i] /= max;
 335         }
 336     }
 337     break;
 338
 339     }
 340 }
 341
 342 void MathUtilities::adaptiveThreshold(std::vector<double> &data)
 343 {
 344     int sz = int(data.size());
 345     if (sz == 0) return;
 346
 347     std::vector<double> smoothed(sz);
 348
 349     int p_pre = 8;
 350     int p_post = 7;
 351
 352     for (int i = 0; i < sz; ++i) {
 353
 354         int first = std::max(0,      i - p_pre);
 355         int last  = std::min(sz - 1, i + p_post);
 356
 357         smoothed[i] = mean(data, first, last - first + 1);
 358     }
 359
 360     for (int i = 0; i < sz; i++) {
 361         data[i] -= smoothed[i];
 362         if (data[i] < 0.0) data[i] = 0.0;
 363     }
 364 }
 365
 366 bool
 367 MathUtilities::isPowerOfTwo(int x)
 368 {
 369     if (x < 2) return false;
 370     if (x & (x-1)) return false;
 371     return true;
 372 }
 373
 374 int
 375 MathUtilities::nextPowerOfTwo(int x)
 376 {
 377     if (isPowerOfTwo(x)) return x;
 378     int n = 1;
 379     while (x) { x >>= 1; n <<= 1; }
 380     return n;
 381 }
 382
 383 int
 384 MathUtilities::previousPowerOfTwo(int x)
 385 {
 386     if (isPowerOfTwo(x)) return x;
 387     int n = 1;
 388     x >>= 1;
 389     while (x) { x >>= 1; n <<= 1; }
 390     return n;
 391 }
 392
 393 int
 394 MathUtilities::nearestPowerOfTwo(int x)
 395 {
 396     if (isPowerOfTwo(x)) return x;
 397     int n0 = previousPowerOfTwo(x), n1 = nearestPowerOfTwo(x);
 398     if (x - n0 < n1 - x) return n0;
 399     else return n1;
 400 }
 401