[Numpy-discussion] How to code functions like a Haar wavelet?

[Paul Barrett]

Roger Flores wrote:

> I've read through the docs but all the filter examples seem to affect 
> one value at a time.  How do I write something like a wavelet that 
> transforms multiple values at the the same time?
>
> I'm hoping for a Python version (as opposed to a C function) for 
> portability and easy exploration.

These algorithms can be done in Python as you might do them in C except 
that the looping is implied.  Assuming 'y' is an array of length N.

The Haar transform can be written as:

s = (y[1::2] + y[0::2])/2.
d = (y[1::2] - y[0::2])

where s is the average (smoothed) values and d is the difference 
(detailed) values.

In the lifting scheme (i.e. using filter banks), it is:

y[1::2] = y[1::2] - y[0::2]
y[0::2] = y[0::2] + y[1::2]/2.

Note that y[1::2] is the difference and y[0::2] is the average. And that 
the order is important, since the operation is being done in-place.

For higher order wavelets, boundary effects must be taken into account. 
For example using a cubic B-spline wavelet, the algorithm is:

s = (y[0:-4] + 4*y[1:-3] + 6*y[2:-2] + 4*y[3:-1] + y[4:])/16.
d = y[2:-2] - s

Note that this expression does not calculate the first and last two 
values.  These must be done as special boundary cases, e.g.

s[0] = (6*y[2:-2] + 4*y[3:-1] + y[4:])/11.
s[1] = (4*y[1:-3] + 6*y[2:-2] + 4*y[3:-1] + y[4:])/15.
s[-2] = (y[0:-4] + 4*y[1:-3] + 6*y[2:-2] + 4*y[3:-1])/15.
s[-1] = (y[0:-4] + 4*y[1:-3] + 6*y[2:-2])/11.

d[0] = y[0] - s[:3]
d[1] = y[1] - s[:4]
d[-2] = y[-2] - s[-4:]
d[-1] = y[-1] - s[-3:]

or by padding the array.

I hope this helps.

 -- Paul

-- 
Paul Barrett, PhD      Space Telescope Science Institute
Phone: 410-338-4475    ESS/Science Software Branch
FAX:   410-338-4767    Baltimore, MD 21218