wavelet.lisp 21 KB
Newer Older
1
2
;; Wavelet transforms.
;; Liam Healy, Mon Nov 26 2007 - 20:43
Liam M. Healy's avatar
Liam M. Healy committed
3
;; Time-stamp: <2014-12-26 13:18:35EST wavelet.lisp>
4
;;
Liam M. Healy's avatar
Liam M. Healy committed
5
;; Copyright 2006, 2007, 2008, 2009, 2011, 2012, 2014 Liam M. Healy
6
7
8
9
10
11
12
13
14
15
16
17
18
19
;; Distributed under the terms of the GNU General Public License
;;
;; 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 3 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, see <http://www.gnu.org/licenses/>.
20
21

(in-package :gsl)
Liam M. Healy's avatar
Liam M. Healy committed
22
(named-readtables:in-readtable :antik)
23

24
;;; /usr/include/gsl/gsl_wavelet.h
25

26
;;; Examples do not agree with C results.
27
28
29
30
31
32
33

;;;;****************************************************************************
;;;; Allocation of wavelets
;;;;****************************************************************************

;;; Utility
(defun forward-backward (symb)
34
  (if (string-equal symb :backward) -1 1))
35

36
37
38
39
;;;;****************************************************************************
;;;; Allocation of wavelets
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
40
(defmobject wavelet "gsl_wavelet"
41
  ((type :pointer) (member :sizet))
Liam Healy's avatar
Liam Healy committed
42
43
  "wavelet"
  :documentation			; FDL
Liam Healy's avatar
Liam Healy committed
44
45
46
47
  "Make and initialize a wavelet object of type 'type.  The
   parameter 'member selects the specific member of the wavelet
   family.  A memory-allocation-failure error indicates either
   lack of memory or an unsupported member requested.")
48

Liam Healy's avatar
Liam Healy committed
49
(defmpar +daubechies-wavelet+ "gsl_wavelet_daubechies"
50
  ;; FDL
51
52
53
54
  "The Daubechies wavelet family of maximum phase with member/2
   vanishing moments.  The implemented wavelets are 
   member=4, 6,..., 20, with member even.")

Liam Healy's avatar
Liam Healy committed
55
(defmpar +daubechies-wavelet-centered+ "gsl_wavelet_daubechies_centered"
56
  ;; FDL
57
58
59
60
  "The Daubechies wavelet family of maximum phase with member/2
   vanishing moments.  The implemented wavelets are 
   member=4, 6,..., 20, with member even.")

Liam Healy's avatar
Liam Healy committed
61
(defmpar +haar-wavelet+ "gsl_wavelet_haar"
62
  ;; FDL
63
64
65
  "The Haar wavelet.  The only valid choice for member for the
   Haar wavelet is member=2.")

Liam Healy's avatar
Liam Healy committed
66
(defmpar +haar-wavelet-centered+ "gsl_wavelet_haar_centered"
67
  ;; FDL
68
69
70
  "The Haar wavelet.  The only valid choice for member for the
   Haar wavelet is member=2.")

Liam Healy's avatar
Liam Healy committed
71
(defmpar +bspline-wavelet+ "gsl_wavelet_bspline"
72
  ;; FDL
73
74
75
76
  "The biorthogonal B-spline wavelet family of order (i,j).  
   The implemented values of member = 100*i + j are 103, 105, 202, 204,
   206, 208, 301, 303, 305 307, 309.")

Liam Healy's avatar
Liam Healy committed
77
(defmpar +bspline-wavelet-centered+ "gsl_wavelet_bspline_centered"
78
  ;; FDL
79
80
81
82
  "The biorthogonal B-spline wavelet family of order (i,j).  
   The implemented values of member = 100*i + j are 103, 105, 202, 204,
   206, 208, 301, 303, 305 307, 309.")

83
(defmfun name ((wavelet wavelet))
84
  "gsl_wavelet_name"
85
86
  (((mpointer wavelet) :pointer))
  :definition :method
87
  :c-return :string
88
  :documentation			; FDL
89
90
  "The name of the wavelet family.")

Liam Healy's avatar
Liam Healy committed
91
92
(defmobject wavelet-workspace
  "gsl_wavelet_workspace"
93
  ((size :sizet))
Liam Healy's avatar
Liam Healy committed
94
95
  "wavelet workspace"
  :documentation			; FDL
Liam Healy's avatar
Liam Healy committed
96
97
98
99
100
101
  "Make a workspace for the discrete wavelet transform.
   To perform a one-dimensional transform on size elements, a workspace
   of size size must be provided.  For two-dimensional transforms of
   size-by-size matrices it is sufficient to allocate a workspace of
   size, since the transform operates on individual rows and
   columns.")
102
103
104
105
106

;;;;****************************************************************************
;;;; Wavelet transforms 1D
;;;;****************************************************************************

liam's avatar
liam committed
107
(defmfun wavelet-transform (wavelet data stride direction workspace)
108
  "gsl_wavelet_transform"
109
  (((mpointer wavelet) :pointer) ((grid:foreign-pointer data) :pointer)
110
   (stride :sizet) ((dim0 data) :sizet) ((forward-backward direction) :int)
111
   ((mpointer workspace) :pointer))
112
  :documentation			; FDL
113
  "Compute in-place forward and inverse discrete wavelet
114
115
116
117
118
   transforms on the array data. The length of the
   transform n is restricted to powers
   of two.  For the transform version of the function the argument
   dir can be either :forward or :backward.  A workspace
   of the same length as data must be provided.
119
120
   For the forward transform, the elements of the original array are 
   replaced by the discrete wavelet
121
   transform f_i -> w_{j,k}
122
123
124
125
126
127
   in a packed triangular storage layout, 
   where j is the index of the level j = 0 ... J-1
   and k is the index of the coefficient within each level,
   k = 0 ... (2^j)-1.  The total number of levels is J = \log_2(n).
   The output data has the following form,
   (s_{-1,0}, d_{0,0}, d_{1,0}, d_{1,1}, d_{2,0},\cdots, d_{j,k},\cdots, d_{J-1,2^{J-1} - 1}) 
128
129
130
   where the first element is the smoothing coefficient
   s_{-1,0}, followed by the detail coefficients d_{j,k} for each level j.
   The backward transform inverts these coefficients to obtain 
131
132
   the original data.")

liam's avatar
liam committed
133
(defmfun wavelet-transform-forward (wavelet data stride workspace)
134
  "gsl_wavelet_transform_forward"
135
  (((mpointer wavelet) :pointer) ((grid:foreign-pointer data) :pointer)
136
   (stride :sizet) ((dim0 data) :sizet) ((mpointer workspace) :pointer))
137
  :documentation			; FDL
138
  "Compute in-place forward and inverse discrete wavelet
139
140
141
   transforms on the array data.  The length of the transform
   is restricted to powers of two.
   A workspace of the same length as data must be provided.
142
   For the forward transform, the elements of the original array are 
143
144
   replaced by the discrete wavelet transform
   f_i -> w_{j,k} in a packed triangular storage layout, 
145
146
147
148
149
   where j is the index of the level j = 0 ... J-1
   and k is the index of the coefficient within each level,
   k = 0 ... (2^j)-1.  The total number of levels is J = \log_2(n).
   The output data has the following form,
   (s_{-1,0}, d_{0,0}, d_{1,0}, d_{1,1}, d_{2,0},\cdots, d_{j,k},\cdots, d_{J-1,2^{J-1} - 1}) 
150
151
152
   where the first element is the smoothing coefficient s_{-1,0},
   followed by the detail coefficients d_{j,k} for each level j.
   The backward transform inverts these coefficients to obtain 
153
154
   the original data.")

liam's avatar
liam committed
155
(defmfun wavelet-transform-inverse (wavelet data stride workspace)
156
  "gsl_wavelet_transform_inverse"
157
  (((mpointer wavelet) :pointer) ((grid:foreign-pointer data) :pointer)
158
   (stride :sizet) ((dim0 data) :sizet) ((mpointer workspace) :pointer))
159
160
161
162
163
  :documentation			; FDL
  "Compute in-place inverse discrete wavelet
   transforms on the array data.  The length of the transform
   is restricted to powers of two.
   A workspace of the same length as data must be provided.
164
   For the forward transform, the elements of the original array are 
165
166
   replaced by the discrete wavelet transform
   f_i -> w_{j,k} in a packed triangular storage layout, 
167
168
169
170
171
   where j is the index of the level j = 0 ... J-1
   and k is the index of the coefficient within each level,
   k = 0 ... (2^j)-1.  The total number of levels is J = \log_2(n).
   The output data has the following form,
   (s_{-1,0}, d_{0,0}, d_{1,0}, d_{1,1}, d_{2,0},\cdots, d_{j,k},\cdots, d_{J-1,2^{J-1} - 1}) 
172
173
174
   where the first element is the smoothing coefficient s_{-1,0},
   followed by the detail coefficients d_{j,k} for each level j.
   The backward transform inverts these coefficients to obtain 
175
176
177
178
179
180
   the original data.")

;;;;****************************************************************************
;;;; Wavelet transforms 2D
;;;;****************************************************************************

liam's avatar
liam committed
181
(defmfun wavelet-2d-transform (wavelet data tda direction workspace)
182
  "gsl_wavelet2d_transform"
183
  (((mpointer wavelet) :pointer) ((grid:foreign-pointer data) :pointer)
184
   (tda :sizet) ((dim0 data) :sizet) ((dim1 data) :sizet)
185
   ((forward-backward direction) :int) ((mpointer workspace) :pointer))
186
187
188
189
  :documentation			; FDL
  "Compute in-place forward and inverse discrete wavelet transforms
  in standard and non-standard forms on the
  array data stored in row-major form with dimensions
liam's avatar
liam committed
190
  and physical row length tda.  The dimensions must
191
  be equal (square matrix) and are restricted to powers of two.  For the
192
193
194
195
  transform version of the function the argument direction can be
  either :forward or :backward.  A
  workspace of the appropriate size must be provided.  On exit,
  the appropriate elements of the array data are replaced by their
196
  two-dimensional wavelet transform.
197
  An error invalid-argument is signalled if the matrix is not square
198
199
200
  with dimension a power of 2, or if insufficient
  workspace is provided.")

liam's avatar
liam committed
201
(defmfun wavelet-2d-transform-forward (wavelet data tda workspace)
202
  "gsl_wavelet2d_transform_forward"
203
  (((mpointer wavelet) :pointer) ((grid:foreign-pointer data) :pointer)
204
   (tda :sizet) ((dim0 data) :sizet) ((dim1 data) :sizet)
205
   ((mpointer workspace) :pointer))
liam's avatar
liam committed
206
  :documentation			; FDL
207
208
  "Compute two-dimensional in-place forward and inverse
  discrete wavelet transforms in standard and non-standard forms on the
liam's avatar
liam committed
209
210
211
212
213
  array data stored in row-major form with dimensions size1
  and size2 and physical row length tda.  The dimensions must
  be equal (square matrix) and are restricted to powers of two.  A
  workspace of the appropriate size must be provided.  On exit,
  the appropriate elements of the array data are replaced by their
214
  two-dimensional wavelet transform.
215
  An error invalid-argument is signalled if the matrix is not square
216
217
218
  with dimension a power of 2, or if insufficient
  workspace is provided.")

liam's avatar
liam committed
219
(defmfun wavelet-2d-transform-inverse (wavelet data tda workspace)
220
  "gsl_wavelet2d_transform_inverse"
221
  (((mpointer wavelet) :pointer) ((grid:foreign-pointer data) :pointer)
222
   (tda :sizet) ((dim0 data) :sizet) ((dim1 data) :sizet)
223
   ((mpointer workspace) :pointer))
224
225
226
227
228
  :documentation			; FDL
  "Compute two-dimensional in-place forward and inverse discrete
   wavelet transforms in standard and non-standard forms on the array
   data stored in row-major form with dimensions size1 and size2 and
   physical row length tda.  The dimensions must be equal (square matrix)
liam's avatar
liam committed
229
230
231
   and are restricted to powers of two.  A workspace of the appropriate
   size must be provided.  On exit, the appropriate elements of the array
   data are replaced by their two-dimensional wavelet transform.  An
232
   error invalid-argument is signalled if the matrix is not square with dimension
liam's avatar
liam committed
233
234
235
   a power of 2, or if insufficient workspace is provided.")

(defmfun wavelet-2d-transform-matrix (wavelet data direction workspace)
236
  "gsl_wavelet2d_transform_matrix"
237
  (((mpointer wavelet) :pointer) ((grid:foreign-pointer data) :pointer)
238
   ((forward-backward direction) :int) ((mpointer workspace) :pointer))
239
  :documentation			; FDL
240
241
  "Compute the two-dimensional in-place wavelet transform on a matrix.")

liam's avatar
liam committed
242
(defmfun wavelet-2d-transform-matrix-forward (wavelet data workspace)
243
  "gsl_wavelet2d_transform_matrix_forward"
244
  (((mpointer wavelet) :pointer) ((grid:foreign-pointer data) :pointer)
245
   ((mpointer workspace) :pointer))
246
  :documentation			; FDL
247
248
  "Compute the two-dimensional in-place wavelet transform on a matrix.")

liam's avatar
liam committed
249
(defmfun wavelet-2d-transform-matrix-inverse (wavelet data workspace)
250
  "gsl_wavelet2d_transform_matrix_inverse"
251
  (((mpointer wavelet) :pointer) ((grid:foreign-pointer data) :pointer)
252
   ((mpointer workspace) :pointer))
253
  :documentation			; FDL
254
255
  "Compute the two-dimensional in-place wavelet transform on a matrix.")

liam's avatar
liam committed
256
(defmfun wavelet-2d-nonstandard-transform (wavelet data tda direction workspace)
liam's avatar
liam committed
257
  "gsl_wavelet2d_nstransform"
258
  (((mpointer wavelet) :pointer) ((grid:foreign-pointer data) :pointer)
259
   (tda :sizet) ((dim0 data) :sizet) ((dim1 data) :sizet) (direction :int)
260
   ((mpointer workspace) :pointer))
261
  :documentation			; FDL
262
263
  "Compute the two-dimensional wavelet transform in non-standard form.")

liam's avatar
liam committed
264
(defmfun wavelet-2d-nonstandard-transform-forward (wavelet data tda workspace)
liam's avatar
liam committed
265
  "gsl_wavelet2d_nstransform_forward"
266
  (((mpointer wavelet) :pointer) ((grid:foreign-pointer data) :pointer)
267
   (tda :sizet) ((dim0 data) :sizet) ((dim1 data) :sizet)
268
   ((mpointer workspace) :pointer))
269
  :documentation			; FDL
270
271
  "Compute the two-dimensional wavelet transform in non-standard form.")

liam's avatar
liam committed
272
(defmfun wavelet-2d-nonstandard-transform-inverse (wavelet data tda workspace)
liam's avatar
liam committed
273
  "gsl_wavelet2d_nstransform_inverse"
274
  (((mpointer wavelet) :pointer) ((grid:foreign-pointer data) :pointer)
275
   (tda :sizet) ((dim0 data) :sizet) ((dim1 data) :sizet)
276
   ((mpointer workspace) :pointer))
277
  :documentation			; FDL
278
279
  "Compute the two-dimensional wavelet transform in non-standard form.")

liam's avatar
liam committed
280
(defmfun wavelet-2d-nonstandard-transform-matrix (wavelet data direction workspace)
liam's avatar
liam committed
281
  "gsl_wavelet2d_nstransform_matrix"
282
  (((mpointer wavelet) :pointer) ((grid:foreign-pointer data) :pointer)
283
   (direction :int) ((mpointer workspace) :pointer))
284
285
286
287
  :documentation
  "Compute the non-standard form of the two-dimensional in-place wavelet
   transform on a matrix.")

liam's avatar
liam committed
288
(defmfun wavelet-2d-nonstandard-transform-matrix-forward (wavelet data workspace)
liam's avatar
liam committed
289
  "gsl_wavelet2d_nstransform_matrix_forward"
290
  (((mpointer wavelet) :pointer) ((grid:foreign-pointer data) :pointer)
291
   ((mpointer workspace) :pointer))
292
  :documentation			; FDL
293
294
295
  "Compute the non-standard form of the two-dimensional in-place wavelet
   transform on a matrix.")

liam's avatar
liam committed
296
(defmfun wavelet-2d-nonstandard-transform-matrix-inverse (wavelet data workspace)
liam's avatar
liam committed
297
  "gsl_wavelet2d_nstransform_matrix_inverse"
298
  (((mpointer wavelet) :pointer) ((grid:foreign-pointer data) :pointer)
299
   ((mpointer workspace) :pointer))
300
  :documentation			; FDL
301
302
303
304
305
306
307
308
309
  "Compute the non-standard form of the two-dimensional in-place wavelet
   transform on a matrix.")

;;;;****************************************************************************
;;;; Example
;;;;****************************************************************************

;;; See GSL manual Section 30.4.
(defparameter *wavelet-sample*
310
 #m(0.0462458471760794d0 0.0462458471760794d0 0.0512458471760794d0 0.0712458471760795d0
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
 0.0712458471760795d0 0.0662458471760795d0 0.0962458471760795d0 0.101245847176079d0
 0.116245847176079d0 0.121245847176079d0 0.116245847176079d0 0.106245847176079d0
 0.0912458471760794d0 0.101245847176079d0 0.0962458471760795d0 0.0962458471760795d0
 0.0962458471760795d0 0.0912458471760794d0 0.0862458471760795d0 0.0812458471760795d0
 0.0862458471760795d0 0.101245847176079d0 0.111245847176079d0 0.116245847176079d0
 0.0762458471760795d0 0.0362458471760795d0 0.0362458471760795d0 0.0212458471760795d0
 0.0112458471760795d0 -0.00875415282392056d0 -0.00875415282392056d0 -0.00375415282392055d0
 0.00624584717607946d0 0.00124584717607945d0 0.00624584717607946d0 -0.00375415282392055d0
 -0.0187541528239206d0 -0.0237541528239205d0 -0.0187541528239206d0 -0.0187541528239206d0
 -0.0287541528239205d0 -0.0237541528239205d0 -0.0337541528239205d0 -0.00875415282392056d0
 -0.0137541528239206d0 -0.00875415282392056d0 0.00124584717607945d0 -0.0237541528239205d0
 -0.0337541528239205d0 -0.0187541528239206d0 -0.00875415282392056d0 -0.00375415282392055d0
 -0.00875415282392056d0 -0.0287541528239205d0 -0.0437541528239205d0 -0.0387541528239205d0
 -0.0587541528239205d0 -0.103754152823921d0 -0.123754152823921d0 -0.153754152823921d0
 -0.188754152823921d0 -0.213754152823921d0 -0.183754152823921d0 -0.0937541528239205d0
 0.0212458471760795d0 0.161245847176079d0 0.306245847176079d0 0.556245847176079d0
 0.81124584717608d0 1.04124584717608d0 1.19624584717608d0 1.26124584717608d0
 1.22624584717608d01 1.07624584717608d0 0.81124584717608d0 0.486245847176079d0
 0.211245847176079d0 0.0512458471760794d0 -0.0687541528239206d0 -0.128754152823921d0
 -0.153754152823921d0 -0.133754152823921d0 -0.103754152823921d0 -0.0687541528239206d0
 -0.0687541528239206d0 -0.0637541528239206d0 -0.0687541528239206d0 -0.0587541528239205d0
 -0.0587541528239205d0 -0.0587541528239205d0 -0.0737541528239206d0 -0.0637541528239206d0
 -0.0637541528239206d0 -0.0637541528239206d0 -0.0537541528239205d0 -0.0737541528239206d0
 -0.0887541528239205d0 -0.0887541528239205d0 -0.0787541528239206d0 -0.0737541528239206d0
 -0.0687541528239206d0 -0.0837541528239206d0 -0.0737541528239206d0 -0.0637541528239206d0
 -0.0537541528239205d0 -0.0687541528239206d0 -0.0687541528239206d0 -0.0837541528239206d0
 -0.0887541528239205d0 -0.0887541528239205d0 -0.0687541528239206d0 -0.0687541528239206d0
 -0.0737541528239206d0 -0.0837541528239206d0 -0.0937541528239205d0 -0.0787541528239206d0
 -0.0887541528239205d0 -0.0837541528239206d0 -0.0887541528239205d0 -0.0937541528239205d0
 -0.0887541528239205d0 -0.0787541528239206d0 -0.0787541528239206d0 -0.0737541528239206d0
 -0.0687541528239206d0 -0.0837541528239206d0 -0.0887541528239205d0 -0.0687541528239206d0
 -0.0687541528239206d0 -0.0637541528239206d0 -0.0637541528239206d0 -0.0887541528239205d0
 -0.0837541528239206d0 -0.0737541528239206d0 -0.0687541528239206d0 -0.0537541528239205d0
 -0.0687541528239206d0 -0.0737541528239206d0 -0.0887541528239205d0 -0.0787541528239206d0
 -0.0687541528239206d0 -0.0687541528239206d0 -0.0637541528239206d0 -0.0837541528239206d0
 -0.0937541528239205d0 -0.0937541528239205d0 -0.0787541528239206d0 -0.0737541528239206d0
 -0.0837541528239206d0 -0.0937541528239205d0 -0.0987541528239205d0 -0.0987541528239205d0
 -0.0887541528239205d0 -0.0937541528239205d0 -0.103754152823921d0 -0.0987541528239205d0
 -0.113754152823921d0 -0.108754152823921d0 -0.108754152823921d0 -0.0987541528239205d0
 -0.108754152823921d0 -0.128754152823921d0 -0.133754152823921d0 -0.128754152823921d0
 -0.113754152823921d0 -0.123754152823921d0 -0.128754152823921d0 -0.133754152823921d0
 -0.148754152823921d0 -0.138754152823921d0 -0.133754152823921d0 -0.128754152823921d0
 -0.133754152823921d0 -0.148754152823921d0 -0.153754152823921d0 -0.138754152823921d0
 -0.128754152823921d0 -0.123754152823921d0 -0.118754152823921d0 -0.113754152823921d0
 -0.118754152823921d0 -0.0887541528239205d0 -0.0737541528239206d0 -0.0487541528239205d0
 -0.0437541528239205d0 -0.0387541528239205d0 -0.0437541528239205d0 -0.0187541528239206d0
 -0.00375415282392055d0 0.00624584717607946d0 0.00124584717607945d0 -0.00875415282392056d0
 -0.00875415282392056d0 0.00124584717607945d0 0.0112458471760795d0 0.0212458471760795d0
 0.0212458471760795d0 0.00124584717607945d0 0.00124584717607945d0 0.00624584717607946d0
 0.0162458471760795d0 0.0162458471760795d0 0.0262458471760795d0 0.00124584717607945d0
 -0.00875415282392056d0 0.0162458471760795d0 0.0112458471760795d0 0.0212458471760795d0
 0.0212458471760795d0 0.00124584717607945d0 -0.00375415282392055d0 0.0112458471760795d0
 0.0162458471760795d0 0.00624584717607946d0 0.0162458471760795d0 0.00624584717607946d0
 0.00624584717607946d0 0.0112458471760795d0 0.0262458471760795d0 0.0312458471760795d0
 0.0162458471760795d0 0.0112458471760795d0 0.00124584717607945d0 0.00624584717607946d0
 0.0212458471760795d0 0.00624584717607946d0 0.00624584717607946d0 0.00624584717607946d0
 -0.00875415282392056d0 0.00624584717607946d0 0.00124584717607945d0 0.00624584717607946d0
 -0.00375415282392055d0 -0.0137541528239206d0 -0.0187541528239206d0 -0.0137541528239206d0
 -0.0137541528239206d0 -0.00875415282392056d0 -0.00375415282392055d0 -0.0237541528239205d0
 -0.0287541528239205d0 -0.0237541528239205d0 -0.0137541528239206d0 -0.00875415282392056d0
 -0.00875415282392056d0 -0.0237541528239205d0 -0.0237541528239205d0 -0.0237541528239205d0
 0.00124584717607945d0 -0.00875415282392056d0 -0.0137541528239206d0 -0.0187541528239206d0
 -0.0337541528239205d0 -0.0137541528239206d0 -0.00875415282392056d0 -0.00875415282392056d0)
  "Data for example wavelet transform from doc/examples/ecg.dat.")

376
377
378
;;; These examples do not agree with their C counterparts; the answers
;;; are completely different.

379
380
381
382
383
(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
   256) using the 20 largest components of the wavelet transform, while
   setting the others to zero.  See GSL manual Section 30.4."
384
  (let* ((n (length cl-data))
385
	 (vector cl-data)
Liam Healy's avatar
Liam Healy committed
386
	 (wavelet (make-wavelet +daubechies-wavelet+ 4))
387
	 (workspace (make-wavelet-workspace n)))
388
    (wavelet-transform-forward wavelet vector 1 workspace)
389
    (let ((absvector (grid:make-foreign-array 'double-float :dimensions n))
390
	  (permutation (make-permutation n)))
391
      (dotimes (i n)
392
	(setf (grid:aref absvector i) (abs (grid:aref vector i))))
393
394
395
      ;; Sort and set to 0 all but the largest 20.
      (sort-vector-index permutation absvector)
      (dotimes (i (- n 20))
396
	(setf (grid:aref vector (grid:aref permutation i))
397
	      0.0d0))) ;; Transform back
398
    (dotimes (i n) (format t "~&~a" (grid:aref vector i)))
399
    (wavelet-transform-inverse wavelet vector 1 workspace)
400
    (grid:copy-to vector)))
401
402
403

(defun wavelet-forward-example (&optional (cl-data *wavelet-sample*))
  "Simpler example, with only a Daubechies wavelet forward transformation."
404
  (let* ((n (length cl-data))
405
	 (vector cl-data)
Liam Healy's avatar
Liam Healy committed
406
	 (wavelet (make-wavelet +daubechies-wavelet+ 4))
407
	 (workspace (make-wavelet-workspace n)))
408
    (wavelet-transform-forward wavelet vector 1 workspace)
409
    (grid:copy-to vector)))