[Bug-gsl] [Patch 2/4] Introduce continuous wavelet transform

[Kevin Bube]

This patch provides the wavelets.


--- /dev/null   2004-04-06 15:27:52.000000000 +0200
+++ gsl-1.6/wavelet/cwt_wavelets.c      2005-04-29 17:24:42.000000000 +0200
@@ -0,0 +1,126 @@
+/* wavelet/cwt_wavelets.c
+ * 
+ * Copyright (C) 2005 Kevin Bube
+ * 
+ * This program is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation; either version 2 of the License, or (at
+ * your option) any later version.
+ * 
+ * This program is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * General Public License for more details.
+ * 
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ */
+
+#include <math.h>
+#include <gsl/gsl_cwt.h>
+#include <gsl/gsl_complex.h>
+#include <gsl/gsl_complex_math.h>
+#include <gsl/gsl_sf_gamma.h>
+
+#define HEAVISIDE_FUNC(arg) ((arg)>0 ? 1 : 0)
+
+static gsl_complex gsl_cwt_dog_wavelet (double scale, double frequency,
+                                       double *parms);
+static gsl_complex gsl_cwt_morlet_wavelet (double scale, double frequency,
+                                          double *parms);
+static gsl_complex gsl_cwt_paul_wavelet (double scale, double frequency,
+                                        double *parms);
+
+/***********************************************************************
+ *
+ * derivative of gaussian wavelet
+ *
+ * One parameter: order of derivative
+ *
+ ***********************************************************************/
+
+static const gsl_cwt_type gsl_cwt_dog_type = {
+  "DOG",
+  &gsl_cwt_dog_wavelet
+};
+
+const gsl_cwt_type *gsl_cwt_dog = &gsl_cwt_dog_type;
+
+static gsl_complex
+gsl_cwt_dog_wavelet (double scale, double frequency, double *parms)
+{
+  double tmp;
+  gsl_complex ret_val;
+
+  tmp = pow(scale*frequency, parms[0]) * exp(-pow(frequency*scale, 2.0)/2);
+  tmp = -tmp / sqrt(gsl_sf_gamma(parms[0] + 0.5));
+    
+  GSL_SET_COMPLEX(&ret_val, 0, 1);
+  ret_val = gsl_complex_pow_real(ret_val, parms[0]);
+
+  return gsl_complex_mul_real(ret_val, tmp);
+}
+
+
+/***********************************************************************
+ *
+ * Morlet wavelet
+ *
+ * One parameter: omega_0
+ *
+ ***********************************************************************/
+
+static const gsl_cwt_type gsl_cwt_morlet_type = {
+  "Morlet",
+  &gsl_cwt_morlet_wavelet
+};
+
+const gsl_cwt_type *gsl_cwt_morlet = &gsl_cwt_morlet_type;
+
+static gsl_complex
+gsl_cwt_morlet_wavelet (double scale, double frequency, double *parms)
+{
+  gsl_complex ret_val;
+  double tmp;
+  
+  tmp = (pow(M_PI, -0.25) * HEAVISIDE_FUNC(frequency) *
+        exp(-(pow(scale*frequency - parms[0], 2.0)/2)));
+
+  GSL_SET_COMPLEX(&ret_val, tmp, 0);
+  
+  return ret_val;
+}
+
+
+/***********************************************************************
+ *
+ * Paul wavelet
+ *
+ * One parameter: order
+ *
+ ***********************************************************************/
+
+static const gsl_cwt_type gsl_cwt_paul_type = {
+  "Paul",
+  &gsl_cwt_paul_wavelet
+};
+
+const gsl_cwt_type *gsl_cwt_paul = &gsl_cwt_paul_type;
+
+static gsl_complex
+gsl_cwt_paul_wavelet (double scale, double frequency, double *parms)
+{
+  gsl_complex ret_val;
+  double tmp;
+  
+  tmp = ((pow(2.0, parms[0]) /
+         sqrt(parms[0] * gsl_sf_fact(2*((unsigned int)parms[0])-1))) *
+        HEAVISIDE_FUNC(frequency) * pow(scale*frequency, parms[0]) *
+        exp(-scale*frequency));
+        
+
+  GSL_SET_COMPLEX(&ret_val, tmp, 0);
+  
+  return ret_val;
+}

pgpqOfAb1FGy1.pgp
Description: PGP signature