帕克斯-麦克莱伦算法 (英语: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
,
f
o
r
n
<
0
a
n
d
n
≥
N
,
N
i
s
a
f
i
n
i
t
e
n
u
m
b
e
r
{\displaystyle h[n]=0,\;for\;n<0\;and\;n\geq N,\;N\,is\,a\,finite\;number}
有限脉冲响应滤波器的优点在于脉冲响应是有限的,使得设计上较为简单。然而如何在有限的点数下,设计出效果最近似于理想目标的滤波器,则是帕克斯-麦克莱伦算法所欲解决的问题。
对于滤波器设计,帕克斯-麦克莱伦算法的精神在于最小化最大误差。在忽略通带与止带之间转换带(transition band)的情况下,最小化通带与止带的最大误差:
Max
f
|
H
(
f
)
−
H
d
(
f
)
|
{\displaystyle {\underset {f}{\operatorname {Max} }}\left|H(f)-H_{d}(f)\right|}
其中
H
(
f
)
=
∑
n
=
−
∞
∞
h
[
n
]
e
−
j
2
π
f
n
{\displaystyle H(f)=\sum _{n=-\infty }^{\infty }h[n]e^{-j2\pi fn}}
为设计滤波器的频率响应,
H
d
(
f
)
{\displaystyle H_{d}(f)}
则为理想目标滤波器的频率响应。
在数码滤波器设计中,常常会将信号的频率做取样,使得频谱具有周期性。设计者即可针对一个周期去做计算就好,减少计算量。所以前两行的最大误差可写成:
Max
F
|
H
(
F
)
−
H
d
(
F
)
|
{\displaystyle {\underset {F}{\operatorname {Max} }}\left|H(F)-H_{d}(F)\right|}
其中
F
{\displaystyle F}
为正规化频率(normalized frequency):
F
=
f
f
s
{\displaystyle F={\frac {f}{f_{s}}}}
滤波器设计时,可利用weighting function将较重要的频带比重放大。如此一来,在利用帕克斯-麦克莱伦算法设计滤波器时,则会较重视比重较大频带的误差。
若在加入weighting function情况下,可将帕克斯-麦克莱伦算法一般化。此时的最大误差则可表示为:
Max
f
|
W
(
f
)
[
H
(
f
)
−
H
d
(
f
)
]
|
{\displaystyle {\underset {f}{\operatorname {Max} }}\left|W(f)\left[H(f)-H_{d}(f)\right]\right|}
另外在数学上,此种将向量取绝对值并找出某个最大的元素的算法,称为取
L
∞
{\displaystyle L_{\infty }}
范数。若能将
H
(
F
)
−
H
d
(
F
)
{\displaystyle H(F)-H_{d}(F)}
离散化写成矩阵的形式,就可以用此方法快速找出最大误差。
下面的文章将说明如何以该算法设计优化滤波器,假设
滤波器长度为N,且N为奇数可表示成
N
=
2
k
+
1
{\displaystyle N=2k+1}
目标滤波器的频率响应
H
d
(
F
)
{\displaystyle H_{d}(F)}
为偶函数
W
(
F
)
{\displaystyle W(F)}
用以表示所指定的权重函数(weighting function)。功用是将特定频段(通常是通带内)的误差调得更小,重视某频段的优化。
此算法共分为6个步骤:
设置极值点起始值
在范围
0
≤
F
≤
0.5
{\displaystyle 0\leq F\leq 0.5}
的范围内,任意选择
k
+
2
{\displaystyle k+2}
点频率
F
0
,
F
1
,
F
2
,
…
,
F
k
+
1
{\displaystyle F_{0},\;F_{1},\;F_{2},\;\ldots ,\;F_{k+1}}
作为极值点(extreme frequency)的起始值。
将此时的最大误差
E
1
{\displaystyle E_{1}}
设为
∞
{\displaystyle \infty }
,但所选择的
k
+
2
{\displaystyle k+2}
点起始值不能落在转换频带(transition band),也不能将所有的起始值设在止带(stop band)上。
极端频率为最后完成设计的滤波器频率响应中,会出现最大误差的频率。一开始所给定的起始值是随机的,会在此算法之后的步骤中逐渐收敛。
此时,令在各点极端频率的误差为
[
H
(
F
m
)
−
H
d
(
F
m
)
]
W
(
F
m
)
=
(
−
1
)
m
+
1
e
,
w
h
e
r
e
m
=
0
,
1
,
2
…
,
k
+
1
{\displaystyle \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
{\displaystyle r[n]=h[n+k],\;k=(N-1)/2}
。此
r
[
n
]
{\displaystyle r[n]}
是设计的滤波器响应
h
[
n
]
{\displaystyle h[n]}
的平移。
r
[
0
]
{\displaystyle r[0]}
为
h
[
n
]
{\displaystyle h[n]}
的正中央项 ( 举例:
r
[
0
]
=
h
[
k
]
,
r
[
1
]
=
h
[
k
+
1
]
{\displaystyle r[0]=h[k],r[1]=h[k+1]}
)。
s
[
0
]
=
r
[
0
]
,
s
[
n
]
=
2
r
[
n
]
,
f
o
r
0
<
n
≤
k
{\displaystyle s[0]=r[0],\ \ s[n]=2r[n],\;for\;0<n\leq k}
。因为
H
d
(
F
)
{\displaystyle H_{d}(F)}
为偶函数 ,所以
r
[
n
]
{\displaystyle r[n]}
也是偶函数 ,则再设计
s
[
n
]
{\displaystyle s[n]}
,计算
r
[
n
]
{\displaystyle r[n]}
的一半范围就好。
如此一来,可将在第1步骤中所得到的误差式表示为:
[
R
(
F
m
)
−
H
d
(
F
m
)
]
W
(
F
m
)
=
(
−
1
)
m
+
1
e
,
w
h
e
r
e
m
=
0
,
1
,
2
…
,
k
+
1
{\displaystyle \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
)
=
∑
n
=
0
k
s
[
n
]
cos
(
2
π
n
F
)
=
e
j
2
π
F
k
H
(
F
)
{\displaystyle R(F)=\sum _{n=0}^{k}s[n]\cos(2\pi nF)=e^{j2\pi Fk}H(F)}
(由于使用
s
[
n
]
{\displaystyle s[n]}
,计算项次从0到k)
经过整理之后可得
∑
n
=
0
k
s
[
n
]
cos
(
2
π
F
m
n
)
+
(
−
1
)
m
W
−
1
(
F
m
)
e
=
H
d
(
F
m
)
{\displaystyle \sum _{n=0}^{k}s[n]\cos(2\pi F_{m}n)+(-1)^{m}W^{-1}(F_{m})e=H_{d}(F_{m})}
上述的误差关系式,可表示为矩阵形式
A
x
=
H
{\displaystyle Ax=H}
[
1
cos
(
2
π
F
0
)
cos
(
4
π
F
0
)
⋯
cos
(
2
π
k
F
0
)
1
/
W
(
F
0
)
1
cos
(
2
π
F
1
)
cos
(
4
π
F
1
)
⋯
cos
(
2
π
k
F
1
)
−
1
/
W
(
F
1
)
⋮
⋮
⋮
⋱
⋮
⋮
1
cos
(
2
π
F
k
)
cos
(
4
π
F
k
)
⋯
cos
(
2
π
k
F
k
)
(
−
1
)
k
/
W
(
F
k
)
1
cos
(
2
π
F
k
+
1
)
cos
(
4
π
F
k
+
1
)
⋯
cos
(
2
π
k
F
k
+
1
)
(
−
1
)
k
+
1
/
W
(
F
k
+
1
)
]
[
s
[
0
]
s
[
1
]
⋮
s
[
k
]
e
]
=
[
H
d
[
F
0
]
H
d
[
F
1
]
⋮
H
d
[
F
k
]
H
d
[
F
k
+
1
]
]
{\displaystyle {\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
{\displaystyle x=A^{-1}H}
计算
s
[
0
]
,
s
[
1
]
,
…
,
s
[
k
]
{\displaystyle s[0],s[1],\ldots ,s[k]}
计算误差函数
计算误差函数
e
r
r
(
F
)
,
f
o
r
0
≤
F
≤
0.5
,
F
∉
t
r
a
n
s
i
t
i
o
n
b
a
n
d
{\displaystyle err(F),\;for0\leq F\leq 0.5\;,F\notin transition\;band}
e
r
r
(
F
)
=
[
R
(
F
)
−
H
d
(
F
)
]
W
(
F
)
=
{
∑
n
=
0
k
s
[
n
]
cos
(
2
π
F
n
)
−
H
d
(
F
)
}
W
(
F
)
{\displaystyle err(F)=\left[R(F)-H_{d}(F)\right]W(F)=\left\{\sum _{n=0}^{k}s[n]\cos(2\pi Fn)-H_{d}(F)\right\}W(F)}
查找极值点
从
e
r
r
(
F
)
{\displaystyle err(F)}
中,找k+2个区域极大或极小值,将区域极大或极小值出现的频率标示为
P
0
,
P
1
,
…
,
P
k
,
P
k
+
1
{\displaystyle P_{0},P_{1},\ldots ,P_{k},P_{k+1}}
区域极大或极小值的判断规则:
不是在边界处的区域高峰(peaks)或低谷(dips)。在此,边界区域即为
F
=
0
,
F
=
0.5
{\displaystyle F=0,F=0.5}
以及频率转换带的两边。
对于在边界区域的频率点,可用下列的标准判断是否为区域极大或极小:
S
i
g
n
(
e
f
f
(
F
d
)
)
{\displaystyle Sign(eff(F_{d}))}
及
S
i
g
n
(
e
r
r
′
(
F
d
)
)
{\displaystyle Sign(err'(F_{d}))}
为同相时,右边界是极值点;反相时,左边界是极值点;其他情况非极值点。
若找到多于
k
+
2
{\displaystyle k+2}
个极值点,选择极值点的优先顺位为:
优先选择不在边界的极值点。
其次选在边界的极值点中,
|
e
r
r
(
F
)
|
{\displaystyle \left|err(F)\right|}
较大的,直到凑足
k
+
2
{\displaystyle k+2}
个极值点为止。
当边界的极值点的
|
e
r
r
(
F
)
|
{\displaystyle \left|err(F)\right|}
一样大时,优先选择转换频带两旁的点。
判断误差是否收敛
计算误差的最大值,令其为
E
0
=
M
a
x
(
|
e
r
r
(
F
)
|
)
{\displaystyle E_{0}=Max(\left|err(F)\right|)}
。
若
E
0
{\displaystyle E_{0}}
为现在的误差最大值,
E
1
{\displaystyle E_{1}}
为前一轮计算的误差最大值,则利用下列规则判断算法的下一步:
若
E
1
−
E
0
>
Δ
,
o
r
E
1
−
E
0
<
0
{\displaystyle E_{1}-E_{0}>\Delta ,\;or\;E_{1}-E_{0}<0}
,设置
F
n
=
P
n
a
n
d
E
1
=
E
0
{\displaystyle F_{n}=P_{n}\;and\;E_{1}=E_{0}}
,回到步骤2。
若
0
≤
E
1
−
E
0
≤
Δ
{\displaystyle 0\leq E_{1}-E_{0}\leq \Delta }
,进行步骤6。
计算所设计滤波器的脉冲响应
h
[
n
]
{\displaystyle \;h[n]}
由先前在步骤2中的关系式,我们可得
h
[
k
]
=
s
[
0
]
,
h
[
k
+
n
]
=
s
[
n
]
/
2
,
h
[
k
−
n
]
=
s
[
n
]
/
2
,
f
o
r
n
=
1
,
2
,
3
,
…
,
k
{\displaystyle 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
]
{\displaystyle h[k]}
即为我们所求的脉冲响应。
权重函数
W
(
F
)
{\displaystyle W(F)}
滤波器响应的影响
编辑
当权重函数在带内设计为1,在带外设计为小于1,会让滤波器较重视通过带通频段的信号。
当权重函数在带外设计为1,在内设计为小于1,会让滤波器具有较好的滤除噪声的能力。
Jian-Jiun Ding (2013), Advanced Digital Signal Processing (页面存档备份 ,存于互联网档案馆 ) [viewed 27/06/2013]
T. W. Parks and J. H. McClellan, “Chebyshev Approximation for Nonrecursive Digital Filter with Linear Phase”, IEEE Trans. Circuit Theory, vol. 19, no. 2, pp. 189-194, March 1972.
J. H. McClellan, T. W. Parks, and L. R. Rabiner “A computer program for designing optimum FIR linear phase digital filter”, IEEE Trans. Audio- Electroacoustics, vol. 21, no. 6, Dec. 1973.
F. Mintz and B. Liu, “Practical design rules for optimum FIR bandpass digital filter”, IEEE Trans. ASSP, vol. 27, no. 2, Apr. 1979.