Mercurial > hg > vamp-fanchirp
comparison FChTransformUtils.cpp @ 10:af59167b3d35 perf
size_t -> int throughout: permits better optimisation for tight for-loops, making the whole thing about 10% faster
author | Chris Cannam |
---|---|
date | Tue, 02 Oct 2018 13:21:15 +0100 |
parents | 7ec763dc767c |
children | 44b86c346a5a |
comparison
equal
deleted
inserted
replaced
9:54dd64b6cfc0 | 10:af59167b3d35 |
---|---|
17 | 17 |
18 | 18 |
19 #include "FChTransformUtils.h" | 19 #include "FChTransformUtils.h" |
20 #include <math.h> | 20 #include <math.h> |
21 | 21 |
22 void cumtrapz(const double *x, const double *y, size_t N, double *accum) | 22 void cumtrapz(const double *x, const double *y, int N, double *accum) |
23 /*Trapezoidal Integrator: 1/2(b-a)(F(a)+F(b))*/ | 23 /*Trapezoidal Integrator: 1/2(b-a)(F(a)+F(b))*/ |
24 { | 24 { |
25 accum[0]=0.0; | 25 accum[0]=0.0; |
26 for (size_t i = 1; i < N; i++) { | 26 for (int i = 1; i < N; i++) { |
27 accum[i]=accum[i-1]+0.5*(x[i]-x[i-1])*(y[i]+y[i-1]); | 27 accum[i]=accum[i-1]+0.5*(x[i]-x[i-1])*(y[i]+y[i-1]); |
28 } | 28 } |
29 } | 29 } |
30 | 30 |
31 | 31 |
32 void interp1(const double *x1, const double *y1, size_t N1, const double *x2, double *y2, size_t N2){ | 32 void interp1(const double *x1, const double *y1, int N1, const double *x2, double *y2, int N2){ |
33 /*1-D linear interpolation*/ | 33 /*1-D linear interpolation*/ |
34 | 34 |
35 for (size_t i = 0; i < N2; i++) { | 35 for (int i = 0; i < N2; i++) { |
36 /*Smaller or equal than the smallest, or larger or equal than the largest.*/ | 36 /*Smaller or equal than the smallest, or larger or equal than the largest.*/ |
37 if ( x2[i] <= x1[0] ) { | 37 if ( x2[i] <= x1[0] ) { |
38 y2[i] = y1[0]; | 38 y2[i] = y1[0]; |
39 } else if ( x2[i] >= x1[N1-1] ) { | 39 } else if ( x2[i] >= x1[N1-1] ) { |
40 y2[i] = y1[N1-1]; | 40 y2[i] = y1[N1-1]; |
41 } else { | 41 } else { |
42 /*Search every value of x2 in x1*/ | 42 /*Search every value of x2 in x1*/ |
43 size_t j = 1; | 43 int j = 1; |
44 size_t salir = 0; | 44 int salir = 0; |
45 while ((j<N1)&&(!salir)) { | 45 while ((j<N1)&&(!salir)) { |
46 if ( x2[i] <= x1[j] ) { | 46 if ( x2[i] <= x1[j] ) { |
47 //y2[i] = ( x2[i]*( y1[j] - y1[j-1] ) + x1[j]*y1[j-1] - x1[j-1]*y1[j] ) / ( x1[j] - x1[j-1] ); | |
48 y2[i] = y1[j-1] + ( ( y1[j] - y1[j-1] )*(x2[i] - x1[j-1] ) )/ ( x1[j] - x1[j-1] ); | 47 y2[i] = y1[j-1] + ( ( y1[j] - y1[j-1] )*(x2[i] - x1[j-1] ) )/ ( x1[j] - x1[j-1] ); |
49 salir = 1; | 48 salir = 1; |
50 } // if | 49 } // if |
51 j++; | 50 j++; |
52 } // for | 51 } // for |
53 } // else | 52 } // else |
54 } // for | 53 } // for |
55 | 54 |
56 } | 55 } |
57 | 56 |
58 void interp1q(const double *y1, const size_t *x2_int, const double *x2_frac, double *y2, size_t N2){ | 57 void interp1q(const double *y1, const int *x2_int, const double *x2_frac, double *y2, int N2){ |
59 | 58 |
60 for(size_t i=0;i<N2;i++){ | 59 for(int i=0;i<N2;i++){ |
61 y2[i] = y1[x2_int[i]]*(1.0-x2_frac[i])+y1[x2_int[i]+1]*x2_frac[i]; | 60 y2[i] = y1[x2_int[i]]*(1.0-x2_frac[i])+y1[x2_int[i]+1]*x2_frac[i]; |
62 } // for | 61 } // for |
63 | 62 |
64 } | 63 } |
65 | 64 |
66 void hanning_window(double *p_window, size_t n, bool normalize) { | 65 void hanning_window(double *p_window, int n, bool normalize) { |
67 | 66 |
68 double accum=0; | 67 double accum=0; |
69 for (size_t i = 0; i < n; i++) { | 68 for (int i = 0; i < n; i++) { |
70 p_window[i] = 0.5*(1.0-cos(2*M_PI*(double)(i+1)/((double)n+1.0))); | 69 p_window[i] = 0.5*(1.0-cos(2*M_PI*(double)(i+1)/((double)n+1.0))); |
71 accum += p_window[i]; | 70 accum += p_window[i]; |
72 } | 71 } |
73 if (normalize) { | 72 if (normalize) { |
74 for (size_t i = 0; i < n; i++) { //window normalization | 73 for (int i = 0; i < n; i++) { //window normalization |
75 p_window[i] = p_window[i]/accum; | 74 p_window[i] = p_window[i]/accum; |
76 } | 75 } |
77 } | 76 } |
78 | 77 |
79 } | 78 } |