N = 6 scaling functions and mother wavelets

introduction

The impetus for this post is figure 1.4 on page 6 of Burrus et al. (Since it’s been a while, that’s Burrus, C. Sidney; Gopinath, Ramesh A.; Guo, Haitao. Introduction to Wavelets and Wavelet Transforms, A Primer. Prentice Hall, 1998.ISBN 0 13 489600 9.)

They offered the figure as an illustration of four different scaling functions; in addition, these scaling functions were parameterized by two angles.

Now I know how to produce the drawings of those scaling functions — and I also know how to produce drawings of the corresponding mother wavelets, which they did not show on page 6. It also turns out that they have interchanged the legends for figures (a) and (b) — and when we can produce the drawings ourselves, that mistake becomes almost irrelevant.

This post also serves to reiterate the calculation sequence for plotting a scaling function and its corresponding mother wavelet.

setup

It turns out that there are equations which must be satisfied by any “reasonable” filter coefficients h. I expect to show these to you in the next post. I have already shown you a description of these solutions of these equations when we have 4 nonzero coefficients. And I have shown you how to compute scaling functions and mother wavelets for N = 4 (D4 only) and for N = 2 (D2 = Haar).

Let me emphasize that what I understand is how to get the equations which the filter coefficients h must satisfy. What I do not yet understand is how to parameterize the solutions in terms of angles.

Let me now hand you a description of the solutions when we have 6 nonzero coefficients. The four examples I am about to work out all come from specific choices of these solutions. That is, for every one of the following examples, N = 6.

The following 4 scaling functions match figure 1.4, although (a) and (b) are swapped. (The legend describes the values “a” and “b” in the equations; the legend for figure a belongs with figure b.)

In addition, they have drawings of the Daubechies D6 (p. 82) and Coiflet C6 (p. 94) mother wavelets, so I have confirmed 6 of the following 8 graphs.

From p. 66 of Burrus et al., we get the following (using H instead of h just because it lets me preserve the equations in H while setting values for h.)

Let me tell you before I start writing that the last two equations are actually

That is, they say that the sum of the even coefficients is equal to the sum of the odd coefficients is equal to half the sum of all the coefficients = . They will look much more complicated than that when I substitute the first four equations.

legend 1.14 a (figure b)

The first choice (the legend for figure 1.14 a) for the angles a and b (in radians) is

(I knew when I saw those that the figure was wrong; by that time I recognized the D6 h’s.)

The resulting h’s are…

Now we form the m0 matrix (in stages for clarity)…

Then we multiply by Sqrt[2]…

Note that since all 4 examples have N = 6, that is the m0 matrix in principle for each of these examples. I won’t show it again.

Now we set the values…

We find the eigenvector of eigenvalue 1. Mathematica returned a unit vector, but I want one whose components add up to 1. I will simply divide the unit vector by the sum of its components.

Here’s the unit vector…

The sum of its components is -0.742756, so I get

This function should be defined on [0,5]. I use that vector for the values of the scaling function at the integers.

Working, I graphed this immediately, but for the post, I ought to put the two drawings — scaling function and mother wavelet — in close proximity. This one turns out to be the Daubechies D6 scaling function.

To get the D6 mother wavelet, we take our filter coefficients h…

… we reverse them, and alternate the signs. The resulting filter coefficients are often denoted by h1 or by g. I’ll use g. The resulting g coefficients are what I expect (the sign of g[0] = h[5] has not been changed), but their Table 6.2 on p. 79 has the negatives of these numbers. And yet, I match their picture on p. 82.

(That’s why there’s a in front of the recipe. Sorry, but you’ll have to look for it in the post.)

Define the mother wavelet…

and plot it…
(These are Daubechies D6 mother wavelet and scaling function, respectively.)

Beautiful. That matches their picture on p. 82.

And here’s the D6 scaling function whose plot I delayed showing you…

and that matches (b), rather than (a), of figure 1.4.

legend 1.14 b (figure a)

Let’s make a different choice for a and b. I believe this is the “Coiflet” of order 6.

The resulting h’s are:

Now we form the m0 matrix exactly as before, and set the new values…

We find the eigenvector of eigenvalue 1.

That’s orthonormal, so I find the sum of its components (0.578913) and divide by it, getting

Set those to be the values of the scaling function at the integers, and define the recursion…

Again, I will delay the plot until I have the mother wavelet.

To get the mother wavelet, we start with the h’s…

(Those agree with p 92)

We reverse them, and alternate the signs. (This time, I need to use their signs! And they have got g[0] = -h[5].)

compute and plot…
(These are the Coiflet C6 mother wavelet and scaling function in order.)

That matches their p. 94 drawing of the C6 mother wavelet.

And here’s the C6 scaling function…

and that matches their figure 1.14 (a).

figure 14.1 c

I do not know if the following pair has a name. We take

The resulting h’s are…

Now we form the m0 matrix, and set the values.

Now we find the eigenvector of eigenvalue 1:

The sum of its components is -1.29645, so we divide by it, to get a vector whose components add up to 1:

That gives us the initial values, and we define the recursion…

Recall the h’s…

We reverse them, and alternate the signs.

Redefine the mother wavelet…

and plot it…

… followed by the scaling function:

The last case for their Figure 1.4.

Here are the angles a and b…

Here are the resulting h’s…

We form the M0 matrix and set its values…

Once again we find the eigenvector of eigenvalue 1. (If you’re using Mathematica, be careful. It’s probably the 2nd one rather than the 1st, in this case.) The normalized eigenvector is…

The sum of its components is 1.39807, so we divide through and get our initial values for the scaling function…

So we define the scaling function

Here are the h’s…

We reverse them, and alternate the signs to get a set of g’s…

Define the mother wavelet…

Compute and plot it…

And here’s the scaling function (which we defined first, as usual)…

Closing

So. We’ve had some practice computing the scaling function and the mother wavelet from the h’s. We know that the h’s are not arbitrary, but are given by solutions to some set of equations, and I’ve shown you the general solutions for N = 4 and N = 6.

Now I really owe you an explanation of where the heck these computations come from. But at least you’ll know why I’m writing equations when you see them.

Next, I think.

Mother wavelet from scaling function: D4 and Haar

D4 dyadic wavelet

We can compute the mother wavelet from the scaling function. Let me show you how to do this for Daubechies’s D4.

First, we need the scaling function.

the D4 scaling function again, quickly

This is review. There is a matrix M0 which has an eigenvector with eigenvalue 1, and that eigenvector gives me the values of the scaling function at the integers. Once I have initial values, I can compute by recursion (and because of the specific form of the recursion, only at points whose denominator is a power of 2).

I start by constructing the matrix M0 in principle…

(This comes from evaluating the dilation equation at the integers, then writing the resulting system of equations as a matrix-vector equation. It turns out to be . In fact, I construct M0 from its pattern, but i know how to derive – hence check – it.)

Then I set the values of the filter coefficients h…

The matrix M0 becomes…

Now I look for an eigenvector having eigenvalue 1. Here are the eigenvalues (first row) and eigenvectors (subsequent rows)…

We see, as we saw before, that we want the first eigenvector from that list. Mathematica has given me a vector of length 1; I want a vector – I’ll explain later – whose components sum to 1.

I take those components as the initial values of the scaling function, i.e. its values at the integers. Then I write a recursive function definition, which includes the requirement that the function is zero outside the interval [0,3].

Frankly, the easiest way to show it to you is a drawing, especially of what Mathematica knows after I make the definition; read the last line.

In order to make a point, I am not going to evaluate that function anywhere else. I have all I need for now.

the D4 mother wavelet

Let me hand you another recipe. The “mother wavelet” satisfies an equation very much like the dilation equation. We have

We do not have to specify that the function is zero outside [0,3]; all we have to do is feed it .

It’s important that it is not a recursion for : the LHS is but the RHS is still the scaling function . And the coefficients have changed. They are often denoted by h1 or by g. I think I’ll use g.

The g(n) are related to the h(n), but not uniquely. In general,

, for L an even integer (I’m pretty sure it has to be even).

An extremely convenient choice is L = N, the number of nonzero h’s, in which case we get

,

which simply says reverse the h’s and alternate the signs.

That is, for our case, N = 4, so we’re looking at h(3-n), and then

g(0) = h(3)
g(1) = -h(2)
g(2) = h(1)
g(3) = -h(0)

(or the negative of all those).

Here are our h’s again…

or

We reverse them, and alternate the signs, in order to get the g’s.

or

Oddly enough, Burrus et al. showed the opposite of those signs. That doesn’t agree with their own equation, and it would lead to the negative of the following graph – but my picture matches their picture!

Now we define in terms of …

It’s time to just plot the D4 mother wavelet.

and that looks good.

NOTE that it will compute whatever it needs of , so in a sense it is recursive. This is the point I wanted to make: I did not have to compute for it, it could do it itself.

Just to put them near each other, here’s a picture of the scaling function:

the Haar wavelet system

Let me show you all of that again for a much simpler case. The Haar wavelets, as they are now called, was invented long before we knew they was a wavelet system. Its scaling function is a step function, defined to be 1 on the half open interval [1,0) and zero elsewhere. I will eventually show you – I hope – that it could be considered the D2 scaling function also.

Oh, here I go again using a cannon when a slingshot would more than suffice. There is no reason to use recursion to graph this scaling function, nor even to graph this mother wavelet. But it’s nice to know that the algorithm works in the simplest case.

Haar (D2) scaling function

We start with the M0 matrix in principle, for two nonzero h’s…

I simply tell you that the h’s have these (two identical) values…

Then the matrix M0 is…

Well, well. By inspection – it’s diagonal! – there are two eigenvalues equal to 1, and the eigenvectors are the columns of that matrix. We might as well choose the eigenvector

The sum of the components is already 1, so we don’t need to scale it to get the initial values of the scaling function.

We define the scaling function as usual…

Let’s plot it.

A unit step function on the half-open unit interval. I’m rather pleased to see that.

D2 wavelet

To get the g’s for the mother wavelet, we reverse the h’s (oh, they’re the same!) and alternate the signs…

We write the definition… ask what we know… and then we plot the mother wavelet…

Yes! That’s exactly what I wanted to see. And if you’ve ever seen the Haar (D2) mother wavelet, it’s what you wanted to see, too.

I think the next post will show you a few more scaling functions and mother wavelets, for six nonzero h’s.

The dyadic expansion of Daubechies D4 scaling function

introduction

edit jun 13: corrected one of the 4 equations at the end.

I have shown you one way to compute the D4 scaling function. I am not entirely comfortable with the method, but I don’t need to be: I can do better.

The method I’ve shown you begins with a terrible approximation — a constant function or signal — and improves it by repeatedly applying a fixed filter to it.

The primer by Burrus et al. also provides MATLAB® code for a quite different method. It turns out that we can find the exact values of the D4 scaling function at the integers. The dilation equation can then be used to find the exact values at the halves… then the exact values at the quarters… and so on. It’s called a dyadic expansion because it works for points whose denominator is a power of 2.

Their code, however, was challenging to figure out.

But I don’t need to figure it out. I can do better.

I can do so much better that it must be illegal, immoral, and fattening.

I don’t have to write any code at all — well, ok, I have to write an “If” clause. Mathematica® can do recursion, and the dilation equation is all I need. Oh, I need initial values, too.

Rather than distract you from this marvel by putting the details first, I’m going to show you the end product first.

the dyadic expansion simply & efficiently

To do a recursion, I need to have initial values. I’ll derive them later in this post, but the values of the scaling function at the integers 0, 1, 2, 3 are components of an eigenvector V:

I set those to be the function values at 4 points…

and then I define the function. (Whoa! I have to use a picture, apparently, because of special characters.)

(Actually, I can make that work. The problem is the < and > symbols, which are apparently being interpreted as containing nonsense HTML. Use & lt; and & gt; without the space after the “&”, which is what I did to make the code display.

The set delayed operator (:=) says don’t do anything until I give it an actual input. The in the section just says “after you compute it, assign it”.

There is a way to see what Mathematica® knows at this point.

Let’s compute it at 1/2. We get

and let’s see what Mathematica® knows now:

He knows a lot, more than when he started. And he will return from the table as soon as he finds any known value, and fall thru to the last line — the recursion — only for unknown values. This is a combined lookup table and recursion.

(I know that both Java and C++ are recursive; it appears that MATLAB® is not.)

Now that is nice. It shows that Mathematica® needed a wide range of values for the computation, and it saved them — and, more importantly, the next time it needs, for example, , it will find it in that list before it gets to the final entry, which is the recursion.

I must emphasize that we are computing exactly, at new points. The values at old points are never changed in this process.

OK, let’s see what we get. There’s absolutely no computational reason to proceed in stages, but I want to make sure you see it grow.

Here I show the initial values ( at the integers)…

We add the values at the halves.

We add the values at the quarters. And I can disdainfully tell it to “recompute” the values at the halves because I know it will look them up rather than actually compute them.

Add the eighths…

Sixteenths…

Thirty-seconds…

Sixty-fourths…

One-twenty-eighths…

I note, from the position of the dots to the right of the maximum, that the function is changing very rapidly there. I wouldn’t be surprised if the convergence in our previous method is not very quick there.

getting the initial values

Let me now back up and show you how I got the scaling function at the integers. Oh, understand that it is zero outside the interval [0,3].

Here is the dilation equation for N=4, that is, for four nonzero filter coefficients h.

First, evaluate the dilation equation at the integers 0 .. 4. If we simply set t = 0, we get

but then we need to zero out the three values of at the negative integers in that equation, getting

For t = 1,

and this time we set , getting

For t = 2…

and we set , getting…

finally, t = 3…

and we set three of the to zero, getting

so, we have the system of equations

and we could solve them.

It is more instructive and useful to look at a matrix representation. That set of equations can be written

with the matrix M0 defined — in stages — as

and then I multiply it by , getting

(I created M0 using the diagonal bands, long before I knew how it was derived. Things were a little clearer to me without that ubiquitous .

Oh, and the vector holds the 4 values of at 0, 1, 2, 3:

But that matrix equation says that is an eigenvector of M0, with eigenvalue 1. And we can find eigenthings.

To solve it, we need the values of h for the D4 scaling function.

For the D4 scaling function, then, M0 becomes

or

Here are its eigenvalues (first row), and corresponding eigenvectors…

We see that the first eigenvalue is 1, so the first eigenvector is the one we want.

That is a vector of length 1. For reasons we will discuss later, I want a vector whose components add up to 1. It’s still an eigenvector.

(The signs changed because I simply divided by the sum, and the sum was negative.) We take those as the exact values of at the integers 0 thru 3. You see that ,
as I have said before.

So: the initial values of the scaling function are the components of an eigenvector of a matrix M0, specifically of an eigenvector with eigenvalue 1. That is how I got them.

the simple recursion

Of course we’ve already set up a recursion combined with a lookup table. I just want to emphasize the distinction between what I actually do, and the unadorned recursion.

We are ready to write the recursion. It looks almost exactly like the dilation equation…

but instead of “==” I use the “set delayed operator “:=”, and the LHS “t” must become a placeholder “t_” ; and the whole definition must be wrapped in an If-test to set all values to zero outside of [0,3] .

Oh, since I’m starting all over again (don’t want to get two sets of definitions mixed up), I need to set the values of h…

Again, I define the function values exactly as before, at 4 points…

And then I define the recursion:

Again, I can code that now, instead of using a picture.

The set delayed operator (:=) says don’t do anything until I give it an actual input. What’s different in this definition is instead of . There is no assignment statement to saying that the result of the computation should be saved.

Let’s compute again. We get the same answer, of course:

and let’s see what Mathematica knows:

That’s exactly what we knew before we computed it at 1/2. This is a pure recursion. He recomputes everything every time he needs it.

Simple, but potentially expensive in computer time.

Summary

From the dilation equation and the h coefficients, we can get the values of the scaling function at the integers. we do this either by directly solving a system of equations, or by recognizing this particular system of equations as an eigenvector problem

.

In either case, if we are at all uncertain what the matrix M0 looks like, we can construct it from the equations.

Note that because the values of at the integers are the components of an eigenvector, they are only determined to with a scalar multiple. That’s ok, we have a reason for normalizing so that the sum of the values is 1. (That, of course, is different from the sum of the squares of the values = 1.) And even that normalization will only get is to within a sign. If you ever plot a scaling function and get the negative of someone else’s picture, all you need to do is change the signs of the inital values at the integers, i.e. change the sign of the eigenvector.

Having the values of at the integers, we can get it at the halves, then the quarters, etc.

As I have said before, the purpose of this algorithm is to see what the scaling function looks like, and it will work for the mother wavelet, too.. In practice, however, we can compute the wavelet coefficients of a signal uisng the filter coefficients h (and the corresponding coefficients for the mother wavelet) — we don’t need the values of the scaling function or the wavelet to get the wavelet coefficients.

(I believe that the previous algorithm, however, serves another purpose. It seems more than plausible that the previous algorithm is used to prove existence and uniqueness of the scaling function, once we have its filter coefficients h.)

For comparison with the next case, let me note that the 4 h’s satisfy the following 4 equations…

edit jun 13: the third and fourth equations were identical. i have corrected the third equation by changing the sign of the cosine term.

wrong:

right:

The Daubechies h’s are found by setting .

All I’ve done is push back the mystery of where the h’s came from. Don’t worry. We’ll get there.

I had thought that I would show you plots of other scaling functions, but it makes more sense to show you how to get the mother wavelet.

Next. Unless something really important comes up.

Wavelets: The D4 scaling function

Introduction

I am going to do something a little risky. I’m going to show you how to approximate a particular function – the Daubechies D4 scaling function – even though I am a little shaky on the justification.

If you have ever looked at wavelets, you have probably seen this picture. It is a quintessential “father wavelet” – more commonly, I think, called a “scaling function”.

I want to compute it; this is not relevant to computing wavelet coefficients, because just as this is computed from its definition, so wavelet coefficients which depend on it and its associated wavelets are computed from its definition rather than from its values.

Talking about this is much like talking about one of one’s favorite movies: what makes something one of my favorite movies is deeply personal, and it doesn’t usually affect my friends that way. To them, at best, it’s just another movie; at worst, they think it is a dud.

I have wanted to know how to compute this function for a long, long time. You have noticed by now, however, that I tend to bounce from one subject to another; my desire to compute this function was put on the back burner.

Let me say again that in practice we do not compute this function: we use its defining characteristics instead of its actual values. But I just had to know how to do it – someday.

Well, it has been “someday” for a little while now; and I’ve been struggling to assemble this post.

It has gotten me to pick up wavelets again, and to look through my books on the subject. This particular computation has gotten me to look at some more challenging Matlab code, and it will lead me to study convolution. (Matlab code would be far less challenging if I had access to Matlab!)

This function is defined on the interval [0, 3]; it has no analytic formula. It is the solution of a difference equation, called variously “the dilation equation”, “the refinement equation”, and the multi-resolution analysis (MRA) equation”:

and in this case, the scaling coefficients h are:

It is called D4 because it was invented by Ingred Daubechies and it has 4 nonzero h’s. As I said, it is called a “scaling function” in contrast to “a wavelet”. We’ll see why, down the road.

If we look at wavelets, we see this function very often. But most books do not explicitly show us how to calculate it; and the calculation is not completely straightforward.

The method I am going to use comes out of DSP, discrete signal processing. It is not a tool I use all that comfortably, but the calculation itself is rather simple. On the other hand, details keep piling up, so I am going to split this post.

The method is iterative. We will start with a constant signal, and apply a filter to it eight times. I printed a picture of the transformed signal after each iteration, and I showed it to two friends who had some familiarity with DSP – but who were not practitioners of it.

They both found the algorithm bewildering. I now understand it in general, but there are still details I wonder about. Nevertheless, what I did to my friends was just wave the pictures in their faces, and mutter about upsampling, downsampling, and convolution.

In this post I want to show you the pictures, the Mathematica code, and the underlying idea. In the next post I will fill in some of the underlying details.

Overview

My starting point was Matlab code in the back of a book. If you have Matlab available to you, duplicating my calculations will be very easy – because their code is also available on the Internet.

The very beginning of the appendix containing the code says

“You are free to use these programs or any derivative of them for any scientific purpose but please reference this book. Updated versions of these programs and others can be found on our webpage at: [http://dsp.rice.edu/software/introduction-wavelets-and-wavelet-transforms-primer]“

except I have corrected the URL (replacing dsp-rice by dsp.rice).

The book is: Burrus, C. Sidney; Gopinath, Ramesh A.; Guo, Haitao. Introduction to Wavelets and Wavelet Transforms, A Primer. Prentice Hall, 1998.ISBN 0 13 489600 9.

If you want to see or to use that code, please go get it; while I have used – and will show you – derivative code of my own, I will not actually publish their code. (After all, they said I could use the programs; they didn’t say I could publish them.) Oh, the specific function you want is their “psa(h, kk)”.

They define two completely straightforward functions: downsample and upsample. Here is my Mathematica code for them.

dnsample[x_List] := Table[x[[i]], {i, 1, Length[x], 2}]

upsample[x_List] := Table[{x[[i]], 0}, {i, Length[x]}] // Flatten

The only subtlety in the downsample code is that it drops the even-indexed terms (rather than the odd ones), counting the first term as index 1. Similarly, the upsample code inserts zeroes after the 1st term.

Then they define an upsample function which takes an additional (integer) argument, S, and instead of inserting a single zero between original values, inserts S-1 zeroes. Yes, I do know why their argument was S instead of S-1. You’ll see why, too. Here is my mathematica code for it.

upsample[x_List, S_Integer] :=
Table[{x[[i]], ConstantArray[0, {S - 1}]}, {i, Length[x]}] // Flatten

(Because Mathematica understands parameter lists, I do not have to change the name of the second upsample function; it distinquishes just fine between the function upsample with two arguments and the function upsample with one argument.)

The tricky point for me in the code is the Matlab function “conv”. As far as I can tell, the following Mathematica code is the equivalent.

conv[x_List, y_List] := ListConvolve[x, y, {1, -1}, 0]

(The challenge I face is understanding the options for the ListConvolve command. I would not bet the farm – I might bet a Dover paperback or two – that my Mathematica command is always equivalent to the Matlab command.)

The algorithm

Now let me show you the code that generates the scaling function.

p = dnsample[conv[hu, p]] // N;
Take[%, 384];

The first line computes the convolution of two sequences, hu and p; it downsamples the result; and then it overwrites the sequence p. (The //N converts from rationals to reals.) The second line says that we are only going to plot the first 384 values in p. To rephrase that, this approximation to the scaling function is the first 384 points of the signal p.

(There are other ways to approximate the D4 scaling function; this is the most straightforward – even though I am not all that comfortable with convolution – because it repeatedly applies the same filter, “convolution followed by downsampling”, to an initial constant signal. I had originally thought to look at a different method in this book, until I realized that all I had to understand in this case was one filter.)

What exactly is hu? And where did the 384 come from?

We start with the h coefficients which define the D4 scaling function:

We multily by , and give the result a new name, h2:

Why? Recall the dilation equation:

All we have done is absorb the into the h’s. Our dilation equation is now

Our function is defined by 4 values of h, and – because of these 4 specific values, I think – our function is nonzero on the interval [0, 3]. We are going to compute it at 128 samples per unit interval, and we have 3 intervals, so we take 384 = 3 * 128 points.

In Mathematica, I used “rules” to set two parameters, K = 3 = 4 – 1, and S = 128, as follows:

ps={K->Length[h2]-1,S->128}

which gives me

For this post, however, I have chosen to work with these specific numbers; my original Mathematica notebook – in contrast to the notebook for this post – is still “rule based” and I can change the parameters very easily.

Now comes a piece of real magic: we upsample the filter coefficients h2, introducing 127 zeros (and so they wrote the code to do the subtraction for them, 127 = 128 – 1). This we call hu (u for upsampled) and the function call is:

hu=upsample[h2,128];

As a result, the four nonzero values of hu are at indices 1, 129, etc.:

There is one of my first questions: since the final filter coefficient is at index 385, why did they not plot the first 385 values of p?

In any case, hu is an upsampled set of filter coefficients. My friends and I were quite bewildered. We had all heard of upsampling a set of data – but upsampling the filter coefficients instead?

This, it turns out, is quite simple. It also has nothing to do with the up and down sampling characteristic of the filter bank approach to wavelets.

By taking 384 or 385 points, we have changed our timescale from [0, 3] to 1 .. 384 or 1 .. 385. That means we have to move the filter coefficients from the integers 0, 1, 2, 3 to the integers 1, 129, 247, 385 (or possibly to 0, 128, 246, 384, but I choose to start indexing at 1 instead off at 0.)

So much for upsampling the filter coefficients h.

What about the iteration itself? Instead of the dilation equation, the authors write (note the superscripts (k) and (k+1) ):

This is reminiscent of a fixed point problem

x = f(x)

which we may try to solve as an iteration: guess an x; compute f(x); now call that x because the equations says f(x) is x; if it converges, we have an answer. That is, the iteration:

might converge. On the other hand, the process might not converge, and even if it does, the answer need not be unique.

We have a slightly different situation because the RHS is at time 2t-n while the LHS is at time t, and it’s the function we’re trying to find, rather than x (or t). Still, if we can find a nice starting function , we can compute , and then , etc.

The idea looks very plausible to me; what threw me for a loop (no pun intended) was digitizing the dilation equation. Understanding that upsampling the filter coefficients was just a change of scale, went a long way toward understanding the algorithm.

I need to show you two things. One is the filter which is applied repeatedly, to get us from to . The other is the initialization, .

Here is the initialization step:

The first line,

p = ConstantArray[1, {1152}]/9;

says we start with 1152 values of p, each of them equal to 1/9. Let me be clear: I do not know why we choose a power of 2 (2^7= 128) points per unit interval – rather than, say, a power of 2 points for the entire interval [0, 3], nor why we initialzied to 1/9; nor why, with 384 points of interest, we initialized our constant signal to 3 times that length (1152 = 3*384 = 9*128).

The second line,

Take[%, 384];

says that we are only interested in the first 384 points, i.e. we are only interested in the answer over the original range [0,3].

And, as I said before, I don’t know why we look at 384 points instead of 385.

The third line,

data = Transpose[{Range[384]/128, %}];

is used only to change the horizontal scale for plotting purposes: I don’t want a horizontal scale of 1 .. 384, but rather [0,3].

The fourth line is the specific graphics commands and options which I used. Unless you have David Park’s “presentation” package, you will need to use the Mathematica “ListPlot” command instead of “ListDraw”.

dwg[0] = Draw2D[{
Thickness[.01], line, ListDraw[%, Joined -> True]
}, Axes -> True, PlotRange -> {All, {-.5, 1.5}}, Background -> back,
AspectRatio -> 1/GoldenRatio]

Although a plot of the initial function seems superfluous, let me err on the side of too many pictures.

I remind you that I have plotted the first 384 values, but the signal p starts out three times that long.

Okay, let’s apply our filter. Here is a snapshot of my code and the plot:

The first line of our code

p = dnsample[conv[hu, p]] // N;

says convolve hu with p, downsample, and replace the signal by that output. The convolution corresponds to the RHS , and the downsampling should have the effect of setting the result to time t rather than time 2t: rather than

Once again, the second line of the code,

Take[%, 384];

says take the first 384 values of the sequence p. Note that p itself was redefined in the first line; the second line selects a subset of p, but the next iteration will use more than the 384 points selected on the second line.

The third line, as before,

data = Transpose[{Range[384]/128, %}];

just changes the scale on the horizontal axis. The fourth line, as before,

dwg[1] = Draw2D[{
Thickness[.01], line, ListDraw[%, Joined -> True]
}, Axes -> True, PlotRange -> {All, {-.5, 1.5}}, Background -> back,
AspectRatio -> 1/GoldenRatio]

plots the first 384 values of the current sequence p.

The plots and summary

The initialization, a constant signal = 1/9, :

After one application (convolve and downsample), :

After the second application, :

After the third…

After the fourth…

After the fifth…

After the sixth…

After the seventh…

Finally, after the eighth (which is the picture I started this post with)…

We could worry about how close the answer is, but – as I said – we don’t need this function numerically; what we need in practice is its 4 filter coefficients h.

The summary is simple enough. I am delighted that there is so little overhead involved in this algorithm: convolve and downsample, repeatedly. This is slick and quick.

I do have more details to show you. Soon.