Comb Filter and delay compensation

[Mirko R.]

Hello,

I am developing a comb filter. The structure is simple. There is just a buffer and a feedback. The filter (comb) frequenzy will set the time for the delay buffer. It works fine. But if I add a simple 6db LP (damping) to the feedback, the comb frequency is not 100% correct anymore. the more I decrease the damping frequency, the more the comb frequency will go down. do someone know how I can compensate this?

Here is the code from the simple LP. If lpFraction is set to 1.0, everything works fine and there are no frequency-shiftings. The damping is off. When I decrease the value, more damping is active but also the frequency will go down.

Code: Select all

//Feedback LP Filter
// lpFraction between 0-1
out = out + lpFraction * (memoryOut - out);
memoryOut = out;

And here is the code to calculate the delay time:

Code: Select all

const double delay = 1.0 / combFrequency * SAMPLERATE;

[Max M.]

The LP (as well as any other filter) has its own (frequency dependent) delay (usually referred to as "group delay" which is basically just another representation of the filter's phase-shift). Thus the resulting delays of your feedback taps become multiples of [delay-delay + filter-delay].
A compensation would be just a matter of calculating the group delay of the LP filter and adjusting the main delay accordingly (I don't recall the exact formula but you'll find it by searching for corresponding keywords). The only problem is that the delay of the filter is frequency dependent so you'll have to choose the frequency the whole module will have the "exact" delay at (obviously the resulting thing is no longer a strict comb filter anymore and its peaks/valleys are no longer at multiples of the base frequency no matter how and where you compensate).

[Mirko R.]

Thanks for your reply. Okay I think I know what I have to do. I need the phase shifting from the filter. For Example:

comb frequency = 2000 hz
damping frequency = 1000 hz

delaytime = 1.0 / combFrequency * 44100 = 22.5

Let's suppose the phase shifting at this damping frequency will be 45,° than I have to calculate it like this:

degFact = 45°/ 360° = 0.125

delaytimeFinal = delaytime - degFact * delaytime = 22.5 - 0.125 * 22.5 = 19.6875

Is that correct? If yes I would like to know how I can get the phase from my 6 db simple filter.

[Max M.]

I think it's much more easy to compute the delay in samples directly w/o any phase shift.

Assuming your filter (it's basically a leaky integrator) direct form is:

Code: Select all

b = [1-k 0]; a = [1 k];

where k is your lpFraction (most likely inverted i.e.: k = 1 - lpFraction), the group delay in samples would be:

Code: Select all

D = k*(cos(w) - k)/(1 - k*(2*cos(w) - k));

where w is the normalized frequency of interest (in radians), i.e.:

Code: Select all

w = 2*pi*f/fs;

where f - frequency of interest
and fs - samplerate
(both in whatever units, e.g. Hz)

I suspect the equation can be simplified (via some trigonometric trickery or so) as it looks a bit complicated for such a simple (just one pole) filter, but this is what I could get at quick glance (right now I don't have any heavy math tools nearby and only tested the equation against octave's grpdelay function).

[Mirko R.]

Great and thanks a lot. Of course, I can take your stuff now and use it but that's not what I really want. I would like to understand what you did there. How should I get the formular for the group delay?

For example when I'm using a filter with the coefficients:

Code: Select all

a0, a1, a2
b1, b2

My transfer function will look like this:

Code: Select all

H(z) = Y/X =  (a0 + a1*z^-1 + a2*z^-2) / (1 + b1*z^-1 + b2*z^-2)

As time domain:

Code: Select all

y[n] = a0*x[n] + a1*x[n-1] + a2*x[n-2] - b1*y[n-1] - b2*y[n-2]

How can I get the group delay with these Informations. Yes I can google it and I found some struff but it looks so academic that I don't understand what is happening.

[kamalmanzukie]

for the things i use comb filters for, i've gotten better results with FIR, and better still using only delays with only an integer number of samples in the buffer

[Max M.]

How should I get the formula for the group delay?

Something like this.

but it looks so academic that I don't understand what is happening.

Well, it is pretty complicated (especially since beside usual `z` -> `e^jw` -> `cos(w) + j*sin(w)` and subsequent trigonometric substitutions you do for magnitude/phase equations, group-delay also involves differentiating and so on. So I'm afraid I would not be able to provide a good step-by-step guide on this).

So... since the H(z) in your recent post is just a generic biquad, I'd rather suggest you to search for a ready-to-use "Biquad Group Delay" equation.
(Unless the goal is to learn all that complex math of course

[Mirko R.]

for the things i use comb filters for, i've gotten better results with FIR, and better still using only delays with only an integer number of samples in the buffer

Why FIR filters and which kind of? I think they need a lot of CPU?!
Delay-Times as integer is in my case not a good idea because the comb frequencies are then not flexible enough. They are steppy. Try to use a float value (or double but float is enough) and interpolate it like this:

Code: Select all

// write you input into the buffer
buffer[write] = in;

// calc delay time
float delayTime = 1.0 / combFrequency * SAMPLERATE;

// calc the read index as float
float index = write - delayTime;

// get the integer indices around the float index
int indexIntA = (int)index;
int indexIntB = indexIntA + 1;

// get the interpolation factor
float interpolationFactor = index - indexIntA;

// read the buffer interpolated
float out = buffer[indexIntA] + interpolationFactor * (buffer[indexIntB] - buffer[indexIntA]);

// increment the write pointer
write++;

The result is fine. try it and don't forget to check the validity of the indices.

Code: Select all

if (index < 0)
   index += buffersize;
....

--------------------------------------------------------------------------------

Well, it is pretty complicated (especially since beside usual `z` -> `e^jw` -> `cos(w) + j*sin(w)` and subsequent trigonometric substitutions you do for magnitude/phase equations, group-delay also involves differentiating and so on. So I'm afraid I would not be able to provide a good step-by-step guide on this).

Of course it is but i know how to handle complex number, poles and zeros, but to understand the mathematics I need examples. The most stuff in the internet is the pure theory without any examples. That's my problem. But thanks anyway. I will try to find a solution. Maybe someone could tell me how to get the phase response from a biquad. That could help me for the first time.

Thanks and,.... sry for my bad english.

[sonigen]

My transfer function will look like this:

Code: Select all

H(z) = Y/X =  (a0 + a1*z^-1 + a2*z^-2) / (1 + b1*z^-1 + b2*z^-2)

if you multiply the transfer function by

(i) z^2/z^2

H(z) = (a0 + a1*z^-1 + a2*z^-2)*z^2 / (1 + b1*z^-1 + b2*z^-2)*z^2

H(z) = (a0*z^2 + a1*z + a2) / (z^2 + b1*z + b2)

That gets rid of the negative powers of z, so next you set z to the complex number...

z = cos(2*PI*f/sr) + i.sin(2*PI*f/sr)

where f is the frequency you're analyzing, and sr is the sample rate.

You can now evaluate H(z), you can calculate with complex number library, or you can sub z into H(z) and simplify. Not sure how much it'll simplify though. The phase & magnitude of the complex result is the phase & magnitude response of the filter.

You just need to unwrap the phase to get group delay. Which AFAIK theres no analytical way, other than just evaluating a few frequencies and searching for any wrapping and then handling it however you need. You know delay at DC is 0, so work from there kinda thing.

Thats all IIRC, im a long time out of practice tbh.

[mystran]

For tuning a specific harmonic (and you can't really tune more than one) you want to use phase delay rather than group delay. Specifically, group delay is the derivative of the phase response with respect to frequency and it's used to predict the delay on the amplitude envelope. That's not really relevant in terms of tuning though, since the resonant peaks always happen at those frequencies where the phase is a multiple of 2*pi.

Calculating phase delay is really simple if your filter is 1st or 2nd order, just calculate the phase response (between 0 and 2*pi), multiply by wavelength and divide by 2*pi. For higher order filters you might have to unwrap the phase, so the most obvious thing to do in most cases is to calculate the delay for each biquad section (or 1st order section) separately and just add them together (which obviously requires factorisation, but often times you're building more complex filters from simple sections anyway, in which case it's a non-issue).

[sonigen]

Yes of course phase delay not group delay, oops.

[kamalmanzukie]

for the things i use comb filters for, i've gotten better results with FIR, and better still using only delays with only an integer number of samples in the buffer

Why FIR filters and which kind of? I think they need a lot of CPU?!
Delay-Times as integer is in my case not a good idea because the comb frequencies are then not flexible enough. They are steppy. Try to use a float value (or double but float is enough) and interpolate it like this:

in this case its just subtracting the delayed signal from the signal, the delay corresponds to the cutoff frequency

i only say it as advice because its my only advice from the one time i ever used comb filters: for attempting to null (or close to) a periodic sound within a signal. the cutoff frequency was the same as the period of the signal i was trying to remove. it actually worked better without feedback, and surprisingly enough it was more effective using integer number of samples than what i assume was some kind of quadratic interpolation, even though the pitch of the signal to kill was not static

admittedly, its not much help for someone who uses comb filters for non deranged reasons. the only thing to enter my head is 'a filter that filters at harmonics? gr8 i will now try removing a single note from a full recording'. but if you ever wanted to, thats what works best

using integer samples, the loss of frequency resolution would mainly be lower frequencies, no?

[kamalmanzukie]

aren't comb filters known for having the harmonics go wonky and out of tune? would having a regular old unit delay without any 'zdf' type of fix contribute to that

[mystran]

aren't comb filters known for having the harmonics go wonky and out of tune? would having a regular old unit delay without any 'zdf' type of fix contribute to that

If we assume there are no damping filters in the loop, then a pure comb should really give you perfectly harmonic partials as long as the interpolator doesn't introduce significant phase non-linearity.

Using all-pass interpolation, you'd expect some detuning of the higher harmonics as the interpolation only really works at low frequencies and still gives you integer number of samples at Nyquist. Using FIR interpolation instead can avoid this issue, at the cost of some high-frequency losses and/or minimum delay time as better kernels need more latency which puts a lower limit on the delay time. Both of these issues can be mitigated by some oversampling though.

Beyond interpolation, minimum-phase damping will introduce some non-linear phase-shift that will detune harmonics somewhat, but this is really the same thing that happens with any physical system as well and if the damping filters are sensible it should usually sound fairly natural rather than "wonky" as such. In extreme cases, if the damping is strong enough to significantly affect the phase at the fundamental it can affect the perceived tuning, so you might still want to compensate.

Intentionally using filters that are arbitrarily more broken usually doesn't help anything.