How do you calculate the group delay in samples of a Biquad Filter at a given frequency from the coefficients?

If I have a BiQuad all-pass filter with a given set of coefficients such as:

g: 0.903228
a1: -1.88142
a2: 0.903228
B1: -2.08299
B2: 1.10714

And my operation per sample to process the filter is:

inputGained_2 = inputGained_1;
inputGained_1 = inputGained;
inputGained = input * g;
output_2 = output_1;
output_1 = output;

output = inputGained + (B1 * inputGained_1) + (B2 * inputGained_2) - (a1 * output_1) - (a2 * output_2);

return output;

And this is running at 44.1 kHz: What equation could I used to calculate the group delay at a given frequency like say 82 hz?

Thanks.

• octave.sourceforge.io/signal/function/grpdelay.html May 13 '20 at 7:13
• Thanks. I see this function: [g, f] = grpdelay(b,a,f,Fs) evaluates the group delay at frequencies f (in Hz). I obviously would have f and Fs. Can you share a few lines of Matlab/Octave code that would let me enter a and b in terms of my variables a1, a2, B1, B2, and g and calculate this? I'm sure it would be simple if you know Matlab/Octave but I'm not good with Matlab/Octave. I've only used it twice before. Also is there any way to look up how the grpdelay(b,a,f,Fs) function works inside Matlab/Octave so I can copy it to C++ for my needs? Thanks.
– mike
May 13 '20 at 7:36
• Wasn't there couple examples on that page I linked? Octave sources can be found from here hg.savannah.gnu.org/hgweb/octave/file/332e644726f9 (check the licensing info - hg.savannah.gnu.org/hgweb/octave/file/332e644726f9/COPYING ). May 13 '20 at 9:44
• Are you sure you need the group delay and not the phase delay? The delay of a sinusoid at frequency $f_0$ is given by the phase delay at $f_0$, not by the group delay. May 13 '20 at 11:18
• I don't know then @MattL. I'm dealing with a situation where I have integrated a series of cascaded all-passes in a Karplus-Strong feedback-delay waveguide synthesis loop. These are creating a certain number of samples of delay in the all-pass cascade at the frequency that the Karplus-strong loop is tuned for. I need to subtract the number of samples delay of the all-passes from the waveguide model or it will be detuned by the extra delay. I need to calculate how many samples to compensate for at that fundamental frequency. Which one do I want in that case? Phase delay or group delay? Thanks.
– mike
May 13 '20 at 15:46

First of all, your filter coefficients are given by

b = g * [1, b1, b2]
a = [1,a1,a2]

And second, I believe that it is the phase delay that you're interested in, not the group delay. The phase delay at frequency $$\omega_0>0$$ is given by

$$\tau_p(\omega_0)=-\frac{\phi(\omega_0)}{\omega_0}\tag{1}$$

where $$\phi(\omega)$$ is the (unwrapped) phase of the filter's frequency response. The frequency $$\omega$$ is the normalized frequency in radians, i.e.,

$$\omega=2\pi\frac{f}{f_s}\tag{2}$$

where $$f$$ is the frequency in Hertz, and $$f_s$$ is the sampling frequency (in Hertz).

So you just need to evaluate the filter's frequency response at the frequency of interest, then compute its phase, and then compute the phase delay from $$(1)$$.

• Thanks Matt L! That was very helpful for clarifying how to enter the coefficients. I did not know what it wanted for numerator or denominator in Matlab. Again, if you did not see my comment above, I have integrated a series of cascaded all-passes in a Karplus-Strong feedback-delay waveguide synthesis loop. These are creating a certain number of samples of delay in the all-pass cascade at the frequency that the Karplus-strong loop is tuned for. I need to subtract the number of samples delay of the all-passes from the waveguide model or it will be detuned. Is that phase delay or group delay?
– mike
May 13 '20 at 16:25
• @mike: I think you should be using the phase delay. May 13 '20 at 16:54
1. Pick two different frequencies slightly above and below your target frequency. The choice depends on what numerical precision you have available and what your sample rate is.
2. Calculate the z-transform at both frequencies
3. Calculate the phase at both points, make sure to check for "phase wrapping"
4. Subtract the phases and divide by the frequency difference
5. Scale properly and multiply with $$-1$$