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 using namespace std;
  24
  25 double MathUtilities::mod(double x, double y)
  26 {
  27     double a = floor( x / y );
  28
  29     double b = x - ( y * a );
  30     return b;
  31 }
  32
  33 double MathUtilities::princarg(double ang)
  34 {
  35     double ValOut;
  36
  37     ValOut = mod( ang + M_PI, -2 * M_PI ) + M_PI;
  38
  39     return ValOut;
  40 }
  41
  42 void MathUtilities::getAlphaNorm(const double *data, int len, int alpha, double* ANorm)
  43 {
  44     int i;
  45     double temp = 0.0;
  46     double a=0.0;
  47
  48     for( i = 0; i < len; i++)
  49     {
  50         temp = data[ i ];
  51
  52         a  += ::pow( fabs(temp), double(alpha) );
  53     }
  54     a /= ( double )len;
  55     a = ::pow( a, ( 1.0 / (double) alpha ) );
  56
  57     *ANorm = a;
  58 }
  59
  60 double MathUtilities::getAlphaNorm( const vector <double> &data, int alpha )
  61 {
  62     int i;
  63     int len = data.size();
  64     double temp = 0.0;
  65     double a=0.0;
  66
  67     for( i = 0; i < len; i++)
  68     {
  69         temp = data[ i ];
  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, int len)
  88 {
  89     if (len == 0) return 0;
  90
  91     vector<double> scratch;
  92     for (int i = 0; i < len; ++i) scratch.push_back(src[i]);
  93     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, int len)
 104 {
 105     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, 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 vector<double> &src,
 130                            int start,
 131                            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, int len, double *min, double *max)
 146 {
 147     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, int Length, double* pMax )
 175 {
 176         int index = 0;
 177         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 vector<double> & data, double* pMax )
 201 {
 202         int index = 0;
 203         int i;
 204         double temp = 0.0;
 205
 206         double max = data[0];
 207
 208         for( i = 0; i < int(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(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 double MathUtilities::getLpNorm(const vector<double> &data, int p)
 321 {
 322     double tot = 0.0;
 323     for (int i = 0; i < int(data.size()); ++i) {
 324         tot += abs(pow(data[i], p));
 325     }
 326     return pow(tot, 1.0 / p);
 327 }
 328
 329 vector<double> MathUtilities::normaliseLp(const vector<double> &data,
 330                                                int p,
 331                                                double threshold)
 332 {
 333     int n = int(data.size());
 334     if (n == 0 || p == 0) return data;
 335     double norm = getLpNorm(data, p);
 336     if (norm < threshold) {
 337         return vector<double>(n, 1.0 / pow(n, 1.0 / p)); // unit vector
 338     }
 339     vector<double> out(n);
 340     for (int i = 0; i < n; ++i) {
 341         out[i] = data[i] / norm;
 342     }
 343     return out;
 344 }
 345
 346 void MathUtilities::adaptiveThreshold(vector<double> &data)
 347 {
 348     int sz = int(data.size());
 349     if (sz == 0) return;
 350
 351     vector<double> smoothed(sz);
 352
 353     int p_pre = 8;
 354     int p_post = 7;
 355
 356     for (int i = 0; i < sz; ++i) {
 357
 358         int first = max(0,      i - p_pre);
 359         int last  = min(sz - 1, i + p_post);
 360
 361         smoothed[i] = mean(data, first, last - first + 1);
 362     }
 363
 364     for (int i = 0; i < sz; i++) {
 365         data[i] -= smoothed[i];
 366         if (data[i] < 0.0) data[i] = 0.0;
 367     }
 368 }
 369
 370 bool
 371 MathUtilities::isPowerOfTwo(int x)
 372 {
 373     if (x < 1) return false;
 374     if (x & (x-1)) return false;
 375     return true;
 376 }
 377
 378 int
 379 MathUtilities::nextPowerOfTwo(int x)
 380 {
 381     if (isPowerOfTwo(x)) return x;
 382     if (x < 1) return 1;
 383     int n = 1;
 384     while (x) { x >>= 1; n <<= 1; }
 385     return n;
 386 }
 387
 388 int
 389 MathUtilities::previousPowerOfTwo(int x)
 390 {
 391     if (isPowerOfTwo(x)) return x;
 392     if (x < 1) return 1;
 393     int n = 1;
 394     x >>= 1;
 395     while (x) { x >>= 1; n <<= 1; }
 396     return n;
 397 }
 398
 399 int
 400 MathUtilities::nearestPowerOfTwo(int x)
 401 {
 402     if (isPowerOfTwo(x)) return x;
 403     int n0 = previousPowerOfTwo(x), n1 = nextPowerOfTwo(x);
 404     if (x - n0 < n1 - x) return n0;
 405     else return n1;
 406 }
 407
 408 double
 409 MathUtilities::factorial(int x)
 410 {
 411     if (x < 0) return 0;
 412     double f = 1;
 413     for (int i = 1; i <= x; ++i) {
 414         f = f * i;
 415     }
 416     return f;
 417 }
 418
 419 int
 420 MathUtilities::gcd(int a, int b)
 421 {
 422     int c = a % b;
 423     if (c == 0) {
 424         return b;
 425     } else {
 426         return gcd(b, c);
 427     }
 428 }
 429