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

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
tiledlayout(4,1)
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?

[William Rose]

@Richard Keniuk,

I'm not sure yet why the "output" signals in plot 4 look so similar to the "input" signals in plot 3. I initially misinterpreted the magnitude spectrum plot (plot 1), because it does not have the label dB on the vertical axis. I do not have access to the Audio Toolbox, so I cannot run your code.

[Richard Keniuk]

Thanks Paul,

That did the trick... I see now the application to group delay and how it damages high order modulation schemes in RF applications through envelope delay. My original demo of these two low frequency signals is not good because there is no envelope to show. The envelope delay is on the order of the sinewave risetime of the first 1/4 of the 40Hz signal so it's kinda hidden. I made a bursted 40Hz signal and ran it through my filter. The group delay shows up and explains the slow amplitude "build" in my 40Hz burst below. Phase timing is identical to the input.

Wouldn't the step response show the envelope shape due to the group delay? The peak is 244 samples or the group delay at 40Hz.

[Paul]

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);
figure
plot(t1,y1,t2,y2),grid
xlim([0 3000])
xline(244)

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);
figure
plot(wd1(1:100)/2/pi*Fs,gd1(1:100),wd2(1:100)/2/pi*Fs,gd2(1:100)),grid 
xline(40)

But the group delays are different. In fact they differ by a factor of 2, which is expected, because phase(BQ2) = 2*phase(BQ1).

Answer 1

[Paul]

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.

[Richard Keniuk]

You are applying group delay to the magnitude and I don't think that's the case. It is a calculation based on the phase response. "Group Delay is defined as the negative derivative (or slope) of the phase response vs. frequency." Courtesy of cjs-labs.com. I was in the understanding that the group delay of 244 samples plus the negative delay in the 200Hz signal would shift (delay) the 40Hz signal by 245 samples relative to the 200Hz signal. Since the signals in the plots are absolute I would think you'd see the delays but you can not (or only slightly). Maybe I've been misinterpreting "group delay". Group delay in FIR filters appears as actual delay by the number of samples.

[Paul]

I should have been clearer that s[n] is low frequency, i.e., low frequency relative to w0.

The steady state output of a stable filter in response to sinusoid input is a sinusoid of the same frequency, with ampliude scaled by the amplitude of the filter and phase shifted by the phase of the filter, both evaluated at the frequency of the input. I believe the results above show exactly this.

It is true that the group delay for a linear phase FIR filter is equal to the phase delay, so if the sinusoidal input signal (or sum of sinusoids I suppose) is in the pass band of the filter, the output will be the input delayed by the group delay, because that's the same as the phase delay for a linear phase filter. But BQ1 isn't a linear phase filter.

The Wikipedia link shows that if the input is the product of low freuqency signal (envelope) and a high frequency sinusoid (carrier), the group delay shifts the envelope and the phase delay shifts the carrier.

[Richard Keniuk]

Paul, I see now. Plotting the phase of my filter at 40 Hz correlates to roughly 0 phase shift which is the plot I've shown. This paper https://www.radio-labs.com/DesignFile/DN004.pdf explaines group delay as you say and shows some pictures. I will need to change my test signal and I will see what you are explaining. Thanks. My results threw me off due to my "distorted" intuition. ;-)

[Paul]

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.

figure
freqhz = logspace(3,5,1000);
hresp = freqz(num,den,freqhz,Fs);
subplot(211);
semilogx(freqhz,abs(hresp));ylabel('Magnitude')
subplot(212);
semilogx(freqhz,angle(hresp));ylabel('Phase (rad)');
xlabel('Frequency (Hz)')
hold on
xline(f0)
xline(f1)

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

figure
plot(tvec,sig1,tvec,envelope(tvec));
legend('sig1','envelope');

Run both signals through the filter

y1 = filter(num,den,sig1);
y2 = filter(num,den,sig2);

Compare the first signal to the input

figure
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.

figure
plot(tvec,sig1,tvec,envelope(tvec),tvec,y1,tvec,envelope(tvec - gdelay(1)/1000));
xlim([.0044 .0056])
ylim([-1 1.1])
xline(.005,'g','Linewidth',2)
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));
figure
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

figure
plot(tvec,sig2,tvec,envelope(tvec),tvec,y2/mag,tvec,envelope(tvec - gdelay(2)/1000));
grid
xlim([.0048 .0052])
ylim([-1 1.1])
xline(.00495,'g','Linewidth',2)
xline(.00495 + pdelay(2)/1000,'g','LineWidth',2)

As expected, the carrier signal is delayed by the phase delay.