qm-dsp  1.8
MathUtilities.cpp
Go to the documentation of this file.
00001 /* -*- c-basic-offset: 4 indent-tabs-mode: nil -*-  vi:set ts=8 sts=4 sw=4: */
00002 
00003 /*
00004     QM DSP Library
00005 
00006     Centre for Digital Music, Queen Mary, University of London.
00007     This file 2005-2006 Christian Landone.
00008 
00009     This program is free software; you can redistribute it and/or
00010     modify it under the terms of the GNU General Public License as
00011     published by the Free Software Foundation; either version 2 of the
00012     License, or (at your option) any later version.  See the file
00013     COPYING included with this distribution for more information.
00014 */
00015 
00016 #include "MathUtilities.h"
00017 
00018 #include <iostream>
00019 #include <algorithm>
00020 #include <vector>
00021 #include <cmath>
00022 
00023 
00024 double MathUtilities::mod(double x, double y)
00025 {
00026     double a = floor( x / y );
00027 
00028     double b = x - ( y * a );
00029     return b;
00030 }
00031 
00032 double MathUtilities::princarg(double ang)
00033 {
00034     double ValOut;
00035 
00036     ValOut = mod( ang + M_PI, -2 * M_PI ) + M_PI;
00037 
00038     return ValOut;
00039 }
00040 
00041 void MathUtilities::getAlphaNorm(const double *data, unsigned int len, unsigned int alpha, double* ANorm)
00042 {
00043     unsigned int i;
00044     double temp = 0.0;
00045     double a=0.0;
00046         
00047     for( i = 0; i < len; i++)
00048     {
00049         temp = data[ i ];
00050                 
00051         a  += ::pow( fabs(temp), double(alpha) );
00052     }
00053     a /= ( double )len;
00054     a = ::pow( a, ( 1.0 / (double) alpha ) );
00055 
00056     *ANorm = a;
00057 }
00058 
00059 double MathUtilities::getAlphaNorm( const std::vector <double> &data, unsigned int alpha )
00060 {
00061     unsigned int i;
00062     unsigned int len = data.size();
00063     double temp = 0.0;
00064     double a=0.0;
00065         
00066     for( i = 0; i < len; i++)
00067     {
00068         temp = data[ i ];
00069                 
00070         a  += ::pow( fabs(temp), double(alpha) );
00071     }
00072     a /= ( double )len;
00073     a = ::pow( a, ( 1.0 / (double) alpha ) );
00074 
00075     return a;
00076 }
00077 
00078 double MathUtilities::round(double x)
00079 {
00080     if (x < 0) {
00081         return -floor(-x + 0.5);
00082     } else {
00083         return floor(x + 0.5);
00084     }
00085 }
00086 
00087 double MathUtilities::median(const double *src, unsigned int len)
00088 {
00089     if (len == 0) return 0;
00090     
00091     std::vector<double> scratch;
00092     for (int i = 0; i < len; ++i) scratch.push_back(src[i]);
00093     std::sort(scratch.begin(), scratch.end());
00094 
00095     int middle = len/2;
00096     if (len % 2 == 0) {
00097         return (scratch[middle] + scratch[middle - 1]) / 2;
00098     } else {
00099         return scratch[middle];
00100     }
00101 }
00102 
00103 double MathUtilities::sum(const double *src, unsigned int len)
00104 {
00105     unsigned int i ;
00106     double retVal =0.0;
00107 
00108     for(  i = 0; i < len; i++)
00109     {
00110         retVal += src[ i ];
00111     }
00112 
00113     return retVal;
00114 }
00115 
00116 double MathUtilities::mean(const double *src, unsigned int len)
00117 {
00118     double retVal =0.0;
00119 
00120     if (len == 0) return 0;
00121 
00122     double s = sum( src, len );
00123         
00124     retVal =  s  / (double)len;
00125 
00126     return retVal;
00127 }
00128 
00129 double MathUtilities::mean(const std::vector<double> &src,
00130                            unsigned int start,
00131                            unsigned int count)
00132 {
00133     double sum = 0.;
00134         
00135     if (count == 0) return 0;
00136     
00137     for (int i = 0; i < (int)count; ++i)
00138     {
00139         sum += src[start + i];
00140     }
00141 
00142     return sum / count;
00143 }
00144 
00145 void MathUtilities::getFrameMinMax(const double *data, unsigned int len, double *min, double *max)
00146 {
00147     unsigned int i;
00148     double temp = 0.0;
00149 
00150     if (len == 0) {
00151         *min = *max = 0;
00152         return;
00153     }
00154         
00155     *min = data[0];
00156     *max = data[0];
00157 
00158     for( i = 0; i < len; i++)
00159     {
00160         temp = data[ i ];
00161 
00162         if( temp < *min )
00163         {
00164             *min =  temp ;
00165         }
00166         if( temp > *max )
00167         {
00168             *max =  temp ;
00169         }
00170                 
00171     }
00172 }
00173 
00174 int MathUtilities::getMax( double* pData, unsigned int Length, double* pMax )
00175 {
00176         unsigned int index = 0;
00177         unsigned int i;
00178         double temp = 0.0;
00179         
00180         double max = pData[0];
00181 
00182         for( i = 0; i < Length; i++)
00183         {
00184                 temp = pData[ i ];
00185 
00186                 if( temp > max )
00187                 {
00188                         max =  temp ;
00189                         index = i;
00190                 }
00191                 
00192         }
00193 
00194         if (pMax) *pMax = max;
00195 
00196 
00197         return index;
00198 }
00199 
00200 int MathUtilities::getMax( const std::vector<double> & data, double* pMax )
00201 {
00202         unsigned int index = 0;
00203         unsigned int i;
00204         double temp = 0.0;
00205         
00206         double max = data[0];
00207 
00208         for( i = 0; i < data.size(); i++)
00209         {
00210                 temp = data[ i ];
00211 
00212                 if( temp > max )
00213                 {
00214                         max =  temp ;
00215                         index = i;
00216                 }
00217                 
00218         }
00219 
00220         if (pMax) *pMax = max;
00221 
00222 
00223         return index;
00224 }
00225 
00226 void MathUtilities::circShift( double* pData, int length, int shift)
00227 {
00228         shift = shift % length;
00229         double temp;
00230         int i,n;
00231 
00232         for( i = 0; i < shift; i++)
00233         {
00234                 temp=*(pData + length - 1);
00235 
00236                 for( n = length-2; n >= 0; n--)
00237                 {
00238                         *(pData+n+1)=*(pData+n);
00239                 }
00240 
00241         *pData = temp;
00242     }
00243 }
00244 
00245 int MathUtilities::compareInt (const void * a, const void * b)
00246 {
00247   return ( *(int*)a - *(int*)b );
00248 }
00249 
00250 void MathUtilities::normalise(double *data, int length, NormaliseType type)
00251 {
00252     switch (type) {
00253 
00254     case NormaliseNone: return;
00255 
00256     case NormaliseUnitSum:
00257     {
00258         double sum = 0.0;
00259         for (int i = 0; i < length; ++i) {
00260             sum += data[i];
00261         }
00262         if (sum != 0.0) {
00263             for (int i = 0; i < length; ++i) {
00264                 data[i] /= sum;
00265             }
00266         }
00267     }
00268     break;
00269 
00270     case NormaliseUnitMax:
00271     {
00272         double max = 0.0;
00273         for (int i = 0; i < length; ++i) {
00274             if (fabs(data[i]) > max) {
00275                 max = fabs(data[i]);
00276             }
00277         }
00278         if (max != 0.0) {
00279             for (int i = 0; i < length; ++i) {
00280                 data[i] /= max;
00281             }
00282         }
00283     }
00284     break;
00285 
00286     }
00287 }
00288 
00289 void MathUtilities::normalise(std::vector<double> &data, NormaliseType type)
00290 {
00291     switch (type) {
00292 
00293     case NormaliseNone: return;
00294 
00295     case NormaliseUnitSum:
00296     {
00297         double sum = 0.0;
00298         for (int i = 0; i < (int)data.size(); ++i) sum += data[i];
00299         if (sum != 0.0) {
00300             for (int i = 0; i < (int)data.size(); ++i) data[i] /= sum;
00301         }
00302     }
00303     break;
00304 
00305     case NormaliseUnitMax:
00306     {
00307         double max = 0.0;
00308         for (int i = 0; i < (int)data.size(); ++i) {
00309             if (fabs(data[i]) > max) max = fabs(data[i]);
00310         }
00311         if (max != 0.0) {
00312             for (int i = 0; i < (int)data.size(); ++i) data[i] /= max;
00313         }
00314     }
00315     break;
00316 
00317     }
00318 }
00319 
00320 void MathUtilities::adaptiveThreshold(std::vector<double> &data)
00321 {
00322     int sz = int(data.size());
00323     if (sz == 0) return;
00324 
00325     std::vector<double> smoothed(sz);
00326         
00327     int p_pre = 8;
00328     int p_post = 7;
00329 
00330     for (int i = 0; i < sz; ++i) {
00331 
00332         int first = std::max(0,      i - p_pre);
00333         int last  = std::min(sz - 1, i + p_post);
00334 
00335         smoothed[i] = mean(data, first, last - first + 1);
00336     }
00337 
00338     for (int i = 0; i < sz; i++) {
00339         data[i] -= smoothed[i];
00340         if (data[i] < 0.0) data[i] = 0.0;
00341     }
00342 }
00343 
00344 bool
00345 MathUtilities::isPowerOfTwo(int x)
00346 {
00347     if (x < 1) return false;
00348     if (x & (x-1)) return false;
00349     return true;
00350 }
00351 
00352 int
00353 MathUtilities::nextPowerOfTwo(int x)
00354 {
00355     if (isPowerOfTwo(x)) return x;
00356     if (x < 1) return 1;
00357     int n = 1;
00358     while (x) { x >>= 1; n <<= 1; }
00359     return n;
00360 }
00361 
00362 int
00363 MathUtilities::previousPowerOfTwo(int x)
00364 {
00365     if (isPowerOfTwo(x)) return x;
00366     if (x < 1) return 1;
00367     int n = 1;
00368     x >>= 1;
00369     while (x) { x >>= 1; n <<= 1; }
00370     return n;
00371 }
00372 
00373 int
00374 MathUtilities::nearestPowerOfTwo(int x)
00375 {
00376     if (isPowerOfTwo(x)) return x;
00377     int n0 = previousPowerOfTwo(x), n1 = nextPowerOfTwo(x);
00378     if (x - n0 < n1 - x) return n0;
00379     else return n1;
00380 }
00381 
00382 double
00383 MathUtilities::factorial(int x)
00384 {
00385     if (x < 0) return 0;
00386     double f = 1;
00387     for (int i = 1; i <= x; ++i) {
00388         f = f * i;
00389     }
00390     return f;
00391 }
00392 
00393 int
00394 MathUtilities::gcd(int a, int b)
00395 {
00396     int c = a % b;
00397     if (c == 0) {
00398         return b;
00399     } else {
00400         return gcd(b, c);
00401     }
00402 }
00403