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

Question

Richard Keniuk 2021-12-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1610470-trying-to-demonstrate-group-delay-but-not-seeing-any-in-dsp-biquadfilter

编辑： Paul 2022-1-16

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?

3 个评论
显示 1更早的评论隐藏 1更早的评论

Richard Keniuk 2022-1-15

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 2022-1-15

编辑：Paul 2022-1-15

在 MATLAB Online 中打开

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 2022-1-13

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1610470-trying-to-demonstrate-group-delay-but-not-seeing-any-in-dsp-biquadfilter#answer_874260

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.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Paul 2022-1-14

编辑：Paul 2022-1-14

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 2022-1-14

编辑：Richard Keniuk 2022-1-14

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 2022-1-15

编辑：Paul 2022-1-16

在 MATLAB Online 中打开

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.

请先登录，再进行评论。

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

3 个评论
显示 1更早的评论隐藏 1更早的评论

采纳的回答

4 个评论
显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

3 个评论 显示 1更早的评论隐藏 1更早的评论

采纳的回答

4 个评论 显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

3 个评论
显示 1更早的评论隐藏 1更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论