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