matlab群延时函数,群延迟函数(group delay function)&群延迟滤波器 | 学步园

最近看了许多介绍Group delay function的论文，文章中大篇幅提到Group delay，group delay of digital filters,对这个方面的知识好像还挺有用的，所以想把它记录下来。然后总结下计算Group delay function的步骤。

假设有N个样本的脉冲响应为h(n)的数字滤波器，n为时间序列标号，对该数字滤波器进行离散时间傅里叶变换(discrete-time Fourier transform (DTFT)),H(ω),用极点的方式表示如下：

      (1)

在公式(1), M(ω) 是滤波器的频幅响应, Φ(ω) 是滤波器的相频响应 ， ω 是连续的角频率，单位为弧度每秒rad/s ,用H(ω)对 ω求导, 得到如下的等式:

      (2)

用公式 (2) 除以M(ω), 写成如下的形式：

     (3)

对公式(3)的右边第一项进行求导展开，得到如下的形式：

   (4)

对等式(4)两边同时除以j/ejΦ(ω) 得到如下等式：

     (5)

从上面的等式(5)我们能得到什么呢？通过观察，我们发现如下熟悉的项：

jd[H(ω)]/dω = n·h(n) 的DTFT

M(ω)•ejΦ(ω) =H(ω) =h(n) DDTFT

–d[Φ(ω)]/dω = 群延迟滤波器(group delay of the filter )

现在可以把上面的公式改成更直观的形式：

     (6)

用DFT代替DTFT可以得到计算群延迟数字滤波器的公式：

   (7)

从上面的公式看出，计算群延迟滤波器的过程。公式(7)是传统群延迟h(n)]值为0的情况。

下面用Matlab演示公式 (7)的计算过程：

clear, clc

Npts = 128; % 画图点数

B = [0.03, 0.0605, 0.121, 0.0605, 0.03];

A = [1, -1.194, 0.436];

Imp_Resp_Length = 40;

[Imp_Resp,n] = impz(B,A,Imp_Resp_Length);

ImpResp_times_Time = Imp_Resp.*n;

[Freq_Resp, W] = freqz(Imp_Resp, 1, Npts, 'whole');

[Deriv_of_Freq_Resp, W] = freqz(ImpResp_times_Time, ...

1, Npts, 'whole');

Grp_Delay = real(Deriv_of_Freq_Resp./Freq_Resp);

%[gdm, fgm] = grpdelay(B,A,Npts,'whole');gdm=Grp_Delay,w=fgm(rad) % Compute Group Delay

Grp_Delay = fftshift(Grp_Delay);% 将序列中后半部分循环移动到前面

Freq = (W-pi)/(2*pi); % 频率轴

figure(1), clf

subplot(2,1,1)

plot(n, Imp_Resp, '-ks', ...

n, ImpResp_times_Time, '-bs', 'markersize', 4)

legend('h(n)','n times h(n)');

ylabel('Amplitude'), xlabel('n'), grid on, zoom on

subplot(2,1,2)

plot(Freq, Grp_Delay,'-rs', 'markersize', 4)

ylabel('Group Delay (samples)')

xlabel('Freq x Fs (Fs = sample rate)')

grid on, zoom on