| 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 | } |