Could use some help with this interpolator for resampling...

sambean · 17 May 2006

If you are looking to resample by a 2^n ratio I think you're best sticking to FIR or IIR polyphase.. I had a read of the first half of the introduction for the paper and it seem to me to be for resampling by an arbitrary amount. The process it suggests is to first oversample using FIR/IIR polyphase and then resample the oversampled wave to get your desired resampling ratio.

You could use this method to 2x oversample by choosing x = 0.5 and calculating the return value from the y[n] adjacent samples.

Lalo made a start at implementing some of the algos from that paper: interpolations : linear, cubic and spline

But, I don't think you could use any of them for downsampling, so you'd still need to revert to a filter solution for this..

Here's a link to an IIR polyphase half band filter in C/C++ which probably can be ported over and used for both resampling stages.. http://www.musicdsp.org/archive.php?classid=3#39

As to the pre-emphasis, I havn't a scooby 😀

Cheers

Sam

sequencer · 17 May 2006

Don't know much about new fangled DSP stuff but, in the olden days Pre-empahsis was used (and still is) for Analogue Audio (especially in Records, remember them...?).

The process basically added 6dB to the Spectral Frequency above a certain frequency when the record master was cut. At the consumer end when playing the record, De-emphasis was applied at the same Frequency point to obtain the original signal back. This process helped with lowering the noise floor so, perhaps the same is applying to the interpolation in some form and reducing the noise floor...er maybe, don't really know on that though 😉

Dunc (showing my age)

Kingston · 17 May 2006

sambean wrote: Here's a link to an IIR polyphase half band filter in C/C++ which probably can be ported over and used for both resampling stages.. http://www.musicdsp.org/archive.php?classid=3#39

Thanks sambean.

*sigh*

I've been reading this for more than an hour now, and I simply cannot get my head around how to implement this in synthmaker. I'll try again tomorrow, but it goes over my head.

which is weird since I did understand your 4x sinc interpolation you posted around here.

I used it with good results, but I suspect this one could be ported to better realtime use as 2X upsampler/downsampler. if I only understood the process..

ie. which of the functions does upsampling, which one downsampling, and how to obtain the two signals (the other one obviously with 1 sample delay) and put them back together.

Kingston · 17 May 2006

I did some reading,

http://users.cscs.wmin.ac.uk/~krukowa/Poly.html

and this is how far I got with synthmaker

[removed buggy file]

I'm not quite sure what to do with it, but I think I got the allpass filters right. I used the 8th order 0.05 transition band version. What is that half-band filter?

I just don't get it. 🙁

sambean · 18 May 2006

It seems to be working.. Subtract one output from the other and pop the result in an analyser, if you look at the FFT you'll see it's a pretty clean cut at half nyquist..

Getting this working in SM may be a little difficult though as there's no way to run at a faster samplerate. you'll need to represent your oversampled data as 2 values, loop the filter code (or just run the same code twice on the 2 values), do whatever you need to do with the oversampled data, then reverse the process..

The filters working, so I'd attempt to do the above with a filter of lower order first, then once it's working, replace it with the higher order one..

If the latency is the problem with the FIR version, you could also look at replacing the linear phase FIR filter coefficients with a minimum phase version..

Kingston · 18 May 2006

As far as I understood it, in synthmaker I could simply process the two streams in parallel and get "double precision" in comparison to non-oversampled material.

I only need to take care of any frequency dependant coefficients in filter etc. processing.

My idea was that I could run 2x oversampled stereo audio (ie. 4 separate streams) "for the price of one" in synthmaker, thanks to SSE.

the only extra processing would come from the actual up/downsampler and the pack/unpack stream functions.

Then downsample back to single stream. I still don't get it, but I'll try a hands-on approach, if I miraculously get it working somehow.

But what is that halfband filter?

rl · 18 May 2006

Kingston wrote: But what is that halfband filter?

The transition region of a half band filter is centered at one quarter of the sampling rate, or fs/4 (or fnyquist/2) to be used in decimators. IIRC its mostly implemented as a FIR filter and has the advantage that half of the coefficients are zero for faster computation.

Kingston · 18 May 2006

I tried some testing,

the filter seems to be working since the steep lowpass is near perfect when summed. This thing is also EXTREMELY efficient with stereo streams.

I can process the whole thing in one pass using packing and unpacking with SSE streams! (2xchannels, 2x parallel oversampled processing=1 stream)

BUT.

It's not working, have a look at this

[removed buggy file]

I don't seem to understand how to reverse the process. Or maybe I cascaded the allpass filters wrong? What about the one sample delays, are they right at all?

[EDIT]

holy shit I think I got it! It was a problem with the reverse cascade. This thing will be extremely useful for us all if it does indeed work as well as I think it does.

I'll do some more testing, but damn this thing is EFFICIENT with parallel processing!

I'll post the final results later today after I'm done making sure it's 100% correct.

Kingston · 18 May 2006

ok here it is,

http://www.michaelkingston.fi/files/pol ... ampler.osm

The quality difference is especially audible in the resonant high frequency end, but you can use anything besides that "moog" filter. There is slight attenuation near nyquist, but IMO it's a good trade off because of greatly increased quality.

play with the coefficients from musicdsp.org if you like, or expand or minimise the filter to 2-12th order version.

I'm using the 10th order version with 86dB stopband attenuation and steep slope.

Basically the stereo input gets packed, then upsampled as a stream, the dual output gets shuffled to 4 streams for efficient parallel processing. It processes stereo 2x oversampled stream in one pass. I can't even believe how efficient this is!

The same happens reversed on the downsampling end.

The whole resampling process only takes about 1.5% CPU on my system, and that's just about the only CPU hit since all oversampled material can be (with few exceptions) parallel processed.

In fact, one could easily create a 4x cascaded oversampler for mono material with this exact same system. Just keep in mind the extra allpass filters for 4X oversampling will use at least another 2-3% CPU. Again, the actual audio processing will be "free" as streams.

I know there are plenty of better coders here than me, so would someone please take a look at the allpass filters, if it would be possible to optimise them further? Both filters A and B are identical apart from coefficients and the same ones are used in both up and downsampling.

Kingston · 18 May 2006

OK I've now implemented this in project and the results are as expected: perfect!

With correct stream routing and packing and unpacking within the project, the overall CPU hit of the resampling is only about 1.5-2%!

Basically my plugin uses about 3.8% CPU in normal mode with all features on, and the 2X oversampling takes it up to only 5.5-6%. Less than half more CPU, even though there's clearly more than double the calculations done since every stereo audio is kind of processed twice. (SSE rocks!)

This makes it extremely useful in my projects and ups the soundquality yet another crucial notch.

This all with the unoptimised allpassfilters!

😍

The coolest thing is that it would likely be impossible to optimise the stream based scheme further in even C++ SSE enabled compilers. I'm finally starting to see how the initial processing overhead of synthmaker can be taken advantage of, making it fully competent with hand optimised SSE code.

Anyway,

like I said before, I could use some help with optimising the filters, but also on optimising the oversampling performance.

How to create an even better oversampler? Less high end attenuation? Minimum phase implementation?

Right now there's some of high end attenuation, maybe more than I'd like. Or maybe I did something wrong afterall?

rl · 18 May 2006

This is interesting. How did you design the filter? Or is this a port of the music-dsp example?

Kingston · 18 May 2006

It's direct a port of the musicdsp filter. You can easily modify it to use 2-12th order allpass cascades.

There in the comments section I found a link to the original research publication on the correct cascade implementation and coefficient design, as it wasn't perfectly clear just looking at the C++ code. AFAIK the c++ code only does upsampling, so I had a bit of difficulty figuring out the correct downsampling.

It's also relatively easy to cascade it further to 4X oversampling.

I think I did it right, but feel free to suggest corrections and synthmaker related optimisations.

At times like this I wish more and more there would be plugin latency reporting to host. My 10th order implementation introduces 10 samples of latency on the signal path.. It's a problem in mixing unfortunately.

rl · 19 May 2006

When feeding a pulse into the schematic and recording the pulse response, there is a discontinuity in the FFT-spectrum at 10kHz. The phase response looks chaotic over the biggest part of the frequency range. Maybe it's caused by calculation inaccuracies. I tested with some filters (Moog, SVF). In any case the non-oversampled version looks cleaner 🙁 . I will check this in more detail later since I really like this oversampler (stereo channels with 2x oversampling in one SSE channel).

Kingston · 19 May 2006

Yeah I've noticed problems with it as well.

Keep in mind that I've at all times I said that I'm not sure if the filter was ever working correctly. Maybe if you have the means you could review the filternetwork as well, using the links provided in this thread?

I'm doing my own testing right now as well, to get to the bottom of this.

The funny thing is, it already *sounds* better than the non-oversampled version on all material I've tested on.

In case I did something wrong and will later correct it, I've saved this possibly "wrong" version for future use as a great sounding process.

Note that the oversampler is linearphase, and it should result in both pre and post "ringing" on the oversampled material, which I've already tested it does. Also note that the attenuation after 10khz could be the normal attenuation the resampling filters introduce, although right now it seems quite a lot.

Just not sure if it does it *right*. 😊

Here's the original implementation and research by the way: http://users.cscs.wmin.ac.uk/~krukowa/pdf/Paper11.pdf

rl · 19 May 2006

Great work, Kingston. From what I see your SM code implementation looks right. Maybe it has to do with float vs. double accuracy since the example code from musicdsp uses doubles. I'm currently looking at Laurent de Soras' HIIR implementation (http://ldesoras.free.fr/prod.html ) which seems to use the same algorithm with floats.

Kingston · 20 May 2006

just a minor update. the code I posted was definitely wrong (well, not the actual code, but the routing). Not much, but wrong enough to introduce the 3dB lowpass point at half nyquist. Sorry about that. I'll post a working version later today. 😉

Nothing wrong with 32bit calculation accuracy in this. This should work just fine even as 16bit integer implementation, which is what the original research did with their specialised chip (there's a whole paragraph about IIRC).

Anyway, back to work. 😀

[edit] ok here's a second order version (I had to simplify to make it more clear) http://www.michaelkingston.fi/files/pol ... mpler4.osm

It doesn't seem to introduce the attenuation. Still not sure wether it's right.

sambean · 20 May 2006

Hey chaps..

Not sure if this is working I'm afraid, I may be wrong though.. I tried to use it to ringmod a sine wave, while analysing the output with a mono FFT. It seems the sine wave is still been aliased back after f-input is greater than fs/4. with proper interpolation this shouldn't be happening, or at least it should be attenuated.. It seems to be doing something to the response though..

I reccomend not using the pack/unpack optimizing until the problem is nailed, this way it's easier to hook up the analyser and spot the aliasing.

I'm going to have a wee peep to see if I can get it working..

sambean · 20 May 2006

OK, think I have it working now, but it's unable to take advantage of any packing.. It works by doing the oversampling the old skool way, by just running each of the half band filters twice in one sample period.. To do this I had to combine all the various allpass codes into one module and and run through the code twice, so the codes quite big...

For the example there's a saw wave going through a tanh waveshaper, you can see that there is virtually no aliasing until you get up to the really high frequencies..

I think this is going to be pretty intensive cpu wise, but the results look pretty good..

http://magalov.myweb.hinet.net/toolkit/IIRPP2.osm

Cheers

Sam

Kingston · 20 May 2006

well what ever even that 10th order version of mine is doing is already sounding better.

Seems to mearure as linearphase response of *something* when it's fed impulses or squarewaves.

Kingston · 20 May 2006

sambean wrote: OK, think I have it working now, but it's unable to take advantage of any packing.. It works by doing the oversampling the old skool way, by just running each of the half band filters twice in one sample period.. To do this I had to combine all the various allpass codes into one module and and run through the code twice, so the codes quite big...

For the example there's a saw wave going through a tanh waveshaper, you can see that there is virtually no aliasing until you get up to the really high frequencies..

I think this is going to be pretty intensive cpu wise, but the results look pretty good..

http://magalov.myweb.hinet.net/toolkit/IIRPP2.osm

Cheers

Sam

Nice,

but I don't see why you couldn't use any packing? it still streams in two doesn't it? just pack before input, and unpack the outputs, justlike I did with my filter (with extra shuffle for stereo)?

Here's a packed version. Seems to be working just fine. (and yes indeed there's phase oddness after nyquist/4 in my own version). Compare.

http://www.michaelkingston.fi/files/pol ... er5sam.osm

The CPU hit is not *that* bad. But surely there are better ways to resample than this?

I just tested this in a project and the CPU hit isn't bad at all.

I know for a fact you did actually manage a correct filter this time, since it has a latency of 5 samples, the correct number for a 10th order linearphase filter.

sambean · 20 May 2006

Here's a way of taking it a little further.. But it adds a one sample delay to the process unfortunately, which I take it you were trying to avoid, by using IIR..

With your stereo version only 2 channels of the SSE were being used for each oversampling section, so what I've done is got rid of the second oversampling section, fedback the output from the 'Effect' into the bottom 2 channels of the first oversampling section, so both upsampling and downsampling filtering is being done in the same block, with an unavoidable sample delay being added before the downsampling.. Now all channels of SSE are being utilized, nothing wasted..

http://magalov.myweb.hinet.net/toolkit/IIRPP3.osm

Kingston · 20 May 2006

Whoha!

Cheers sambean. What a great optimisation.

latency is not a problem at all with any process I use, when it stays somewhere below 100 samples.

and hopefully we will be getting a VST plugin latency reporting to hosts for synthmaker somewhere in the future.

Maybe this oversampling thing is a great incentive for malc to implement it? 😍

Oh and now I *finally* get the filtercascade. That's the bit I was struggling with and ended up with a half working solution. Luckily I get to keep the CPU usage of my original filter. It seems about equal thanks to your new SSE shuffle.

sambean · 22 May 2006

Hey Kingston, some time ago on kvr you posted a rather tasty sounding stereo field extraction example that used the wavelet transform. I seem to remember you saying that it relied on the fft to perform the transform and was thus out-of-scope for a real-time application..

Just thinking aloud, what we have here is a high quality half-band filter, which can be easily configured to give you the upper and lower frequency band.. http://users.cscs.wmin.ac.uk/~krukowa/pdf/Paper27.pdf (see fig a) Would it not be possible to do an alias free wavelet decomposition using this half-band filter..

eg Split Input into highband and lowband Downsample lowband by 2 using hop (alias free) Split Lowband into highband and lowband Downsample lowband by 2 using hop etc etc

[Do the wavelet processing]

Reconstruct signal using reverse process...

Whether it would be fisable or not would depend on the number of bands.. Also, I'm not sure if this method would be equivalent to the method you used..

i may have a try this evening at getting the wavelet transform and reconstructing the signal.. I'd tried this before with the FIR interpolator filter but had given up as the cutoff wasn't strong enough and the combined latencies became a hassle to manage..

Any thoughts?

Sam

Kingston · 22 May 2006

You know, I've been thinking about doing a synthmaker version of the stereo extraction gadget of mine. The first problem that comes to mind is that we would seriously need the isoffline() for it to work accurately enough. If it uses good sounding true wavelets, it might be feasible as a realtime process.

You see, I used a FFT based pseudo-wavelet approach. Yes indeed it looks like this paper might help with realtime implementation.

SSE will also heavily decrease the burden on the CPU. Obviously my Matlab script wasn't using SSE in the most optimal way, if at all.

Wether the signal is represented in wavelet, or FFT form, every single bin will have to be separately compared for the stereo image extraction. Basically for 256 sample FFT that means 256 comparisons in frequency domain, with the obvious FFT artifacts. Do that with optimal pseudo-wavelet FFT and the comparisons will increase exponentially.

SSE wavelets might be the key to higher performance and low CPU. I'll be interested to see your signal de- and reconstruction with this.

Don't worry about latency by the way, at least with the stereo extraction it's not an issue, since it's hardly a realtime process.

[edit]

It looks like this paper might be the answer to the realtime implementation of the stereo extraction.

1. create wavelet filterbank using the method in the paper. 2. depending on implementation feasibility, creating one big wavelet bank might work, or maybe simply use (for example) 16 separate STFT (short-time FFT) comparison processors on each wavelet bank process those, or process the wavelet banks directly.

3. recompose the signal.

It might just work, but I'll have to first see how to use the stereo extraction directly on the wavelet filter banks. It it proves impossible, at least STFT is fast, and in this case reasonably accurate.

Oh and having the same separation of the original resampler allpass filter on the bands is quite awesome! Great find sambean. Of course it would still be very CPU intensive due to the sheer number of the halfband filters needed for wavelet construction. Also having to do STFT using the code component is a chore and a half.

Has anyone already done FFT/STFT in synthmaker?

rl · 22 May 2006

Kingston wrote: Has anyone already done FFT/STFT in synthmaker?

I did some experiments, but gave up according to:

reason #1: no iFFT (constructing a 512 bin sine-oscillator bank in S|M is no fun and makes no sense CPU-wise).

@Developers: the FFT implementations I know are symmetrical and can be controlled by a single parameter to do FFT or iFFT. So maybe it's not too expensive to provide a iFFT version of S|M's FFT module.

other reasons: I'm unsure how the FFT module internally works in terms of window size, windowing function, phase correction (My guess is: The input graph data is filled up to the next 2^n boundary (but not above 1024) and no window function is applied).

For STFT obviously some overlap logic has to be build using sample accurate triggers for the FFT process and I do not yet know how to build this using the available module set.

Kingston · 22 May 2006

Well for this particular task it's not as impossible as it sounds. Something like 16sample STFT/iSTFT code module could be constructed.

Not *quite* as evil as having to do 4096 sample FFT with windowing in SM code module.

no function calls...

sambean · 22 May 2006

Green iFFT would be very nice indeed.. As for stream FFT/iFFT that would be nice too along with some sort of block processing and series oversampling... but i can wait...

Here's the first take on the Half Band Wavelet Transform.. It seems to work fine but there are problems.. As the IIR filters that are being used are such high order, they have quite long impulse responses, and this is interfering with the sound at the transition bands in the form of ringing, this is more prominent at the lower registers.. It may be fixable by adjusting the phases, which have probably been screwed up through the process.. But overall the principle is working, try deleting some of the connections between the modules and you should see bands dissapearing from the frequency display.. If you can import a sample it may sound better than the oscillator as the ringing seems to only affect when sound starts and stops or a harmonic lies exactly on the cusp of one of the bands..

The efficieny could probably be improved lots, but everything is laid out as it is for clarity..

Here it is... http://magalov.myweb.hinet.net/toolkit/IIRPPWavelet.osm

Sam

sambean · 22 May 2006

A quick point about processing in between the transform pairs:

Each consecutive link is oversampled by a factor of 2, so any processing you do in a code module can be put within a hop(), hop(2) for the first band, hop(4) for the next etc etc, and hop (32) for the last 2..

This not only improves cpu but means any spectral processing you do applies to the whole frequency band.. So if you do write a short 16 bin FFT routine, you can pop it inside a hop, and this will give you 16 bins per frequency band..

Just need to see if it's possible to eliminate the ringing, as I think this could be pretty useful...

Sam

Kingston · 22 May 2006

odd ringing indeed.

I did some testing and it looks (and sounds) as if the phases of the filters are reconstructed wrong.

Right now it kind of sounds like the high bands are arriving first, making for a weird group/band delay series.

I'll see what can be done with the stereo extractor. It might not be such a simple implementation, without having FFTs written for synthmaker already. I mean I can *use* FFTs, but writing my own routine..*shrudder*

But thanks to the downsampling I can at least put 16bin, or similar STFTs to every single wavelet band. How this differs from my pseudo wavelet approach, is that I had to process every single one of the bands on full rate... quite the waste!

So instead of processing the lowband with a small 16bin FFT, I had to process it with 4096, and only the highest band could use the small FFT size. I feel stupid now, for not using wavelets in the first place, although when I originally did it realtime use was hardly any concern to me.

Kingston · 22 May 2006

Just a though as I haven't got the time do fix it myself,

when reconstructing the bands, the unit delay that accumulates through the reconstruction in each successive band is most definitely the reason for the ringing. The low bands are now always being delayed one unit further than the high bands.

That's an easy fix luckily. 😀

sambean · 23 May 2006

Here's a link to a 16-point FFT:

http://www.csua.berkeley.edu/~emin/sour ... index.html

A problem implementing it may be that nested hops are required, which SM doesnt support ATM eg

hop(DOWNSAMPLE_RATE) { inbuffer[n] = in

hop(16) { ... FFT code ... } n = n + 1; n = (n<16)&n }

A way round this may be to write the main code with no downsampling using the code module, and then using asm incorporate the downsample hop..

Kingston · 23 May 2006

thanks. that looks fairly straightforward. I'll have a go at it when I find the time. Too bad with all these exciting new developments I have to leave the studio for a week and can't work.

cAPSLOCK · 12 Jul 2006

For what it's worth that code now lives here: http://web.mit.edu/~emin/www//source_co ... index.html

I sure would like to see fft and ifft implemented in SM (green ones).