linear phase oversampling

[dasdeck]

Hi!

I'm looking to improve my oversampling algorithms. for now I either use a pretty long FIR (257) or the polyphase filter from musicDSP (which is not linear phase of course). I looked a little bit into other implementations and found out that for example the Drop by Cytomic does all of its oversampling linear phase (at least in the audible/visible range) with a delay of 52 samples! and the aliasing is already very low at 8x oversampling, better then my approach at 63 samples delay.

I use three FIRs (x2,x4,x8), with each 37 taps. The FIRs are designed with a squared (sqrt()) blackman-harris window. The frequency response of my "outer" filter seems to be matching the Drop's. But I had to lower the cutoff frequency to about 20k to match an equal amount of aliasing rejection.

What I did in the past is to use one filter for the up sampling, meaning that I used a different cutoff (nyquist/4) and stuffed with three zeros instead of one. Since the filter is not as steep as it gets farther away from the nyquist, I suppose its better to split the oversampling into 2x cascades.

Right now i'm thinking of only making the first cascade linear phase and from there use any other lowpass filter that has a "no" delay and is mostly linear phase in the lower area of the passband.
Since the upper frequencies will get thrown away anyways, I should be ok with weird artefacts in the 44k+ area right?.

I have some ideas to tryout, but I'd like to have your input. Maybe I'm totally of track here, and whats going on in the Drop is totally trivial and I'm just missing it. Or maybe there are other, better ways of archiving what I want. Basically I want an oversampler with a transparent phase response and a lot of rejection at like 8x that is usable in realtime.

Cheers!

JM

[dasdeck]

I'm not quite sure what the name of the design method that I used is, but I'm pretty sure it is not raised cosine. I suppose its some variation of a sinc function. However, will there be much improvements in filter steepness compared to the current method with applied windowing? I suppose not? Is the "trick" somewhere else?

here is the method for the coefficients (only one half, they get "mirrored after i apply the window") I use.

Code: Select all

h[0] = 0.5;
    for(int i = 1; i < numTaps ; i++)
    {
        h[i] = (sin(  i * double_Pi * 0.5) / (i * double_Pi);
    }

then I apply a squared blackman-nuttall window

Code: Select all

const double a = 0.3636819;
    const double b = 0.4891775;
    const double c = 0.1365995;
    const double d = 0.0106411;
    const double pi = double_Pi;
    const double n1 = length*2-1;
    
    for(int i = 0 ; i < length ; ++i)
    {
        int n = i+length;
        const double a1 = (2 * pi * n) / n1;
        const double a2 = (4 * pi * n) / n1;
        const double a3 = (6 * pi * n) / n1;
        const double factor = a - b * cos (a1) + c * cos (a2) - d * cos (a3);
        data[i] *= sqrt(factor);
    }

the resulting plot in the frequency domain is this (pink is phase response, yellow frequency)

plot.png

Any ideas on this topic? maybe a reading suggestion?

cheers!

edit: I just realised that my method is effectively the sinc function, as freq0 is always 0 , thus the divisor is always zero. I guess i should stop talking and go test the raised cosine method then? At least then its in my codebase.

edit: I simplified the formula for the halfband coefficients.

JM

[dasdeck]

Ok, so the question boils down to this:

I use a 37 tap windowed-sinc filter cascade for up-sampling to x8 resulting in a delay of 63 samples. how to improve from there in terms of aliasing rejection vs. delay?
From what I've seen more rejection should be possible with even less delay (e.g., Cytomic's the Drop @ 53 samples and seemingly more rejection than my model ).

cheers!

jm

[sonigen]

Code: Select all

h[0] = 0.5;
    for(int i = 1; i < numTaps ; i++)
    {
        h[i] = (sin(  i * double_Pi * 0.5) / (i * double_Pi);
    }

That looks wrong to me. Firstly h[0] would be the center tap, at x=0, based on the way you calculate the other taps, but in that case it should be 1. If on the other hand you are mirroring but duplicating h[0], then the mapping of the other taps is wrong. Secondly the sinc cutoff is fixed at nyquist. When you design a windowed sinc filter it's the cutoff of the sinc that changes, and the window that truncates it to a finite length. So it should be something like this...

Code: Select all

for (int i = 0; i < numtaps i++)
{
    int x = i-numtaps/2;
    if (x == 0) fir[i] = 1; else fir[i] = sin(2*PI*i*fc/fs)/x;
}

;

then you apply a window function to it. Personally i think the kaiser window is the best bang for buck.

[2DaT]

Ok, so the question boils down to this:

I use a 37 tap windowed-sinc filter cascade for up-sampling to x8 resulting in a delay of 63 samples. how to improve from there in terms of aliasing rejection vs. delay?
From what I've seen more rejection should be possible with even less delay (e.g., Cytomic's the Drop @ 53 samples and seemingly more rejection than my model ).

cheers!

jm

Its always a sacrifice between length, high frequency passage , aliasing and phase response.

[dasdeck]

I know , but as I stated, the over-sampler in "the Drop" does crystal clear at 8x oversampling with a delay of 53 samples. with my current implementation I can not get there. I "measured" the anti aliasing performance by sine sweeping, and the difference to my implementation is obvious. of course it may be that his algorithm produces less harmonics to begin with :O. on the other the saturation's aliasing at no oversamplin in the Drop sounds just as bad as with a simple tanh approximation. I can get my saturator to not alias with my algorithm, but i have to set the cutoff frequency a little more below the nyquist. I wanted to know if there is a better way.

thanks

JM

[dasdeck]

Code: Select all

h[0] = 0.5;
    for(int i = 1; i < numTaps ; i++)
    {
        h[i] = (sin(  i * double_Pi * 0.5) / (i * double_Pi);
    }

That looks wrong to me. Firstly h[0] would be the center tap, at x=0, based on the way you calculate the other taps, but in that case it should be 1. If on the other hand you are mirroring but duplicating h[0], then the mapping of the other taps is wrong. Secondly the sinc cutoff is fixed at nyquist. When you design a windowed sinc filter it's the cutoff of the sinc that changes, and the window that truncates it to a finite length. So it should be something like this...

Code: Select all

for (int i = 0; i < numtaps i++)
{
    int x = i-numtaps/2;
    if (x == 0) fir[i] = 1; else fir[i] = sin(2*PI*i*fc/fs)/x;
}

;

then you apply a window function to it. Personally i think the kaiser window is the best bang for buck.

thanks !

I'll try the kaiser window again. It is variable isn't it? so what setup would i use for a good half-band? I will try your algorithm, but it is quite different. for example it does not fall of like mine does, so that is left to the window too.
also the mirroring I do takes place with h[0] as mirror axis, so no, h[0] does not get used twice.

Code: Select all

for (int j=0; j < tBuffer.getSize(); ++j)
    {
        if(j < buffer.getSize())
            t[j] = h[buffer.getSize()-j-1];
        else
            t[j] = h[j-buffer.getSize()+1];
    }

I end up with a filter that is symmetrical around the middle coefficient (linear phase) . as you can see from the plot I posted it actually works quite well.

You have a 1 in the middle, which makes sense. hmm, I'll give your version a try.

thanks

JM

[aciddose]

Wow those filters are terrible.

I'm too busy to graph mine at the moment, but I use a combination of windows.

First generate the sinc impulse and apply a kaiser-bessel with B = 6.84999~

Then take minimum phase and apply hann.

You will end up with an impulse with a far greater slope and far lower stop band by taking advantage of the ripples generated by either window.

The thing I don't know however is how you would modify the hann window in order to get the same effect when applying it linear phase immediately after the kaiser-bessel.

The slope is greater than 100db per octave, stop-band level is below -100db, delay is one sample although the performance is barely affected if you simply drop that sample giving a delay of zero. Of course, not linear-phase.

[dasdeck]

Wow those filters are terrible.

Care to elaborate ? This is a very simple implementation. but what is so terrible about it?

I'm too busy to graph mine at the moment, but I use a combination of windows.

First generate the sinc impulse and apply a kaiser-bessel with B = 6.84999~

been there, will try again (with that B value).

Then take minimum phase and apply hann.

you mean transforming the resulting filter?
I'm not sure on how to a apply a window on a non symmetrical FIR thought.
Did not do the transforming either, as I want to go linear phase

You will end up with an impulse with a far greater slope and far lower stop band by taking advantage of the ripples generated by either window.

The thing I don't know however is how you would modify the hann window in order to get the same effect when applying it linear phase immediately after the kaiser-bessel.

The slope is greater than 100db per octave, stop-band level is below -100db, delay is one sample although the performance is barely affected if you simply drop that sample giving a delay of zero. Of course, not linear-phase.

Sounds good, I will give it a try. But what is the phase response like? I generally prefer the sound of linear phase filters and I think I can hear the phase distortion form IIR over-samplers and I don't like it. Not sure I tested non-linear-phase FIRs thought.

Thanks a lot for you input, I hope you can clear it up some more.

Cheers!

jm

[aciddose]

I can just post the filter itself and you can try using it.

The transformation is:

Code: Select all

	DFT(&data);
	RealCepstrum(&data);
	IDFT(&data);

	PhaseCrop(&data);

	DFT(&data);
	MinimumPhase(&data);
	IDFT(&data);

with:

Code: Select all

void RealCepstrum(FTDATA *data)
{
	for (int i = 0; i < data->length; i++) {
		data->RF[i] = logf(cabs(data->RF[i], data->IF[i]));
		data->IF[i] = 0.0f;
	}
}

void PhaseCrop(FTDATA *data)
{
	data->RT[0] = 0.0f;

	for (int i = 0; i < data->length/2+1; i++) {
		data->RT[i] *= 2.0f;
	}

	for (int i = data->length/2+1; i < data->length; i++) {
		data->RT[i] = 0.0f;
	}
}

void MinimumPhase(FTDATA *data)
{
	for (int i = 0; i < data->length; i++) {
		cexp(data->RF[i], data->IF[i], &data->RF[i], &data->IF[i]);
	}
}

And extra stuff... hopefully have it all.

Code: Select all

void DFT(FTDATA *data)
{
	int k, i;

	const float _2pi = __PI * 2.0f;
	float a = 0.0f;
	float d = (2.0f / (float)data->length) * __PI;
	for (k = 0; k < data->length; k++)
	{
		data->RF[k] = 0.0f;
		data->IF[k] = 0.0f;
		float p = 0.0f;
		for (i = 0; i < data->length; i++)
		{
			float sr =  cosf(p);
			float si = -sinf(p);
			data->RF[k] += ((data->RT[i] * sr) - (data->IT[i] * si));
			data->IF[k] += ((data->RT[i] * si) + (data->IT[i] * sr));
			p += a;
			if (p > _2pi) {
				p -= _2pi;
			}
		}
		a += d;
	}
}

void IDFT(FTDATA *data)
{
	int k, i;

	const float _2pi = __PI * 2.0f;
	float i_n = 1.0f / (float)data->length;
	float a = 0.0f;
	float d = (2.0f / (float)data->length) * __PI;
	for (k = 0; k < data->length; k++)
	{
		data->RT[k] = 0.0f;
		data->IT[k] = 0.0f;
		float p = 0.0f;
		for (i = 0; i < data->length; i++)
		{
			float sr =  cosf(p);
			float si = -sinf(p);
			data->RT[k] += ((data->RF[i] * sr) + (data->IF[i] * si));
			data->IT[k] += ((data->RF[i] * si) - (data->IF[i] * sr));
			p += a;
			if (p > _2pi) {
				p -= _2pi;
			}
		}
		data->RT[k] *= i_n;
		data->IT[k] *= i_n;
		a += d;
	}
}

float cabs(float x, float y)
{
	return sqrtf(x*x + y*y);
}

void cexp(float x, float y, float *zx, float *zy)
{
	float expx = expf(x);
	*zx = expx * cosf(y);
	*zy = expx * sinf(y);
}

FTDATA is just a struct like:

Code: Select all

struct FTDATA
{
 float *RT;
 float *IT;
 float *RF;
 float *IF;
};

Horrible of course but I haven't bothered to clean up my code.

http://xhip.net/temp/sinc_b_minphase.wav
http://xhip.net/temp/blep_c_minblep.wav

Now the place where you can minimize the latency is with the integrated - minblep - version.

You can simply drop the first 16 samples without introducing too much error. This sample is 16x oversampled and typically you'd need a 1-sample delay to get the phase to line up.

With the minimum phase sinc sample however I believe if you drop 16 samples off the start you will introduce quite a lot of error.

In any case experiment and see if it works for you.

The point of what I posted was just to say that you can get a lot better results with higher slope and lower stopband with far shorter filters and lower latency by using a little trickery with the windows. If you modify the B coefficient I posted even slightly (0.00001) it will produce vastly inferior results. If you look very closely at the spectrum you get from the impulses you should notice that the combs/teeth of the stop band ripple can be made up line up just perfectly between the windows so as to minimize their level while maximizing slope.

I've tried to graph it but unfortunately none of the tools I have contain windows accurate enough to graph the filter itself without the important part being masked by sidebands/noise from the window.

Maybe using matlab or whatever other tool to set it up correctly might help.

As for how the filter sounds, because I'm using such an ultra-short filter (for processing cost) there is a significant amount of passband "droop" near cutoff, but I find this can be cured if desired with a tiny bit of high-shelf EQ. Generally only something you'd do to the end product, although it may also be possible to further transform the filter itself to get this result without compromising the slope/stopband, when I tried I found it wasn't worth the effort.

[dasdeck]

Thanks a lot aciddose

Cheers!

jm

[mystran]

One quick note: if you're going for higher than 2x, there's not necessarily much point in doing cascades with FIRs, as you can design for more branches directly instead. This way you can also do 3x or 5x or 7x or whatever else you might find useful and you only need to run one pass (=convolve the input with each of the branches, then interleave the outputs for up-sampling, or convolve the input with the respective branch and sum for down-sampling).

[aciddose]

The whole layered system doesn't make much sense when you can do a single convolution instead.

The convolution may be more expensive, but if you integrate ahead of time (blep) the convolution is transformed into a simple addition.

You can also drop any samples with a delta lower than some limit... say you want to allow aliasing noise up to -80db, simply skip any delta lower than 1/10000.


See http://sourceforge.net/projects/protracker/ for a naive example of this in action.

BTW: somebody needs to contribute an accurate sallen-key filter to the project. Those I've been able to find online so far with claims they are sallen-key don't actually provide a sallen-key implementation, for example independent gain/cutoff/Q. The implementation itself should be super simple but I don't have the resources to do it myself.

[sonigen]

I'll try the kaiser window again. It is variable isn't it? so what setup would i use for a good half-band? I will try your algorithm, but it is quite different. for example it does not fall of like mine does, so that is left to the window too.
also the mirroring I do takes place with h[0] as mirror axis, so no, h[0] does not get used twice.

h[0] as mirror axis makes sense but i think it should be 1.

I dont know what setup exactly, you have to make the decision yourself, how much rolloff you can live with, how much rejection you need, etc..

To be honest I use a program called ScopeFIR to design what FIRs i need, but, i bought it about 6 years ago when it was only about 70$, it's now $500 or something insane, but IIRC the demo will let you design filters it just wont let you export the coeficients. So you could try out the various window functions and sinc cutoffs, parameters etc.. to see what they give you. And then write your own code to generate the coefs.

[dasdeck]

Thanks a lot for your input guys!

I know have a much better understanding of what "variables" to tweak to get closer to the results I want. I find these concept surprisingly easy to understand when there are no mathematical proofs all over the place in the explanation. I'm already very close to what i'm aiming for.

So thanks again for your precious time.
Let me know if I can return the favour some how.
And just PM me if you want to see the end results of your effort

Cheers!

jm

edit: i feel confused again but nevermind , I'll get it at some point. my oversampler has 63 samples delay, the frequency response is like flat until 18k, the cpu consumption is totally within reason and the rejection sufficient. I'll try to learn more about these concepts thought.