Merge branch 'wavelet-distance' into 'master'

added function to compute wavelet distance

See merge request lscsoft/bayeswave!154

src/BayesWaveWavelet.c

@@ -134,6 +134,46 @@ void SineGaussianTime(double *hs, double *sigpar, int N, int flag, double Tobs)
   }
 }
 
+double SineGaussianDistance(double *params0, double *params1)
+{
+  int i;
+  
+  double t[2];
+  double f[2];
+  double Q[2];
+  double tau[2];
+  
+  t[0] = params0[0];
+  t[1] = params1[0];
+
+  f[0] = params0[1];
+  f[1] = params1[1];
+
+  Q[0] = params0[2];
+  Q[1] = params1[2];
+
+  double tau_min = 2.e-3;
+  double tau_max = 0.5;
+  
+  for(i=0; i<2; i++)
+  {
+    tau[i] = Q[i]/(LAL_TWOPI*f[i]);
+    if(tau[i] < tau_min) tau[i] = tau_min;
+    if(tau[i] > tau_max) tau[i] = tau_max;
+  }
+  
+  //principle axes of wavelet
+  double x = tau[0]*tau[1];
+  double y = tau[0]*tau[0] + tau[1]*tau[1];
+  
+  //metric distance between wavelets
+  double dt = t[1]-t[0];
+  double df = M_PI*x*(f[1]-f[0]);
+  double d  = (dt*dt + df*df)/y + 0.5*log(y/(2.0*x));
+  
+  return d;
+}
+
 void SineGaussianFourier(double *hs, double *sigpar, int N, int flag, double Tobs)
 {
    double f0, t0, Q, Amp;

src/BayesWaveWavelet.h

@@ -25,6 +25,7 @@
 /* ********************************************************************************** */
 
 double SineGaussianSNR(double *params, double *Snf, double Tobs);
+double SineGaussianDistance(double *params0, double *params1);
 void SineGaussianBandwidth(double *sigpar, double *fmin, double *fmax);
 void SineGaussianTime(double *hs, double *sigpar, int N, int flag, double Tobs);
 void SineGaussianFourier(double *hs, double *sigpar, int N, int flag, double Tobs);