code/src/Algorithm/compute_coefficients.c

   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 }