← Back to branch summary

~ubuntu-branches/ubuntu/utopic/ardour3/utopic

« back to all changes in this revision

Viewing changes to libs/qm-dsp/maths/MathUtilities.cpp

  • Committer: Package Import Robot
  • Author(s): Felipe Sateler
  • Date: 2013-09-21 19:05:02 UTC
  • Revision ID: package-import@ubuntu.com-20130921190502-8gsftrku6jnzhd7v
Tags: upstream-3.4~dfsg
Import upstream version 3.4~dfsg

Show diffs side-by-side

added added

removed removed

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