discrete wavelet transform in reaktor

Very interesting... I'll have to figure out why it seems to work in MATLAB but not here... I know it's NBD since you got it working, but I'm still really curious...

Great job getting it working!

 

i mean, they are all the same values in a different order, which is a problem we've seen with MATLAB lists before.
those guys seem pretty smart so i'd say the error is in my implementation...
except that my implementation was originally based on MATLAB code, so i don't know!

anyhow i'm really quite happy about this, thanks for your help! now i will fix up the display.

 

I think the discrepancy lies with the wonky way (if it is, in fact, wonky) the reconstruction coefficients are configured (that pattern you said you'd never be able to figure out on your own). That didn't come from MATLAB originally - that came from a Java program that I found on the web. It'll have to be sussed out at some point, or at least I need to identify the translation pattern from something like MATLAB, so that people can use their own wavelets, but I'm more excited about doing stuff right now... But I know that you're supposed to just be able to reverse these filters and they'll work for reconstruction.

Speaking of doing stuff, I was curious how the CPU usage scaled (I know it's not optimized yet), so I just made a 5-layer (32 channel) version by simply expanding the busses and copy/pasting the setups that were there (noting that first Modify module). But I don't get any sound. I'm somewhat certain I haven't forgotten anything stupid, but who knows... Anyway, I posted it here to save the effort of having to do it yourself (now you can just copy/paste), and I figured it should have been pretty straightforward, so I was worried about any problems in the implementation, too. Obviously not urgent.

I'm really glad I could help. This has been a really fun project and it's got me back into some things that I thought I would never be able to use again. So thank you, too! And I still have to figure some things out, like these coefficients, so there's that, too.

 

CPU - i believe we've gone from O(1) (xeno's paradox) to O(n), but i'm not sure.

if you aren't familiar with Big O notation - https://en.wikipedia.org/wiki/Big_O_notation

there isn't currently a clock for the 5th layer, so that's probably the issue there. i didn't bother expanding it further because the current scheme is rather inefficient.

there's quite a bit of work to clean this up. we can cut the analyze and reconstruction CPU down by 50%. everything else we're stuck with.

 

fun fact
if every wavelet has the same number of layers, you can delete the delay entirely. i mean, functionally speaking the signal will still be delayed, but we don't need to add any more at that point.

 

No, I put in a clock for a 5th layer...

So I tried duplicating the Dual Tree cell and its MD in both the 16ch and 32ch versions, here are the results on my system (i know there's no formatting):
x..|16ch|32ch
x1| 7% | 7.5%
x2|12%|18%
x3|20%|37%
x4|29%|54%
x5|39%|71%

So in both cases they're eventually linear, with the 16ch version adding 9-10% per copy, and the 32ch version adding 17% per copy. So that lends creedence to your supposition that it's O(n) (if I remember my big O notation properly, that is), if the Reaktor CPU usage indicator is directly (and linearly) related to actual CPU usage (I guess that doesn't really matter since in the system it's the guage that matters).

A single 16ch and a single 32ch together use 13.5% on my system. And one of each of 32ch,16ch,8ch,4ch uses 20% total. More complex setups like peeling off each set of g[n]s and h[n]s from the same tree, just at different levels, and sending them to different Modification Blocks - these I have not tried yet. But I think something like that is what it'll take for the type of multi resolution analysis (MRA) we were talking about earlier (a Modification Block between each level of the Analysis tree).

 

did you also edit the clock receivers to explicity receive them? it's really a bad system right now, i'm just trying to think of a better way that's easy to patch.

 

Ah, found it... Yeah that's a lot of copying and pasting - I'll let you come up with something better before I go ahead with all that...

 

Ok, here are some other fun facts:
For a dual-tree complex discrete wavelet packet transform (DTCWPT), given a set of scaling function coefficients ('h' coeffs):

- For the first tree (i.e. the "real part"), the wavelet function coefficients (the 'g' coefficients) are the scaling function coefficients (the 'h' coefficients) in reverse order with every second one negated (starting with the second one, ending on the last).
- For the second tree (i.e. the "imaginary part"), the scaling function coefficients (h) are the reverse order of the scaling function coeffs in the first tree, and the wavelet function coeffs (g) are the second-tree scaling function ones reversed, with every second term negated.

If my first-tree scaling function (rh0-3 in our parlance) coefficients are: a, b, c, d then my first-tree wavelet coefficients (rg0-3) are: d, -c, b, -a, my second-tree scaling function coefficients (ih0-3) are: d, c, b, a, and my second-tree wavelet coefficients (ig0-3) are: a, -b, c, -d.

I think.

And then there's nothing wonky or unique about how the reconstruction coefficients are arranged - but I haven't found anything that lays out that pattern explicitly (the alternating h/g with evens and odds like they are).

There are some exceptions for dual-tree complex packet:

The First Stage: can use a "simple" scaling function/wavelet, and must be different from the ones used in the rest of the trees. Additionally, the second-tree scaling function coefficients should be shifted one coefficient forward, and the second-tree wavelet coefficients should be shifted one coefficient back so they are two coefficients apart (two filter coefficients equals one sample). That's why the FSfarras filter is so great for this - it has zeroes in front and behind, so you don't need to worry about circular effects (I think). So if my first-tree FS scaling function filter coefficients are: 0, a, b, c, d, 0, then my first-tree wavelet function filter coeffs will be 0, -d, c, -b, a, 0, and my second-tree FS scaling function coefficients will be: d, c, b, a, 0, 0, and my second-tree FS wavelet function coefficients will be: 0, 0, -d, c, -b, a.

The "simple" filters: most filters in both trees come out to be the "simple" filters - the only requirement on these is PR. The first-tree scaling function and wavelet function coefficients have the same relationship as above (reversed, negate every 2nd), but on the second-tree they should be exactly the same, i.e., no reversal of the scaling function coefficients. If my first-tree scaling function coefficients are: a, b, c, d, then my wavelet function coefficients will be d, -c, b, -a and these will be identical for the second tree.

The "non-simple" filters have the same relationships as described at first (before the "exceptions" part).

As far as which are which, well, that's hard to do without a drawing, so let me mspaint it real quick ("q"s are the non-simple filters (e.g. qshift), and "s" are the simple ones, "FS" is first stage):



I've attached a file "filters.txt" that contains all the coefficients used in this project. That includes: "FSfarras" for first stage, "db5" for "s"/simple, and "qshift10" for "q"/non-simple, as well as the relationships for the inverses, so anyone can create their own wavelet-based project starting with the numbers in "filters.txt".

 

Ok, so I made a test application (attached) with 16 knobs, each of which is tied to a Modify macro input which is multiplied by the output of the Delay module tied to the top of xy2polar. The output of the Multiply is the definition of "magnitude".

So basically you can adjust the knob and change the magnitude inside the Modify module.

Some notes:
1) First of all, it works! I can cut out or amplify frequency bands at will!
2) The frequency ranges seem to be powers-of-2 wide, so the first knob cuts out the whole left half, the second knob gets the third quarter, etc.
3) It aliases. Not a whole lot, but yeah, actually quite a bit. Here's a pic:

Except for the speeding hump on the left (attenuated due to knob 16 being turned to 0), the rest of that nonsense is aliasing, or some other form of distortion.

I don't know if maybe filters with larger numbers of points or what is the solution here - that and tighter frequency bands are really the only places left to go. We're supposed to be able to get around this, right?

edit: we could also do the double-density version - the one with 3 bands (or more?) instead of 2. I mean, they invented this stuff for a reason, apparently...

Attached "knobs1" version.

 

heh, i just made a video about doing that with the regular ol' dual tree a few days ago.

i think the aliasing is just a reality. before i worry too much about going down a more technical route i want to:

1) create a system that detects the incoming wavelets and determines automatically which wavelet to send to the g[n] and h[n] for use in the next macro, thus eliminating the need for manual wavelet connections. there seem to be a very simple set of rules to this.

2) same for the clocks.

3) fix the display

4) cut the analyze/reconstruction CPU in half

obviously if you have a sweet new double density wavelet to try out by the time i'm done with all that, that'd be great, but this is what i'll be focusing on for the next little bit.

 

If I understand you correctly, that is what I just posted earlier this evening: the rules for creating the various parts of the wavelet from the original scaling function.

That would be awesome. Please keep me in the loop when you finish working on this or other goals.

I actually did go ahead and flesh out the 32channel version, complete with the proper Clocks in every Analyze and Reconstruct module. I couldn't get it to compile - it would always hang. I left it alone for 10 minutes and still no good. I don't know if your optimizations will address that sort of thing, but I thought it might be a useful note, in any case.

I'm gonna read up on them and see what the bennies are before I go all in.

 

if i had to guess, it's probably the clock setup giving you problems. it creates a highly variable signal path (not really, but in theory). core has a really hard time with that unfortunately.

this one dropped from about 12% to about 8% on my machine, and i haven't even optimized the reconstruction yet (plus there's an FFT analyzer in there as well, eating up around 4%).

i couldn't get the clocks the way i wanted. you'll have to route them; and the wavelets by hand, unfortunately. this is due to some limitations with bundles that i can't figure out how to work around.

 

No big deal, really. BTW, what's the maximum number of fibers a Bundle can carry?

The more I look at the Reconstruction setup, the more I think we're only working on half as many samples as we should be each cycle... I don't think this would necessarily interfere with PR, but it could be why the aliasing is so bad (according to the papers I've looked at there shouldn't be so much, in my interpretation). Don't worry about it, though - the current setup isn't wrong (clearly - it works, after all), just maybe incomplete.

I'll post more on the topic if/when I get something concrete like a working PR that uses 10 points from each sub-band, or at least some better reasoning than "...the more I think...".

 

So I made 0.937 into a 32-channel setup. The good news is that I'm cruising at 3% usage, compared to 2.5% on a 16-channel setup. The bad news is it's no longer PR - there's distortion.


I've attached the ensemble if you're interested.

 

So in MATLAB, using the DTCWPT, on the same chirp I've been using, on level=4 (16 channels), I zeroed out the bottom-most channel and reconstructed the signal and got this:


which is more along the lines of the level of aliasing I think we should expect. I think the difference has something to do with the convolution steps... I'm not sure how the typical convolution process (the one that yields 2N-1 values) compares to what is going on in the Analyze and Reconstruct modules. E.g., If I were to convolve [2 4 6] with [1 2 1], I would expect [2 8 16 16 6] as the result, but it seems that we're only doing the middle-most step, the one in which the filter and the samples line up exactly.

Plus, I still have only a very limited idea of what is going on in the Reconstruction modules wrt multiplying the low-pass inverse constants by both the High and Low-pass values (same for hi-pass inverse constants)... Why does *that* work? Everything else I've read online shows the inverse lo-pass filter being applied to the lo-pass values, and inverse hi-pass filter being applied to the hi-pass values. My only guess is that, since the algorithm calls for adding up the results of the lo-pass multiplication and the hi-pass multiplication anyway, that we're just doing that all in one step...

Is it possible that there's a more involved convolution step?

 

do you have 32 modfiy macros in that setup? they need to be there for the delay. i'll check it out on my work computer soon.

as for convolution, i really don't think that's the problem, but i could be wrong. right now we are doing linear convolution which is simply an FIR filter (IE the analyze macros). the more complex setups for circular convolution only work for periodic signals. plus they essentially fill the filters with 'garbage' so i don't think they will reduce aliasing, but add to it.

 

okay, here is my contention:
perfect reconstruction (of an unmodified signal) works.
if it ever doesn't the fault is solely ours, at least with these wavelets. the only exception can be numerical error, which we've already proven with the dual tree, should be infallible at least up to 12 layers of deconstruction.

you can of course prove PR to yourself simply by doing one-layer transforms, then recognizing that for any other setup
a) the last set of filters will perfectly reconstruct the last sets of inputs.
b) the next set of filters will perfectly reconstruct the perfectly reconstructed values it receives.
c) and so on and so forth.

so your problem was that the 'max' input to one of the modifies is disconnected. connect it and everything will work.

that's one reason why i want each signal to contain it's own layer info, so the modify macros can calculate delay without external assistance, reducing the number of errors like these, which are quite simple to make and overlook. another reason is that each signal may have a unique number so the global constant doesn't work anymore.

 

Ok, lesson learned... My misgivings come from not understanding exactly why this setup is working when it is so different from almost every other piece of code on the web (on the Reconstruction side). But that's my problem to worry about and I'll keep it out of the discussion until I have something to actually say about it.

I was really surprised by the "failure" at 5 levels because of this exactly. I thought there might be some small errors propagating through that were getting magnified for some reason - crazy stuff, because I thought I had double-checked all the connections.

D'oh! Jeez, thanks - man, I'm really sorry to have wasted your time on that... I can't believe I missed that over and over. I swear I remember hooking that guy up - I must have missed the spot when dropping the connection on there...

 

presumably we could rearrange the 'first' wavelet to look the way you think it should, the re-arranage all the values in reconstruction accordingly (again, the values are the same, in a different order). maybe if you did that you'd get something that looks more sensible?

i have been meaning to do this, but there are so many improvements to work on, constantly.

don't worry about it, i'm glad the 5 layers works after all, and it only took a few minutes. i've gotten very good at scanning large reaktor documents for small errors, people send me questions like that all the time.