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
24double 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
32double 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
41void 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
59double 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
78double 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
87double 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
103double 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
116double 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
129double 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
145void 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
174int 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
200int 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
226void 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
245int MathUtilities::compareInt (const void * a, const void * b)
246{
247 return ( *(int*)a - *(int*)b );
248}
249
250void 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
289void 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
320void 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
344bool
346{
347 if (x < 1) return false;
348 if (x & (x-1)) return false;
349 return true;
350}
351
352int
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
362int
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
373int
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
382double
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
393int
395{
396 int c = a % b;
397 if (c == 0) {
398 return b;
399 } else {
400 return gcd(b, c);
401 }
402}
403
static int nearestPowerOfTwo(int x)
Return the nearest integer power of two to x, e.g.
static void getAlphaNorm(const double *data, unsigned int len, unsigned int alpha, double *ANorm)
static int compareInt(const void *a, const void *b)
static double sum(const double *src, unsigned int len)
Return the sum of the values in the given array of the given length.
static void normalise(double *data, int length, NormaliseType n=NormaliseUnitMax)
static double mod(double x, double y)
Floating-point division modulus: return x % y.
static int nextPowerOfTwo(int x)
Return the next higher integer power of two from x, e.g.
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 mean(const double *src, unsigned int len)
Return the mean of the given array of the given length.
static bool isPowerOfTwo(int x)
Return true if x is 2^n for some integer n >= 0.
static int gcd(int a, int b)
Return the greatest common divisor of natural numbers a and b.
static void adaptiveThreshold(std::vector< double > &data)
Threshold the input/output vector data against a moving-mean average filter.
static double factorial(int x)
Return x!
static double round(double x)
Round x to the nearest integer.
static int getMax(double *data, unsigned int length, double *max=0)
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 void circShift(double *data, int length, int shift)
static int previousPowerOfTwo(int x)
Return the next lower integer power of two from x, e.g.
static double princarg(double ang)
The principle argument function.