Yulewalk滤波器设计：MA系数确定及numf双边用法解析

一、MA系数（分子B）的推导公式与确定逻辑 #

Yule-Walk算法通过最小二乘拟合目标幅频响应，MA系数的推导基于线性系统卷积关系与最小二乘求解：

核心推导公式 #

已知滤波器传递函数：
$$
\frac{B(z)}{A(z)} = H(z) \implies B(z) = H(z)A(z)
$$
转换为时域脉冲响应的卷积关系：
$$
b[n] = h[n] * a[n]
$$
其中：

• $b[n]$ 为MA系数序列（长度为 $n_b+1$）
• $h[n]$ 是滤波器 $H(z)=B(z)/A(z)$ 的脉冲响应
• $a[n]$ 是AR系数序列（分母，固定 $a[0]=1$）

\begin{bmatrix}
a[0] & 0 & \dots & 0 \
a[1] & a[0] & \dots & 0 \
\vdots & \vdots & \ddots & \vdots \
a[n_b] & a[n_{b}-1] & \dots & a[0]
\end{bmatrix}
\begin{bmatrix}
b[0] \ b[1] \ \vdots \ b[n_b]
\end{bmatrix}
$$
通过最小二乘求解该方程组，即可得到MA系数 $b[n]$。

二、numf函数使用双边输入的原因 #

代码中两次调用numf的输入分别为：

1. Qh = numf([R(1)/2, R(2:nr)], A, na)
2. B = real(numf(hh(1:nr), A, na))

双边输入的本质 #

这里的输入是自相关序列的双边形式（或脉冲响应的截断双边序列），原因在于：

• Yule-Walk算法中，幅频响应平方的IFFT得到的自相关序列 $R[k]$ 是双边的（包含正、负延迟分量），R(1)/2是对直流分量的修正（对称序列的0延迟项仅存在一次，需折半匹配双边序列的定义）
• 频域处理后通过IFFT得到的脉冲响应 $h[n]$ 是周期延拓的双边序列，截断前nr个点可保留完整的有效时域信息，满足线性方程组构建的采样需求

numf函数内部逻辑 #

numf的核心是构建Toeplitz矩阵并求解最小二乘：

function b = numf(h, a, nb)
nh = max(size(h)); 
impr = filter(1, a, [1 zeros(1, nh-1)]); % 生成A(z)的脉冲响应
b = h/toeplitz(impr, [1 zeros(1, nb)])'; % 最小二乘求解B系数

• impr是分母 $A(z)$ 的逆系统脉冲响应，用于构建符合卷积关系的Toeplitz矩阵
• 双边输入的h包含了足够的时域采样点，保证最小二乘解的准确性与稳定性

三、代码中的数值处理技巧详解 #

1. 幅频响应的插值与对称化 #

npt = 512 + 1;
Ht(1:end-1) = interp1(ff, aa, linspace(0, 1, npt-1), 'linear');
Ht = [Ht Ht(npt-1:-1:2)]; % 对称化生成全频域响应

• 用interp1将用户给定的离散频点插值到512个采样点，提升频域采样密度，保证后续IFFT的精度
• 对称化操作是为了生成实值信号的幅频响应（实信号FFT满足共轭对称性），确保后续IFFT得到实值自相关序列

2. 自相关序列的窗函数截断 #

R = real(ifft(Ht .* Ht)); % 幅频平方的IFFT得到自相关序列
R = R(1:nr) .* (0.54 + 0.46 * cos(pi * nt / (nr - 1))); % 汉宁窗加窗

• $|H(e^{j\omega})|2$ 的IFFT对应滤波器的自相关序列，取前nr=4*na个点是为了截断到有效自相关分量长度
• 汉宁窗用于抑制自相关序列的旁瓣，减少频域混叠，提升AR系数求解的稳定性

3. 对数幅度与相位恢复 #

Ss = 2 * real(freqz(Qh, A, n, 'whole'))';
hh = ifft(exp(fft(Rwindow .* ifft(log(Ss)))));

• Ss是AR分量重构的功率谱，取对数后做IFFT得到复倒谱，再通过Rwindow提取因果部分，最后指数变换恢复相位信息
• 这一步从幅频响应中恢复因果脉冲响应，保证求解的MA系数对应因果稳定的滤波器

4. 多项式稳定性修正 #

A = polystab(denf(R, na)); % 修正分母多项式，保证滤波器稳定

• denf通过Modified Yule-Walker方法求解AR系数，可能出现单位圆外的不稳定极点，polystab将这些极点反射到单位圆内，确保滤波器的稳定性

总结 #

MA系数基于时域卷积的最小二乘求解确定，numf的双边输入是为了匹配自相关序列与脉冲响应的双边特性，保证线性方程组的完整性。代码中的数值处理围绕频域采样密度、自相关截断、稳定性修正、相位恢复四个核心，最终实现目标幅频响应的最优拟合。

内容的提问来源于stack exchange，提问作者feilong du