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

3 ビュー (過去 30 日間)
Richard Keniuk
Richard Keniuk 2021 年 12 月 14 日
編集済み: 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 件のコメント
Richard Keniuk
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
Paul 2022 年 1 月 15 日
編集済み: Paul 2022 年 1 月 15 日
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).

サインインしてコメントする。

採用された回答

Paul
Paul 2022 年 1 月 13 日
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 件のコメント
Richard Keniuk
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
Paul 2022 年 1 月 15 日
編集済み: Paul 2022 年 1 月 16 日
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.

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeDigital Filtering についてさらに検索

製品


リリース

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by