bug-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Bug-gsl] [Patch 0/4] Introduce continuous wavelet transform


From: Kevin Bube
Subject: Re: [Bug-gsl] [Patch 0/4] Introduce continuous wavelet transform
Date: Mon, 02 May 2005 16:23:57 +0200
User-agent: Gnus/5.1006 (Gnus v5.10.6) Emacs/21.3 (gnu/linux)

Hi all,

during the weekend I had second glance at the patches and I think now it
is a better idea to change the cwt functions to accept complex data
instead of only real ones. The following two patches do that. Furthermore
the functions require less memory and should run a bit faster, as less
data is copied.

Regards,

Kevin

--- gsl-1.6/wavelet/gsl_cwt.old 2005-04-29 17:24:42.000000000 +0200
+++ gsl-1.6/wavelet/gsl_cwt.h   2005-05-02 09:54:30.000000000 +0200
@@ -22,6 +22,7 @@
 #include <gsl/gsl_types.h>
 #include <gsl/gsl_errno.h>
 #include <gsl/gsl_fft_complex.h>
+#include <gsl/gsl_complex.h>
 
 #undef __BEGIN_DECLS
 #undef __END_DECLS
@@ -59,7 +60,6 @@
 {
   gsl_fft_complex_wavetable *wavetable;
   gsl_fft_complex_workspace *workspace;
-  double *packed;
   size_t n;
 }
 gsl_cwt_workspace;
@@ -75,13 +75,13 @@
 extern void gsl_cwt_workspace_free (gsl_cwt_workspace * work);
 
 extern int gsl_cwt_transform (const gsl_cwt_wavelet *wlet,
-                             double *data, double *imag,
+                             gsl_complex *data,
                              size_t stride, size_t n,
                              gsl_cwt_direction dir,
                              gsl_cwt_workspace *w);
 
 extern int gsl_cwt_transform_forward (const gsl_cwt_wavelet *wlet,
-                                     double *data, double *imag,
+                                     gsl_complex *data,
                                      size_t stride, size_t n,
                                      gsl_cwt_workspace *w);
 
--- gsl-1.6/wavelet/cwt.old     2005-04-29 17:24:42.000000000 +0200
+++ gsl-1.6/wavelet/cwt.c       2005-05-02 09:58:41.000000000 +0200
@@ -32,71 +32,53 @@
 #define ELEMENT(a,stride,i) ((a)[(stride)*(i)])
 
 int
-gsl_cwt_transform (const gsl_cwt_wavelet *wlet, double *data, double *imag,
+gsl_cwt_transform (const gsl_cwt_wavelet *wlet, gsl_complex *data,
                   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);
+    return gsl_cwt_transform_forward (wlet, data, 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)
+gsl_cwt_transform_forward (const gsl_cwt_wavelet *wlet, gsl_complex *data,
+                          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);
+  gsl_fft_complex_forward((gsl_complex_packed_array) data, 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,
+    data[i] = gsl_complex_mul(data[i],
+                             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);
-
+    data[i] = gsl_complex_mul_real(data[i], norm);
   }
 
   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);
-
+    data[i] = gsl_complex_mul(data[i],
+                             gsl_complex_conjugate
+                             (wlet->type->function(wlet->scale,
+                                                   omega,
+                                                   wlet->parms)));
+    data[i] = gsl_complex_mul_real(data[i], norm);
   }
 
 
-  gsl_fft_complex_inverse(w->packed, stride, n, w->wavetable, w->workspace);
+  gsl_fft_complex_inverse((gsl_complex_packed_array) data, 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;
 }
 
@@ -117,19 +99,11 @@
     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);

Attachment: pgpAGE4JZFBbL.pgp
Description: PGP signature


reply via email to

[Prev in Thread] Current Thread [Next in Thread]