[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] [Patch 2/4] Introduce continuous wavelet transform
From: |
Kevin Bube |
Subject: |
[Bug-gsl] [Patch 2/4] Introduce continuous wavelet transform |
Date: |
Fri, 29 Apr 2005 18:38:37 +0200 |
User-agent: |
Gnus/5.1006 (Gnus v5.10.6) Emacs/21.3 (gnu/linux) |
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
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Bug-gsl] [Patch 2/4] Introduce continuous wavelet transform,
Kevin Bube <=