Trying to demonstrate group delay but not seeing any in dsp.BiquadFilter.

3 views (last 30 days)
Richard Keniuk
Richard Keniuk on 14 Dec 2021
Edited: Paul on 16 Jan 2022
In the following code I am using part of the parametric equalizer example to show group delay in signals. I have generated a 2 channel signal of 40Hz @ -3.02dB FS on channel 1 and a signal of 200Hz @ -3.02dB FS on channel 2. In the plot ax1 is the magnitude response from 20-300Hz, and ax2 is the group delay over the same frequency range. In ax3 is the input signal, and ax4 shows the output signal after the filter has been applied. You can see the 3dB increase in the 40Hz signal. But... I don't see the ~244 sample delay in the 40Hz output. If I look closely at the 200Hz signal I can see the approximately -1 sample delay in the output.
format compact;
Fs = 44.1e3;
Samples = [0:10*Fs];
yInput(Samples+1,1) = 1/sqrt(2)*sin(2*pi*40*Samples/Fs);
yInput(Samples+1,2) = 1/sqrt(2)*sin(2*pi*200*Samples/Fs);
F0 = 40; N = 2; Q1 = 2; G = 3;
Wo = F0/(Fs/2); BW1 = Wo/Q1;
[B1,A1] = designParamEQ(N,G,Wo,BW1);
BQ1 = dsp.BiquadFilter('SOSMatrix',[B1.',[1,A1.']]);
% hfvt = fvtool(BQ1,'Fs',Fs,'FrequencyScale','Log');
% set(hfvt, 'OverlayedAnalysis', 'grpdelay');
yOutput = BQ1(yInput);
% determine magnitude and group delay
[h,fh] = freqz(BQ1,Fs,'whole',Fs);
[d,fd] = grpdelay(BQ1,Fs,'whole',Fs);
MagnitudeAt40Hz = 20*log10(abs(h(41))), DelayAt40Hz = d(41)
MagnitudeAt40Hz = 3.0000
DelayAt40Hz = 243.6285
MagnitudeAt200Hz = 20*log10(abs(h(201))), DelayAt200Hz = d(201)
MagnitudeAt200Hz = 0.0328
DelayAt200Hz = -1.3291
% plot
ax1 = nexttile; semilogx(ax1,fh(21:301),20*log10(abs(h(21:301))));
xlim([20 300]); xlabel('Freq (Hz)'); ylabel('Magnitude'); grid on;
ax2 = nexttile; semilogx(ax2,fd(21:301),d(21:301));
xlim([20 300]); xlabel('Freq (Hz)'); ylabel('Group Delay'); grid on;
ax3 = nexttile; plot(yInput(43350:44850,:));
xlim([0 1500]); xticks(0:250:1500); xlabel('Time (samples)'); ylabel('Input Waveform'); grid on;
ax4 = nexttile; plot(yOutput(43350:44850,:));
xlim([0 1500]); xticks(0:250:1500); xlabel('Time (samples)'); ylabel('Output Waveform'); grid on;
What am I doing wrong here? Shouldn't the blue trace be shifted 244 samples to the right in the Output Waveform?
Paul on 15 Jan 2022
The peak of the step response does very closely approximate (maybe it's exact? IDK) the group delay for this example, but I'm not so sure that would be true in general. I suspect that the quality of that approximation is going to be very dependent on the form of the filter. For example, let's compare BQ1 and BQ2 = BQ1^2
Fs = 44.1e3;
F0 = 40; N = 2; Q1 = 2; G = 3;
Wo = F0/(Fs/2); BW1 = Wo/Q1;
[B1,A1] = designParamEQ(N,G,Wo,BW1);
BQ1 = dsp.BiquadFilter('SOSMatrix',[B1.',[1,A1.']]);
BQ2 = dsp.BiquadFilter('SOSMatrix',repmat([B1.',[1,A1.']],2,1));
[y1,t1] = stepz(BQ1);
[y2,t2] = stepz(BQ2);
xlim([0 3000])
The time to the first peak of the step response is the same for BQ1 and BQ2.
[gd1,wd1] = grpdelay(BQ1,1e4);
[gd2,wd2] = grpdelay(BQ2,1e4);
But the group delays are different. In fact they differ by a factor of 2, which is expected, because phase(BQ2) = 2*phase(BQ1).

Sign in to comment.

Accepted Answer

Paul on 13 Jan 2022
I thought group delay shifts the low frequency envelope of the signal and phase delay (phase / frequency) shifts the high frequency carrier of the signal. That is
u[n] = s[n]*cos(w0*n) -> y[n] = s[n - grpdelay]*cos(w0*(n - phasedelay)
But the example code doesn't include the low frequency envelope on the input sine wave. Because the phase of the filter is nearly exactly zero at 40 Hz and pretty close to zero at 200 Hz, the phase delay at those frequencies is also essentially zero, so we basically see zero phase shift in both outputs relative to the inputs.
Paul on 15 Jan 2022
Example from: Group delay and phase delay example. The results here differ slightly from the linked page.
Define a simple, second order filter with resonant frequency at f0
f0 = 10e3; % Hz
Fs = 1e6; % Hz, sampling frequency
[num,den] = bilinear(2*5000*[1 0],[1 2*5000 (2*pi*f0)^2],Fs);
We will test the filter at f0 and f1
f1 = 15e3; % Hz
Bode plot of the filter.
freqhz = logspace(3,5,1000);
hresp = freqz(num,den,freqhz,Fs);
semilogx(freqhz,angle(hresp));ylabel('Phase (rad)');
xlabel('Frequency (Hz)')
hold on
At f0, the phase of the filter is 0, so there should be zero phase delay for an input at f0, but the slope of the phase is high, so the goup delay will be large. At f1, the phase delay should be large, but the group delay is small because the phase curve is nearly flat.
Compute the phase and group delays.
pdelay = phasedelay(num,den,2*pi*[f0 f1]/Fs);
gdelay = grpdelay(num,den,2*pi*[f0 f1]/Fs);
Convert phase and group delay from samples to milliseconds
pdelay = pdelay(:).'/Fs*1000
pdelay = 1×2
0.0001 0.0147
gdelay = gdelay(:).'/Fs*1000
gdelay = 1×2
0.2001 0.0051
Define two signals at f0 at f1 enveloped by a slowly changing signal
Tf = 10e-3;
tvec = 0:1/Fs:Tf;
envelope = @(t) exp(-5e5*(t-Tf/2).^2);
sig1 = envelope(tvec) .* sin(2*pi*f0*tvec);
sig2 = envelope(tvec) .* sin(2*pi*f1*tvec);
Plot sig1 and its enevelope
Run both signals through the filter
y1 = filter(num,den,sig1);
y2 = filter(num,den,sig2);
Compare the first signal to the input
plot(tvec,sig1,tvec,envelope(tvec),tvec,y1,tvec,envelope(tvec - gdelay(1)/1000));
It looks like the envelope is delayed.
Plot again, zoom in around the peak and draw two vertical lines separated by the group delay.
plot(tvec,sig1,tvec,envelope(tvec),tvec,y1,tvec,envelope(tvec - gdelay(1)/1000));
xlim([.0044 .0056])
ylim([-1 1.1])
xline(.005 + gdelay(1)/1000,'g','LineWidth',2)
It looks like the peaks of the input and output envelopes are separated by the group delay. Note also that the carrier signal is not shifted from input to output, as expected because the phase delay is zero.
Repeat for the second signal, scale the output by the filter gain so we can see what's going on
mag = abs(interp1(freqhz,hresp,f1));
plot(tvec,sig2,tvec,envelope(tvec),tvec,y2/mag,tvec,envelope(tvec - gdelay(2)/1000));
It doesn't look like the envelope has shifted, as expected because gdelay(2) is small.
Plot again, zoom in and draw two vertical lines separated by the phase delay
plot(tvec,sig2,tvec,envelope(tvec),tvec,y2/mag,tvec,envelope(tvec - gdelay(2)/1000));
xlim([.0048 .0052])
ylim([-1 1.1])
xline(.00495 + pdelay(2)/1000,'g','LineWidth',2)
As expected, the carrier signal is delayed by the phase delay.

Sign in to comment.

More Answers (0)




Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!