Commit cc7f1b72 authored by liam's avatar liam
Browse files

Finish up wavelet example, which does not agree with the C

counterpart, and implement Hankel transforms.


git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3254 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent 878682c3
......@@ -3,7 +3,7 @@
; description: Definition of GSLL system
; date:
; author: Liam Healy
; modified: Mon Nov 26 2007 - 23:33
; modified: Sat Dec 8 2007 - 18:41
;********************************************************
;;; $Id: $
......@@ -173,4 +173,5 @@
(:file "numerical-differentiation" :depends-on (init numerical-integration))
(:file "chebyshev" :depends-on (init numerical-integration))
(:file "series-acceleration" :depends-on (init))
(:file "wavelet" :depends-on (init data))))
(:file "wavelet" :depends-on (init data))
(:file "hankel" :depends-on (init data))))
;********************************************************
; file: hankel.lisp
; description: Discrete Hankel Transforms.
; date: Sat Dec 8 2007 - 16:50
; author: Liam Healy
; modified: Sat Dec 8 2007 - 18:38
;********************************************************
;;; $Id: $
(in-package :gsl)
;;; Everything compiles, but not tested -- need example.
(defun-gsl allocate-hankel (size)
"gsl_dht_alloc"
((size :size))
:c-return :pointer
:documentation
"Allocate a Discrete Hankel transform object of size @var{size}.")
(defun-gsl init-hankel (hankel nu xmax)
"gsl_dht_new"
((hankel :pointer) (nu :double) (xmax :double))
:documentation
"Initialize the transform @var{t} for the given values of
@var{nu} and @var{x}.")
(defun-gsl new-hankel (size nu xmax)
"gsl_dht_new"
((size :size) (nu :double) (xmax :double))
:c-return :pointer
:documentation
"Allocate a Discrete Hankel transform object of size
@var{size} and initializes it for the given values of @var{nu} and
@var{x}.")
(defun-gsl free-hankel (hankel)
"gsl_dht_free"
((hankel :pointer))
:c-return :void
:documentation
"Free the Hankel object.")
(defun-gsl apply-hankel (hankel array-in array-out)
"gsl_dht_apply"
((hankel :pointer) ((gsl-array array-in) :pointer)
((gsl-array array-out) :pointer))
:documentation
"Apply the transform @var{t} to the array @var{f_in}
whose size is equal to the size of the transform. The result is stored
in the array @var{f_out} which must be of the same length.")
(defun-gsl sample-x-hankel (hankel n)
"gsl_dht_x_sample"
((hankel :pointer) (n :int))
:c-return :double
:documentation
"The points where the function @math{f(t)} is assumed to be sampled.")
(defun-gsl sample-k-hankel (hankel n)
"gsl_dht_k_sample"
((hankel :pointer) (n :int))
:c-return :double
:documentation
"The value of the @var{n}-th sample point in k-space,
j_{\nu,n+1}/X}.")
......@@ -298,11 +298,11 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>
</tr>
<tr>
<td><a href="Wavelet-Transforms.html">Wavelet Transforms</a></td>
<td>Not started</td>
<td>Done but example disagrees with C</td>
</tr>
<tr>
<td><a href="Discrete-Hankel-Transforms.html">Discrete Hankel Transforms</a></td>
<td>Not started</td>
<td>Done but untested (no example available)</td>
</tr>
<tr>
<td><a
......@@ -358,7 +358,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>
<address><a href="mailto:">Liam Healy</a></address>
<!-- Created: Feb 25 2005 -->
<!-- hhmts start -->
Last modified: April 11, 2007
Last modified: December 9, 2007
<!-- hhmts end -->
</div>
......
......@@ -3,12 +3,14 @@
; description: Wavelet transforms.
; date: Mon Nov 26 2007 - 20:43
; author: Liam Healy
; modified: Mon Nov 26 2007 - 23:26
; modified: Sun Dec 2 2007 - 17:10
;********************************************************
;;; $Id: $
(in-package :gsl)
;;; Examples do not agree with C results.
;;;;****************************************************************************
;;;; Allocation of wavelets
;;;;****************************************************************************
......@@ -364,6 +366,9 @@
-0.0337541528239205d0 -0.0137541528239206d0 -0.00875415282392056d0 -0.00875415282392056d0)
"Data for example wavelet transform from doc/examples/ecg.dat.")
;;; These examples do not agree with their C counterparts; the answers
;;; are completely different.
(defun wavelet-example (&optional (cl-data *wavelet-sample*))
"Demonstrates the use of the one-dimensional wavelet transform
functions. It computes an approximation to an input signal (of length
......@@ -375,14 +380,32 @@
(let ((wavelet (allocate-wavelet *daubechies-wavelet* 4))
(workspace (allocate-wavelet-workspace n)))
(wavelet-transform-forward wavelet vector 1 workspace)
;; Sort and set to 0 all but the largest 20.
(with-data (permutation permutation n)
(sort-vector-index permutation vector)
(loop for i from 20 below n do
(setf (gsl-aref vector (gsl-aref permutation i))
0.0d0))) ;; Transform back
(with-data (absvector vector-double n)
(dotimes (i n)
(setf (gsl-aref absvector i) (abs (gsl-aref vector i))))
;; Sort and set to 0 all but the largest 20.
(with-data (permutation permutation n)
(sort-vector-index permutation absvector)
(dotimes (i (- n 20))
(setf (gsl-aref vector (gsl-aref permutation i))
0.0d0)))) ;; Transform back
(dotimes (i n) (format t "~&~a" (gsl-aref vector i)))
(cl-invalidate vector)
(wavelet-transform-inverse wavelet vector 1 workspace)
(prog1 (data vector)
(free-wavelet-workspace workspace)
(free-wavelet wavelet))))))
(defun wavelet-forward-example (&optional (cl-data *wavelet-sample*))
"Simpler example, with only a Daubechies wavelet forward transformation."
(let ((n (length cl-data)))
(with-data (vector vector-double n)
(setf (data vector) cl-data)
(let ((wavelet (allocate-wavelet *daubechies-wavelet* 4))
(workspace (allocate-wavelet-workspace n)))
(wavelet-transform-forward wavelet vector 1 workspace)
(cl-invalidate vector)
(prog1 (data vector)
(free-wavelet-workspace workspace)
(free-wavelet wavelet))))))
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment