Commit | Line | Data |
---|---|---|
81923e5c BA |
1 | #include <stdlib.h> |
2 | #include "Util/types.h" | |
3 | #include <math.h> | |
4 | #include <gsl/gsl_wavelet.h> | |
5 | #include <gsl/gsl_spline.h> | |
6 | ||
7 | // compute rows of the matrix of reduced coordinates | |
8 | void compute_coefficients(PowerCurve* powerCurves, uint32_t nbSeries, uint32_t nbValues, | |
9 | Real* reducedCoordinates, uint32_t index, uint32_t nbReducedCoordinates) | |
10 | { | |
11 | uint32_t D = (1 << nbReducedCoordinates); | |
12 | Real* x = (Real*) malloc(nbValues*sizeof(Real)); | |
13 | for (uint32_t i=0; i<nbValues; i++) | |
14 | x[i] = i; | |
15 | Real* y = (Real*) malloc(nbValues*sizeof(Real)); | |
16 | Real* interpolatedTranformedCurve = (Real*) malloc(D*sizeof(Real)); | |
17 | gsl_interp* linearInterpolation = gsl_interp_alloc(gsl_interp_linear, nbValues); | |
18 | gsl_interp_accel* acc = gsl_interp_accel_alloc(); | |
19 | gsl_wavelet_workspace* work = gsl_wavelet_workspace_alloc(D); | |
20 | //gsl_wavelet* w = gsl_wavelet_alloc(gsl_wavelet_bspline, 206); //used for starlight | |
21 | gsl_wavelet* w = gsl_wavelet_alloc(gsl_wavelet_daubechies, 6); //used for power curves | |
22 | //gsl_wavelet* w = gsl_wavelet_alloc(gsl_wavelet_haar, 2); | |
23 | for (uint32_t i = 0; i < nbSeries; i++) | |
24 | { | |
25 | //Spline interpolation to have D = 2^u sample points | |
26 | for (uint32_t j=0; j<nbValues; j++) | |
27 | y[j] = powerCurves[i].values[j]; | |
28 | gsl_interp_init(linearInterpolation, x, y, nbValues); | |
29 | for (uint32_t j=0; j<D; j++) | |
30 | { | |
31 | interpolatedTranformedCurve[j] = | |
32 | gsl_interp_eval(linearInterpolation, x, y, j*((Real)(nbValues-1)/(D-1)), acc); | |
33 | } | |
34 | //DWT transform (in place) on interpolated curve [TODO: clarify stride parameter] | |
35 | gsl_wavelet_transform_forward(w, interpolatedTranformedCurve, 1, D, work); | |
36 | //Fill reducedCoordinates with energy contributions | |
37 | uint32_t t0 = 1; | |
38 | uint32_t t1 = 1; | |
39 | for (uint32_t j=0; j<nbReducedCoordinates; j++) | |
40 | { | |
41 | t1 += (1 << j); | |
42 | //reducedCoordinates[(index+i)*nbReducedCoordinates+j] = sqrt( sum( x[t0:t1]^2 ) ) | |
43 | Real sumOnSegment = 0.0; | |
44 | for (uint32_t u=t0; u<t1; u++) | |
45 | sumOnSegment += interpolatedTranformedCurve[u]*interpolatedTranformedCurve[u]; | |
46 | reducedCoordinates[(index+i)*nbReducedCoordinates+j] = sqrt(sumOnSegment); | |
47 | t0 = t1; | |
48 | } | |
49 | } | |
50 | gsl_interp_free(linearInterpolation); | |
51 | gsl_interp_accel_free(acc); | |
52 | free(x); | |
53 | free(y); | |
54 | free(interpolatedTranformedCurve); | |
55 | gsl_wavelet_free(w); | |
56 | gsl_wavelet_workspace_free(work); | |
57 | } |