qm-dsp  1.8
MathUtilities.cpp
Go to the documentation of this file.
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
346 {
347  if (x < 1) return false;
348  if (x & (x-1)) return false;
349  return true;
350 }
351 
352 int
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
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
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
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 
static void adaptiveThreshold(std::vector< double > &data)
Threshold the input/output vector data against a moving-mean average filter.
static void circShift(double *data, int length, int shift)
static double mod(double x, double y)
Floating-point division modulus: return x % y.
static double princarg(double ang)
The principle argument function.
static double mean(const double *src, unsigned int len)
Return the mean of the given array of the given length.
static double sum(const double *src, unsigned int len)
Return the sum of the values in the given array of the given length.
static double median(const double *src, unsigned int len)
Return the median of the values in the given array of the given length.
static double factorial(int x)
Return x!
static void getFrameMinMax(const double *data, unsigned int len, double *min, double *max)
Return through min and max pointers the highest and lowest values in the given array of the given len...
static int nextPowerOfTwo(int x)
Return the next higher integer power of two from x, e.g.
static int compareInt(const void *a, const void *b)
static void getAlphaNorm(const double *data, unsigned int len, unsigned int alpha, double *ANorm)
static bool isPowerOfTwo(int x)
Return true if x is 2^n for some integer n >= 0.
static int previousPowerOfTwo(int x)
Return the next lower integer power of two from x, e.g.
static int gcd(int a, int b)
Return the greatest common divisor of natural numbers a and b.
static int nearestPowerOfTwo(int x)
Return the nearest integer power of two to x, e.g.
static double round(double x)
Round x to the nearest integer.
static int getMax(double *data, unsigned int length, double *max=0)
static void normalise(double *data, int length, NormaliseType n=NormaliseUnitMax)