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