revised version of wavelet transform

[Ivo Alxneit]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

On Sunday 25 April 2004 22:21, you wrote:
> I've had a look at this now.  It is a good improvement.
>
> A couple of points I noticed,
>
> - in the initialisation of h1, h2, g1, g2 in daubechies and bspline
> there are non-constant static arrays, which are not thread-safe.  To
> avoid this, I'd recommend supplying all the coefficients as
> precomputed tables (even though they only differ in sign).

ok, i'll do this.

>
> - there could be some value in putting the allowed wavelet lengths in
> an array in the wavelet_type struct, to allow a calling routine to
> inspect them. The array could be terminated with 0.

let's see.
you mean something like a call gsl_allowed_members(gsl_wavelet_daubechies) 
returning an array [4, 6, 8 , 10, 12, 14, 16, 18, 20, 0] so that a user can 
check if a certain member is available or not before trying to allocate the 
structure.

i dont think that this is really usefull. maybe your thought comes from the 
fact that the "size_t n" in gsl_wavelet_alloc() is really the member 
selected. this CAN BE the number of coefficients (haar, daubechies) or just 
an acronym (bspline, 204 means 2,4 spline of the original publication). i did 
implement all members of those families that were given in the original 
publication. so it is not like e.g. in some special functions where the 
implementation only works for a certain sub-range of an argument. also you 
usually have good reasons to select a certain member of a family, so you know 
beforehand what you need. for the daubechies "the user knows" that you select 
the member with 2N coefficients to apply a wavelet with N non-vanishing 
moments. of course the documentation will state what members of a certain 
family are implemented.

if you still consider your idea usefull i will implement it, no problem. 
(maybe i did not understand your point completely)

one more rfc from my side. please have a look at the three postscript files in 
files.tar.gz.
ping.ps (the signal),
h_ping.ps (heisenberg box plot of the wavelet transform) and
h_ping-c.ps (heisenberg box plot of the wavelet transform, wavelets are 
centered).

heisenberg box plot: "signal axis" (x or t) horizontal, running left to right. 
"frequency axis" (cm-1, Hz) vertical, low frequency at the bottom high 
frequency at the top). the magnitude of the coefficients are plotted (dark 
values equals large coefficient). so it tells you what frequency component 
appears at which location in the signal. 
in the case of an fft the heisenberg plot consists of stripes of equal 
dimension as the frequency range is equally spaved and each component spans 
the whole signal-axis (no localization in space for fft). sliceing this plot 
vertically at some place gives the spectrum of the input signal.

for a dwt this is different because the frequancy range is logarithmic 
(doubling the frequency going from one band to the next) and the wavelets are 
localized. you will notice that the frequency bands in h_ping-c.ps nicely 
line up. This illustrates that the ping signal is only "on" during that time 
and that it contains some contribution for all frequency bands. in h_ping.ps, 
however, the bands do not line up.

the current implementation results in h_ping.ps because the wavelets are not 
centered i.e. the position on the wavelet is associated with its trailing 
(left) edge. this is the way all implementations i checked work. if you 
center the wavelet i.e. its position is associated with its center and you 
will obtain the easier-to-read heisenberg box plot (h_ping-c.ps).
leaving the choice (centered/not centered) to the user is easy to implement:
1) add a boolean to the gsl_wavelet structure.
2) allocate according to gsl_wavelet_allocate (family, member, 
centered/not-centered).
3) add if(centered) nmod -= (w->nc >>1) in gsl_dwt_step().

or (instead of 1 and 2 above) add a parameter to the different calls 
gsl_dwt_....(..., centered/not centered) 

should i add this ? i prefer the second solution

p.s. there is a typo in dwt.c. i misnamed the standart/non-standart 2D 
transform. gsl_dwt_2D should be named gsl_dwt_2Dns and vice versa !
sorry  :-(
- -- 
Dr. Ivo Alxneit
Laboratory for Solar Technology   phone: +41 56 310 4092
Paul Scherrer Institute             fax: +41 56 310 2688
CH-5232 Villigen                   http://solar.web.psi.ch
Switzerland                   gnupg key: 0x515E30C7






-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.2.3 (GNU/Linux)

iD8DBQFAjRzpAd7CE1FeMMcRApjSAKCk35wEjJSS8XnSXs5y2lqK2Kl4bACePoct
qDx4v1rE54VX5y3EcjqPZUQ=
=2Mhf
-----END PGP SIGNATURE-----
-------------- next part --------------
A non-text attachment was scrubbed...
Name: files.tar.gz
Type: application/x-tgz
Size: 32045 bytes
Desc: not available
URL: <http://sourceware.org/pipermail/gsl-discuss/attachments/20040426/8fd2f48e/attachment.bin>