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