wavelet.lisp 19.9 KB
Newer Older
1 2
;; Wavelet transforms.
;; Liam Healy, Mon Nov 26 2007 - 20:43
Liam Healy's avatar
Liam Healy committed
3
;; Time-stamp: <2009-02-16 09:59:18EST wavelet.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5 6 7

(in-package :gsl)

8
;;; /usr/include/gsl/gsl_wavelet.h
9

10
;;; Examples do not agree with C results.
11 12 13 14 15 16 17

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

;;; Utility
(defun forward-backward (symb)
18
  (if (string-equal symb :backward) -1 1))
19

20 21 22 23
;;;;****************************************************************************
;;;; Allocation of wavelets
;;;;****************************************************************************

Liam Healy's avatar
Liam Healy committed
24 25
(defmobject wavelet "gsl_wavelet"
  ((type :pointer) (member sizet))
Liam Healy's avatar
Liam Healy committed
26 27
  "wavelet"
  :documentation			; FDL
Liam Healy's avatar
Liam Healy committed
28 29 30 31
  "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.")
32

Liam Healy's avatar
Liam Healy committed
33
(defmpar +daubechies-wavelet+ "gsl_wavelet_daubechies"
34
  ;; FDL
35 36 37 38
  "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
39
(defmpar +daubechies-wavelet-centered+ "gsl_wavelet_daubechies_centered"
40
  ;; FDL
41 42 43 44
  "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
45
(defmpar +haar-wavelet+ "gsl_wavelet_haar"
46
  ;; FDL
47 48 49
  "The Haar wavelet.  The only valid choice for member for the
   Haar wavelet is member=2.")

Liam Healy's avatar
Liam Healy committed
50
(defmpar +haar-wavelet-centered+ "gsl_wavelet_haar_centered"
51
  ;; FDL
52 53 54
  "The Haar wavelet.  The only valid choice for member for the
   Haar wavelet is member=2.")

Liam Healy's avatar
Liam Healy committed
55
(defmpar +bspline-wavelet+ "gsl_wavelet_bspline"
56
  ;; FDL
57 58 59 60
  "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
61
(defmpar +bspline-wavelet-centered+ "gsl_wavelet_bspline_centered"
62
  ;; FDL
63 64 65 66
  "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.")

67
(defmfun name ((wavelet wavelet))
68
  "gsl_wavelet_name"
69 70
  (((mpointer wavelet) :pointer))
  :definition :method
71
  :c-return :string
72
  :documentation			; FDL
73 74
  "The name of the wavelet family.")

Liam Healy's avatar
Liam Healy committed
75 76 77
(defmobject wavelet-workspace
  "gsl_wavelet_workspace"
  ((size sizet))
Liam Healy's avatar
Liam Healy committed
78 79
  "wavelet workspace"
  :documentation			; FDL
Liam Healy's avatar
Liam Healy committed
80 81 82 83 84 85
  "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.")
86 87 88 89 90

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

liam's avatar
liam committed
91
(defmfun wavelet-transform (wavelet data stride direction workspace)
92
  "gsl_wavelet_transform"
93
  (((mpointer wavelet) :pointer) ((c-pointer data) :pointer)
Liam Healy's avatar
Liam Healy committed
94
   (stride sizet) ((dim0 data) sizet) ((forward-backward direction) :int)
95
   ((mpointer workspace) :pointer))
96
  :documentation			; FDL
97
  "Compute in-place forward and inverse discrete wavelet
98 99 100 101 102
   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.
103 104
   For the forward transform, the elements of the original array are 
   replaced by the discrete wavelet
105
   transform f_i -> w_{j,k}
106 107 108 109 110 111
   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}) 
112 113 114
   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 
115 116
   the original data.")

liam's avatar
liam committed
117
(defmfun wavelet-transform-forward (wavelet data stride workspace)
118
  "gsl_wavelet_transform_forward"
119 120
  (((mpointer wavelet) :pointer) ((c-pointer data) :pointer)
   (stride sizet) ((dim0 data) sizet) ((mpointer workspace) :pointer))
121
  :documentation			; FDL
122
  "Compute in-place forward and inverse discrete wavelet
123 124 125
   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.
126
   For the forward transform, the elements of the original array are 
127 128
   replaced by the discrete wavelet transform
   f_i -> w_{j,k} in a packed triangular storage layout, 
129 130 131 132 133
   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}) 
134 135 136
   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 
137 138
   the original data.")

liam's avatar
liam committed
139
(defmfun wavelet-transform-inverse (wavelet data stride workspace)
140
  "gsl_wavelet_transform_inverse"
141 142
  (((mpointer wavelet) :pointer) ((c-pointer data) :pointer)
   (stride sizet) ((dim0 data) sizet) ((mpointer workspace) :pointer))
143 144 145 146 147
  :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.
148
   For the forward transform, the elements of the original array are 
149 150
   replaced by the discrete wavelet transform
   f_i -> w_{j,k} in a packed triangular storage layout, 
151 152 153 154 155
   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}) 
156 157 158
   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 
159 160 161 162 163 164
   the original data.")

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

liam's avatar
liam committed
165
(defmfun wavelet-2d-transform (wavelet data tda direction workspace)
166
  "gsl_wavelet2d_transform"
167
  (((mpointer wavelet) :pointer) ((c-pointer data) :pointer)
Liam Healy's avatar
Liam Healy committed
168
   (tda sizet) ((dim0 data) sizet) ((dim1 data) sizet)
169
   ((forward-backward direction) :int) ((mpointer workspace) :pointer))
170 171 172 173
  :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
174
  and physical row length tda.  The dimensions must
175
  be equal (square matrix) and are restricted to powers of two.  For the
176 177 178 179
  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
180
  two-dimensional wavelet transform.
181
  An error invalid-argument is signalled if the matrix is not square
182 183 184
  with dimension a power of 2, or if insufficient
  workspace is provided.")

liam's avatar
liam committed
185
(defmfun wavelet-2d-transform-forward (wavelet data tda workspace)
186
  "gsl_wavelet2d_transform_forward"
187
  (((mpointer wavelet) :pointer) ((c-pointer data) :pointer)
Liam Healy's avatar
Liam Healy committed
188
   (tda sizet) ((dim0 data) sizet) ((dim1 data) sizet)
189
   ((mpointer workspace) :pointer))
liam's avatar
liam committed
190
  :documentation			; FDL
191 192
  "Compute two-dimensional in-place forward and inverse
  discrete wavelet transforms in standard and non-standard forms on the
liam's avatar
liam committed
193 194 195 196 197
  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
198
  two-dimensional wavelet transform.
199
  An error invalid-argument is signalled if the matrix is not square
200 201 202
  with dimension a power of 2, or if insufficient
  workspace is provided.")

liam's avatar
liam committed
203
(defmfun wavelet-2d-transform-inverse (wavelet data tda workspace)
204
  "gsl_wavelet2d_transform_inverse"
205
  (((mpointer wavelet) :pointer) ((c-pointer data) :pointer)
Liam Healy's avatar
Liam Healy committed
206
   (tda sizet) ((dim0 data) sizet) ((dim1 data) sizet)
207
   ((mpointer workspace) :pointer))
208 209 210 211 212
  :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
213 214 215
   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
216
   error invalid-argument is signalled if the matrix is not square with dimension
liam's avatar
liam committed
217 218 219
   a power of 2, or if insufficient workspace is provided.")

(defmfun wavelet-2d-transform-matrix (wavelet data direction workspace)
220
  "gsl_wavelet2d_transform_matrix"
221 222
  (((mpointer wavelet) :pointer) ((c-pointer data) :pointer)
   ((forward-backward direction) :int) ((mpointer workspace) :pointer))
223
  :documentation			; FDL
224 225
  "Compute the two-dimensional in-place wavelet transform on a matrix.")

liam's avatar
liam committed
226
(defmfun wavelet-2d-transform-matrix-forward (wavelet data workspace)
227
  "gsl_wavelet2d_transform_matrix_forward"
228 229
  (((mpointer wavelet) :pointer) ((c-pointer data) :pointer)
   ((mpointer workspace) :pointer))
230
  :documentation			; FDL
231 232
  "Compute the two-dimensional in-place wavelet transform on a matrix.")

liam's avatar
liam committed
233
(defmfun wavelet-2d-transform-matrix-inverse (wavelet data workspace)
234
  "gsl_wavelet2d_transform_matrix_inverse"
235 236
  (((mpointer wavelet) :pointer) ((c-pointer data) :pointer)
   ((mpointer workspace) :pointer))
237
  :documentation			; FDL
238 239
  "Compute the two-dimensional in-place wavelet transform on a matrix.")

liam's avatar
liam committed
240
(defmfun wavelet-2d-nonstandard-transform (wavelet data tda direction workspace)
liam's avatar
liam committed
241
  "gsl_wavelet2d_nstransform"
242
  (((mpointer wavelet) :pointer) ((c-pointer data) :pointer)
Liam Healy's avatar
Liam Healy committed
243
   (tda sizet) ((dim0 data) sizet) ((dim1 data) sizet) (direction :int)
244
   ((mpointer workspace) :pointer))
245
  :documentation			; FDL
246 247
  "Compute the two-dimensional wavelet transform in non-standard form.")

liam's avatar
liam committed
248
(defmfun wavelet-2d-nonstandard-transform-forward (wavelet data tda workspace)
liam's avatar
liam committed
249
  "gsl_wavelet2d_nstransform_forward"
250
  (((mpointer wavelet) :pointer) ((c-pointer data) :pointer)
Liam Healy's avatar
Liam Healy committed
251
   (tda sizet) ((dim0 data) sizet) ((dim1 data) sizet)
252
   ((mpointer workspace) :pointer))
253
  :documentation			; FDL
254 255
  "Compute the two-dimensional wavelet transform in non-standard form.")

liam's avatar
liam committed
256
(defmfun wavelet-2d-nonstandard-transform-inverse (wavelet data tda workspace)
liam's avatar
liam committed
257
  "gsl_wavelet2d_nstransform_inverse"
258
  (((mpointer wavelet) :pointer) ((c-pointer data) :pointer)
Liam Healy's avatar
Liam Healy committed
259
   (tda sizet) ((dim0 data) sizet) ((dim1 data) sizet)
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-matrix (wavelet data direction workspace)
liam's avatar
liam committed
265
  "gsl_wavelet2d_nstransform_matrix"
266 267
  (((mpointer wavelet) :pointer) ((c-pointer data) :pointer)
   (direction :int) ((mpointer workspace) :pointer))
268 269 270 271
  :documentation
  "Compute the non-standard form of the two-dimensional in-place wavelet
   transform on a matrix.")

liam's avatar
liam committed
272
(defmfun wavelet-2d-nonstandard-transform-matrix-forward (wavelet data workspace)
liam's avatar
liam committed
273
  "gsl_wavelet2d_nstransform_matrix_forward"
274 275
  (((mpointer wavelet) :pointer) ((c-pointer data) :pointer)
   ((mpointer workspace) :pointer))
276
  :documentation			; FDL
277 278 279
  "Compute the non-standard form of the two-dimensional in-place wavelet
   transform on a matrix.")

liam's avatar
liam committed
280
(defmfun wavelet-2d-nonstandard-transform-matrix-inverse (wavelet data workspace)
liam's avatar
liam committed
281
  "gsl_wavelet2d_nstransform_matrix_inverse"
282 283
  (((mpointer wavelet) :pointer) ((c-pointer data) :pointer)
   ((mpointer workspace) :pointer))
284
  :documentation			; FDL
285 286 287 288 289 290 291 292 293
  "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*
294
 #m(0.0462458471760794d0 0.0462458471760794d0 0.0512458471760794d0 0.0712458471760795d0
295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 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
 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.")

360 361 362
;;; These examples do not agree with their C counterparts; the answers
;;; are completely different.

363 364 365 366 367
(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."
368
  (let* ((n (length cl-data))
369
	 (vector cl-data)
Liam Healy's avatar
Liam Healy committed
370
	 (wavelet (make-wavelet +daubechies-wavelet+ 4))
371
	 (workspace (make-wavelet-workspace n)))
372
    (wavelet-transform-forward wavelet vector 1 workspace)
373 374
    (let ((absvector (make-marray 'double-float :dimensions n))
	  (permutation (make-permutation n)))
375
      (dotimes (i n)
376
	(setf (maref absvector i) (abs (maref vector i))))
377 378 379
      ;; Sort and set to 0 all but the largest 20.
      (sort-vector-index permutation absvector)
      (dotimes (i (- n 20))
380
	(setf (maref vector (maref permutation i))
381
	      0.0d0))) ;; Transform back
382
    (dotimes (i n) (format t "~&~a" (maref vector i)))
383
    (wavelet-transform-inverse wavelet vector 1 workspace)
Liam Healy's avatar
Liam Healy committed
384
    (cl-array vector)))
385 386 387

(defun wavelet-forward-example (&optional (cl-data *wavelet-sample*))
  "Simpler example, with only a Daubechies wavelet forward transformation."
388
  (let* ((n (length cl-data))
389
	 (vector cl-data)
Liam Healy's avatar
Liam Healy committed
390
	 (wavelet (make-wavelet +daubechies-wavelet+ 4))
391
	 (workspace (make-wavelet-workspace n)))
392
    (wavelet-transform-forward wavelet vector 1 workspace)
Liam Healy's avatar
Liam Healy committed
393
    (cl-array vector)))