帕克斯-麦克莱伦算法

帕克斯-麦克莱伦算法（Parks–McClellan algorithm,又称为Remez-exchange algorithm、Mini-max algorithm)，为一个用以设计最佳化有限脉冲响应滤波器（finite impulse response filter）的迭代算法，由James McClellan和Thomas Parks于1972年的著作中提出。

此算法的主要精神，在于利用迭代的方式最小化滤波器在通带（pass band）和止带（stop band）的最大误差，因此有时也称为最小化最大误差算法（Mini-max filter design）。由于帕克斯-麦克莱伦算法也属于Remez-exchange algorithm为了设计有限脉冲响应滤波器而产生的一种变形，因此也有人以Remez-exchange algorithm代称。

## 有限脉冲响应滤波器设计

有限脉冲响应滤波器（finite impulse response filter）利用有限的点数来表示滤波器的脉冲响应，对于N点有限脉冲响应滤波器

$$h[n] = 0,\; for\; n<0\; and\; n\ge N, \;N\, is\, a\, finite\; number$$

有限脉冲响应滤波器的优点在于脉冲响应是有限的，使得设计上较为简单。然而如何在有限的点数下，设计出效果最近似于理想目标的滤波器，则是帕克斯-麦克莱伦算法所欲解决的问题。

对于滤波器设计，帕克斯-麦克莱伦算法的精神在于最小化最大误差。在忽略通带与止带之间转换带(transition band)的情况下，最小化通带与止带的最大误差：$\underset{f}{\operatorname{Max}} \left| H(f) - H_d(f) \right|$

其中$H(f) = \sum_{n=-\infty}^{\infty } h[n] e^{-j2\pi fn}$为设计滤波器的频率响应，$H_d(f)$则为理想目标滤波器的频率响应。

在数位滤波器设计中，常常会将信号的频率做取样，使得频谱具有周期性。设计者即可针对一个周期去做计算就好,减少计算量。所以前两行的最大误差可写成：

$$\underset{F}{\operatorname{Max}} \left| H(F) - H_d(F) \right|$$

其中$F$为正规化频率(normalized frequency)：$F=\frac{f}{f_s}$

滤波器设计时，可利用weighting function将较重要的频带比重放大。如此一来，在利用帕克斯-麦克莱伦算法设计滤波器时，则会较重视比重较大频带的误差。

若在加入weighting function情况下，可将帕克斯-麦克莱伦算法一般化。此时的最大误差则可表示为：$\underset{f}{\operatorname{Max}} \left| W(f) \left[H(f) - H_d(f) \right] \right|$

另外在数学上，此种将向量取绝对值并找出某个最大的元素的算法，称为取$L_{\infty}$范数。若能将$H(F) - H_d(F)$离散化写成矩阵的形式，就可以用此方法快速找出最大误差。

## 帕克斯-麦克莱伦算法

下面的文章将说明如何以该算法设计最佳化滤波器，假设

- 滤波器长度为N，且N为奇数可表示成$N = 2k+1$
- 目标滤波器的频率响应$H_d(F)$为偶函数
- $W(F)$用以表示所指定的权重函数(weighting function)。功用是将特定频段(通常是通带内)的误差调得更小，重视某频段的最佳化。

此算法共分为6个步骤：

- 设定极值点起始值
- 在范围$0\le F\le 0.5$的范围内，任意选择$k+2$点频率$F_0 ,\; F_1,\; F_2,\; \ldots,\; F_{k+1}$作为极值点(extreme frequency)的起始值。
- 将此时的最大误差$E_1$设为$\infty$，但所选择的$k+2$点起始值不能落在转换频带(transition band)，也不能将所有的起始值设在止带(stop band)上。
- 极端频率为最后完成设计的滤波器频率响应中，会出现最大误差的频率。一开始所给定的起始值是随机的，会在此算法之后的步骤中逐渐收敛。
- 此时，令在各点极端频率的误差为$\left[ H(F_m)-H_d(F_m) \right] W(F_m) = (-1)^{m+1} e,\; where\;m=0,1,2 \ldots , k+1$。其中e为设计滤波器响应式与理想滤波器响应式在相对应频率点的误差值。
- 计算目前的频率响应
- 为了方便算法运算之后的进行，我们可稍微整理误差的表示方式。若令
- $r[n]=h[n+k],\; k=(N-1)/2$。此$r[n]$是设计的滤波器响应$h[n]$的平移。$r[0]$为$h[n]$的正中央项 ( 举例：$r[0]=h[k], r[1]=h[k+1]$ )。
- $s[0]=r[0],\ \ s[n]=2r[n],\; for\; 0 < n \le k$。因为$H_d(F)$为偶函数，所以$r[n]$也是偶函数，则再设计$s[n]$，计算$r[n]$的一半范围就好。
- 如此一来，可将在第1步骤中所得到的误差式表示为：
- $$\left[ R(F_m)-H_d(F_m) \right] W(F_m) = (-1)^{m+1} e,\; where\;m=0,1,2 \ldots , k+1$$
- 其中。
- $R(F) = \sum_{n=0}^k s[n] \cos (2 \pi nF) = e^{j2 \pi Fk} H(F)$ (由于使用$s[n]$,计算项次从0到k)
- 经过整理之后可得
- $$\sum_{n=0}^k s[n] \cos (2 \pi F_m n) + (-1)^m W^{-1} (F_m) e = H_d(F_m)$$
- 上述的误差关系式，可表示为矩阵形式$Ax = H$
- $$\begin{bmatrix} 1 & \cos (2 \pi F_0) &\cos (4 \pi F_0) & \cdots & \cos (2 \pi kF_0)& 1/W(F_0) \\ 1 & \cos (2 \pi F_1) &\cos (4 \pi F_1) & \cdots & \cos (2 \pi kF_1)& -1/W(F_1) \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ 1 & \cos (2 \pi F_k) &\cos (4 \pi F_k) & \cdots & \cos (2 \pi kF_k)& (-1)^k/W(F_k) \\ 1 & \cos (2 \pi F_{k+1}) &\cos (4 \pi F_{k+1}) & \cdots & \cos (2 \pi kF_{k+1})& (-1)^{k+1}/W(F_{k+1}) \\ \end{bmatrix} \begin{bmatrix} s[0] \\ s[1] \\ \vdots \\ s[k] \\e \end{bmatrix} = \begin{bmatrix} H_d[F_0] \\ H_d[F_1] \\ \vdots \\ H_d[F_k] \\ H_d[F_{k+1}] \end{bmatrix}$$
- 因此，我们可由$x = A^{-1} H$计算$s[0], s[1], \ldots, s[k]$
- 计算误差函数
- 计算误差函数$err(F), \; for 0 \le F \le 0.5 \; , F \notin transition \; band$
- $$err(F)= \left[ R(F)-H_d(F) \right] W(F) = \left\{ \sum_{n=0}^k s[n] \cos (2 \pi F n) - H_d(F) \right\} W(F)$$
- 寻找极值点
- 从$err(F)$中，找k+2个区域极大或极小值，将区域极大或极小值出现的频率标示为$P_0, P_1, \ldots, P_k, P_{k+1}$
- 不是在边界处的区域高峰(peaks)或低谷(dips)。在此，边界区域即为$F=0, F=0.5$以及频率转换带的两边。
- 对于在边界区域的频率点，可用下列的标准判断是否为区域极大或极小：$Sign(eff(F_d))$及$Sign(err'(F_d))$为同相时，右边界是极值点；反相时，左边界是极值点；其他情况非极值点。
- 若找到多于$k+2$个极值点，选择极值点的优先顺位为：
- 优先选择不在边界的极值点。
- 其次选在边界的极值点中，$\left| err(F) \right|$较大的，直到凑足$k+2$个极值点为止。
- 当边界的极值点的$\left| err(F) \right|$一样大时，优先选择转换频带两旁的点。
- 判断误差是否收敛
- 计算误差的最大值，令其为$E_0 = Max(\left| err(F) \right|)$。
- 若$E_0$为现在的误差最大值，$E_1$为前一轮计算的误差最大值，则利用下列规则判断算法的下一步：
- 若$E_1 - E_0 > \Delta, \; or \; E_1 - E_0 < 0$，设定$F_n = P_n \; and \; E_1 = E_0$，回到步骤2。
- 若$0 \le E_1 - E_0 \le \Delta$，进行步骤6。
- 计算所设计滤波器的脉冲响应$\;h[n]$
- 由先前在步骤2中的关系式，我们可得
- $$h[k] = s[0], h[k+n] = s[n]/2, h[k-n]=s[n]/2, \; for \; n=1,2,3, \ldots, k$$
- 此$h[k]$即为我们所求的脉冲响应。

### 权重函数$W(F)$滤波器响应的影响

当权重函数在带内设计为1，在带外设计为小于1，会让滤波器较重视通过带通频段的讯号。

当权重函数在带外设计为1，在内设计为小于1，会让滤波器具有较好的滤除噪声的能力。

## 特征

用此方法设计出来的滤波器，一定会满足以下两个情况：

- 有k+2个以上的极值点(极大点与极小点)。
- 在极点上，$W(F_m)|R(F_m)-H_d(F_m)|$

### 过渡频段

过渡频段(transition band)对设计滤波器的误差也会有影响。将过渡频段设计地窄一些，$R(F)$和$H_d$的误差就会比较大；将过渡频段设计地宽一些，$R(F)$和$H_d$的误差就会比较小。再设计上可以适当的增加过渡频段宽度，让通带和止带地响应更趋近于理想值。

假设我们想要带通段的涟波小于等于$\delta_1$，带止段的涟波小于等于$\delta_2$,过渡频段小于$\Delta F$ ($\Delta F=\frac{|f_1-f_2|}{f_s}$ $f_1$，$f_2$为过渡频段的上、下界)。则要设计滤波器长度$N$为：

$$N\geq\frac{2}{3\Delta F}\log_{10}(\frac{1}{\delta_1 \delta_2})$$

移项可得：$\delta_1 \delta_2=10^{\frac{-3N\Delta F}{2}-1}$

若要设计的频段中有多的过渡频段，则$\Delta F$取最小长度的过渡频段带入计算。

对于一固定长度的数位滤波器，再设计上可以牺牲一些频段留给过渡频段，将涟波缩小。但要注意不可将过渡频段设计过长，因为过渡频段是无法使用的。