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

[Kevin Bube]

This patch provides the actual user interface functions.

-- 
publickey 2048R/0AFDFB19: http://www.icbm.de/~bube/publickey.asc
fingerprint: 542B 1378 04AA AF1F 572E  78BF 1BF5 5C71 0AFD FB19

--- /dev/null   2004-04-06 15:27:52.000000000 +0200
+++ gsl-1.6/wavelet/cwt.c       2005-04-29 17:24:42.000000000 +0200
@@ -0,0 +1,136 @@
+/* wavelet/cwt.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.
+ */
+
+/* Reference: C. Torrence and G. P. Compo: A Practical Guide to Wavelet
+ *            Analysis; Bull. Amer. Meteor. Soc. 79, 61--78 (1998)
+ */
+
+#include <stdlib.h>
+#include <config.h>
+#include <gsl/gsl_errno.h>
+#include <gsl/gsl_cwt.h>
+#include <gsl/gsl_fft_complex.h>
+#include <gsl/gsl_complex.h>
+#include <gsl/gsl_complex_math.h>
+
+#define ELEMENT(a,stride,i) ((a)[(stride)*(i)])
+
+int
+gsl_cwt_transform (const gsl_cwt_wavelet *wlet, double *data, double *imag,
+                  size_t stride, size_t n, gsl_cwt_direction dir,
+                  gsl_cwt_workspace *w)
+
+{
+  if (dir == gsl_cwt_forward)
+    return gsl_cwt_transform_forward (wlet, data, imag, stride, n, w);
+  else
+    return GSL_EINVAL;  /* UNIMPLEMENTED */
+}
+
+int
+gsl_cwt_transform_forward (const gsl_cwt_wavelet *wlet, double *data,
+                          double *imag, size_t stride, size_t n,
+                          gsl_cwt_workspace *w)
+{
+  size_t i;
+  double norm, omega;
+  gsl_complex buf;
+
+  for (i = 0; i < n; i++){
+    w->packed[2*i]   = data[i];
+    w->packed[2*i+1] = 0.0;
+  }
+
+  gsl_fft_complex_forward(w->packed, stride, n, w->wavetable, w->workspace);
+
+  /* convolution with wavelet */
+  norm  = sqrt(2*M_PI * wlet->scale);
+  for (i = 0; i <= n/2; i++){
+    omega = 2*M_PI*i / n;
+    GSL_SET_COMPLEX(&buf, w->packed[2*i], w->packed[2*i+1]);
+    buf = gsl_complex_mul(buf,
+                         gsl_complex_conjugate
+                         (wlet->type->function(wlet->scale,
+                                               omega,
+                                               wlet->parms)));
+    buf = gsl_complex_mul_real(buf, norm);
+    w->packed[2*i]   = GSL_REAL(buf);
+    w->packed[2*i+1] = GSL_IMAG(buf);
+
+  }
+
+  for (i = i; i < n; i++){
+    omega = -2*M_PI*(n-i) / n;
+    GSL_SET_COMPLEX(&buf, w->packed[2*i], w->packed[2*i+1]);
+    buf = gsl_complex_mul(buf,
+                         gsl_complex_conjugate
+                         (wlet->type->function(wlet->scale,
+                                               omega,
+                                               wlet->parms)));
+    buf = gsl_complex_mul_real(buf, norm);
+    w->packed[2*i]   = GSL_REAL(buf);
+    w->packed[2*i+1] = GSL_IMAG(buf);
+
+  }
+
+
+  gsl_fft_complex_inverse(w->packed, stride, n, w->wavetable, w->workspace);
+  
+  for (i = 0; i < n; i++){
+    data[i] = w->packed[2*i];
+    imag[i] = w->packed[2*i+1];
+  }
+
+  return GSL_SUCCESS;
+}
+
+gsl_cwt_workspace *gsl_cwt_workspace_alloc (size_t n)
+{
+  gsl_cwt_workspace *cwt_workspace;
+  
+  if (!(cwt_workspace = malloc(sizeof(gsl_cwt_workspace))))
+    return NULL;
+
+  if (!(cwt_workspace->workspace = gsl_fft_complex_workspace_alloc(n))){
+    free(cwt_workspace);
+    return NULL;
+  }
+
+  if (!(cwt_workspace->wavetable = gsl_fft_complex_wavetable_alloc(n))){
+    free(cwt_workspace->workspace);
+    free(cwt_workspace);
+    return NULL;
+  }
+
+  if (!(cwt_workspace->packed = malloc(2 * n * sizeof(double)))){
+    free(cwt_workspace->workspace);
+    free(cwt_workspace->wavetable);
+    free(cwt_workspace);
+  }
+
+  return cwt_workspace;
+}
+
+void gsl_cwt_workspace_free (gsl_cwt_workspace * work)
+{
+  free(work->packed);
+  free(work->workspace);
+  free(work->wavetable);
+  free(work);
+}

pgpAhshDgrr0R.pgp
Description: PGP signature