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 }