1. 四种FIR滤波器对称结构
设FIR滤波器的输入和输出序列分别为
x
[
n
]
x[n]
x[n]和
y
[
n
]
y[n]
y[n],滤波器系数为
ω
m
,
0
≤
m
<
M
\omega_m,0\leq m<M
ωm,0≤m<M,则滤波过程可用如下差分方程表示
y
[
n
]
=
∑
m
=
0
M
−
1
ω
m
x
[
n
−
m
]
(1)
y[n]=\sum_{m=0}^{M-1}\omega_m x[n-m] \tag{1}
y[n]=m=0∑M−1ωmx[n−m](1)
此外,可将输出序列写成系统冲激响应
h
[
n
]
h[n]
h[n]与输入序列的卷积,如下
y
[
n
]
=
∑
m
=
0
M
−
1
h
[
m
]
x
[
n
−
m
]
(2)
y[n]=\sum_{m=0}^{M-1}h[m]x[n-m] \tag{2}
y[n]=m=0∑M−1h[m]x[n−m](2)
从上面的表达式可以看出,当以样点为单位进行延时时,FIR滤波器的抽头系数和其冲激响应本质上是一致的。
FIR滤波器最吸引人的性质是可以实现线性相位(保证信号过滤波器后不失真)。具有线性相位的FIR滤波器的冲激响应满足如下对称条件:
h
[
n
]
=
±
h
[
M
−
1
−
n
]
(4)
h[n]=\pm h[M-1-n] \tag{4}
h[n]=±h[M−1−n](4)
上面的表达式包含4种情况,分别称为1/2/3/4型FIR滤波器,对应关系如下:
{
1
型:
h
[
n
]
=
h
[
M
−
1
−
n
]
,
M
为奇数
2
型:
h
[
n
]
=
h
[
M
−
1
−
n
]
,
M
为偶数
3
型:
h
[
n
]
=
−
h
[
M
−
1
−
n
]
,
M
为奇数
3
型:
h
[
n
]
=
−
h
[
M
−
1
−
n
]
,
M
为偶数
\left\{\begin{matrix} 1 型:h[n]=h[M-1-n], M为奇数 \\ 2 型:h[n]=h[M-1-n], M为偶数 \\ 3 型:h[n]=-h[M-1-n], M为奇数 \\ 3 型:h[n]=-h[M-1-n], M为偶数 \end{matrix}\right.
⎩
⎨
⎧1型:h[n]=h[M−1−n],M为奇数2型:h[n]=h[M−1−n],M为偶数3型:h[n]=−h[M−1−n],M为奇数3型:h[n]=−h[M−1−n],M为偶数
下面以1型FIR滤波器为例,推导其频率响应的表达式,其他形式的FIR滤波器可按类似方式进行推导。
H
(
e
j
ω
)
=
∑
m
=
0
N
−
1
h
[
n
]
e
−
j
n
ω
=
(
∑
n
=
0
(
N
−
3
)
/
2
h
[
n
]
e
−
j
n
ω
)
+
h
[
N
−
1
2
]
e
−
j
(
N
−
1
)
ω
/
2
+
(
∑
n
=
(
N
+
1
)
/
2
N
−
1
h
[
n
]
e
−
j
n
ω
)
H(e^{j\omega})=\sum_{m=0}^{N-1}h[n]e^{-jn\omega}=(\sum_{n=0}^{(N-3)/2}h[n]e^{-jn\omega})+h[\frac{N-1}{2}]e^{-j(N-1)\omega /2}+(\sum_{n=(N+1)/2}^{N-1}h[n]e^{-jn\omega})
H(ejω)=m=0∑N−1h[n]e−jnω=(n=0∑(N−3)/2h[n]e−jnω)+h[2N−1]e−j(N−1)ω/2+(n=(N+1)/2∑N−1h[n]e−jnω)
由于
N
N
N为奇数,所以上式以中间点为对称重新,将求和操作分成前面部分+中间部分+后面部分。对后面部分,令
m
=
N
−
1
−
n
m=N-1-n
m=N−1−n,则
∑
n
=
(
N
+
1
)
/
2
N
−
1
h
[
n
]
e
−
j
n
ω
=
∑
m
=
0
(
N
−
3
)
/
2
h
[
N
−
1
−
m
]
e
−
j
(
N
−
1
−
m
)
ω
=
∑
m
=
0
(
N
−
3
)
/
2
h
[
m
]
e
−
j
(
N
−
1
−
m
)
ω
\sum_{n=(N+1)/2}^{N-1}h[n]e^{-jn\omega}=\sum_{m=0}^{(N-3)/2}h[N-1-m]e^{-j(N-1-m)\omega}=\sum_{m=0}^{(N-3)/2}h[m]e^{-j(N-1-m)\omega}
∑n=(N+1)/2N−1h[n]e−jnω=∑m=0(N−3)/2h[N−1−m]e−j(N−1−m)ω=∑m=0(N−3)/2h[m]e−j(N−1−m)ω,将其带入上式得到:
H
(
e
j
ω
)
=
(
∑
n
=
0
(
N
−
3
)
/
2
h
[
n
]
e
−
j
n
ω
)
+
h
[
N
−
1
2
]
e
−
j
(
N
−
1
)
ω
/
2
+
(
∑
n
=
0
(
N
−
3
)
/
2
h
[
n
]
e
−
j
(
N
−
1
−
n
)
ω
)
=
(
∑
n
=
0
(
N
−
3
)
/
2
h
[
n
]
(
e
−
j
n
ω
+
e
−
j
(
N
−
1
−
n
)
ω
)
)
+
h
[
N
−
1
2
]
e
−
j
(
N
−
1
)
ω
/
2
=
e
−
j
(
N
−
1
)
ω
/
2
[
∑
n
=
0
(
N
−
3
)
/
2
h
[
n
]
(
e
−
j
ω
[
n
−
(
N
−
1
)
/
2
]
+
e
j
ω
[
n
−
(
N
−
1
)
/
2
]
)
+
h
[
N
−
1
2
]
]
=
e
−
j
(
N
−
1
)
ω
/
2
[
2
∑
n
=
0
(
N
−
3
)
/
2
h
[
n
]
c
o
s
(
N
−
1
2
−
n
)
ω
+
h
[
N
−
1
2
]
]
H(e^{j\omega})=(\sum_{n=0}^{(N-3)/2}h[n]e^{-jn\omega})+h[\frac{N-1}{2}]e^{-j(N-1)\omega /2}+(\sum_{n=0}^{(N-3)/2}h[n]e^{-j(N-1-n)\omega})\\=(\sum_{n=0}^{(N-3)/2}h[n](e^{-jn\omega}+e^{-j(N-1-n)\omega}))+h[\frac{N-1}{2}]e^{-j(N-1)\omega /2}\\=e^{-j(N-1)\omega/2}[\sum_{n=0}^{(N-3)/2}h[n](e^{-j\omega[n-(N-1)/2]}+e^{j\omega[n-(N-1)/2]})+h[\frac{N-1}{2}]]\\=e^{-j(N-1)\omega/2}[2\sum_{n=0}^{(N-3)/2}h[n]cos(\frac{N-1}{2}-n)\omega+h[\frac{N-1}{2}]]
H(ejω)=(n=0∑(N−3)/2h[n]e−jnω)+h[2N−1]e−j(N−1)ω/2+(n=0∑(N−3)/2h[n]e−j(N−1−n)ω)=(n=0∑(N−3)/2h[n](e−jnω+e−j(N−1−n)ω))+h[2N−1]e−j(N−1)ω/2=e−j(N−1)ω/2[n=0∑(N−3)/2h[n](e−jω[n−(N−1)/2]+ejω[n−(N−1)/2])+h[2N−1]]=e−j(N−1)ω/2[2n=0∑(N−3)/2h[n]cos(2N−1−n)ω+h[2N−1]]
上面的表达式中,令
m
=
(
N
−
1
2
−
n
m=(\frac{N-1}{2}-n
m=(2N−1−n,则可进一步得到:
H
(
e
j
ω
)
=
e
−
j
(
N
−
1
)
ω
/
2
[
2
∑
n
=
1
(
N
−
1
)
/
2
h
[
N
−
1
2
−
n
]
c
o
s
(
n
ω
)
+
h
[
N
−
1
2
]
]
H(e^{j\omega})=e^{-j(N-1)\omega/2}[2\sum_{n=1}^{(N-1)/2}h[\frac{N-1}{2}-n]cos(n\omega)+h[\frac{N-1}{2}]]
H(ejω)=e−j(N−1)ω/2[2n=1∑(N−1)/2h[2N−1−n]cos(nω)+h[2N−1]]
令
a
[
n
]
=
{
h
[
N
−
1
2
]
,
n
=
0
2
h
[
N
−
1
2
−
n
]
,
0
<
n
≤
N
−
1
2
a[n]=\left\{\begin{matrix} h[\frac{N-1}{2}],n=0 \\ 2h[\frac{N-1}{2}-n], 0<n\leq \frac{N-1}{2} \end{matrix}\right.
a[n]={h[2N−1],n=02h[2N−1−n],0<n≤2N−1
则
H
(
e
j
ω
)
=
e
−
j
(
N
−
1
)
ω
/
2
∑
n
=
0
(
N
−
1
)
/
2
a
[
n
]
c
o
s
(
n
ω
)
H(e^{j\omega})=e^{-j(N-1)\omega/2}\sum_{n=0}^{(N-1)/2}a[n]cos(n\omega)
H(ejω)=e−j(N−1)ω/2∑n=0(N−1)/2a[n]cos(nω)。显然,此时
H
(
e
j
ω
)
H(e^{j\omega})
H(ejω)具有线性相位。令
ϕ
(
ω
)
=
a
r
g
(
H
(
e
j
ω
)
)
=
−
(
N
−
1
)
ω
/
2
H
r
(
e
j
ω
)
=
∑
n
=
0
(
N
−
1
)
/
2
a
[
n
]
c
o
s
(
n
ω
)
\phi(\omega)=arg(H(e^{j\omega}))=-(N-1)\omega/2\\ H_r(e^{j\omega})=\sum_{n=0}^{(N-1)/2}a[n]cos(n\omega)
ϕ(ω)=arg(H(ejω))=−(N−1)ω/2Hr(ejω)=n=0∑(N−1)/2a[n]cos(nω)
则
H
(
e
j
ω
)
=
e
ϕ
(
ω
)
H
r
(
e
j
ω
)
H(e^{j\omega})=e^{\phi(\omega)}H_r(e^{j\omega})
H(ejω)=eϕ(ω)Hr(ejω)。另外3种FIR滤波器的频响可以按类似的方式进行推导,直接将推导结果总结在下面:
滤波器类型 | ϕ ( ω ) \phi(\omega) ϕ(ω) | H r ( e j ω ) H_r(e^{j\omega}) Hr(ejω) |
---|---|---|
1型 | − N − 1 2 ω -\frac{N-1}{2}\omega −2N−1ω | 2 ∑ n = 1 ( N − 1 ) / 2 h [ N − 1 2 − n ] c o s ( n ω ) + h [ N − 1 2 ] = ∑ n = 0 ( N − 1 ) / 2 a [ n ] c o s ( n ω ) 2\sum_{n=1}^{(N-1)/2}h[\frac{N-1}{2}-n]cos(n\omega)+h[\frac{N-1}{2}]=\sum_{n=0}^{(N-1)/2}a[n]cos(n\omega) 2∑n=1(N−1)/2h[2N−1−n]cos(nω)+h[2N−1]=∑n=0(N−1)/2a[n]cos(nω) |
2型 | − N − 1 2 ω -\frac{N-1}{2}\omega −2N−1ω |
2
∑
n
=
0
(
N
−
1
)
/
2
h
[
n
]
c
o
s
ω
(
N
−
1
2
−
n
)
=
∑
n
=
1
N
/
2
b
[
n
]
c
o
s
ω
(
n
−
1
/
2
)
2\sum_{n=0}^{(N-1)/2}h[n]cos\omega(\frac{N-1}{2}-n)=\sum_{n=1}^{N/2}b[n]cos\omega(n-1/2)
2∑n=0(N−1)/2h[n]cosω(2N−1−n)=∑n=1N/2b[n]cosω(n−1/2) b [ n ] = 2 h [ N 2 − n ] b[n]=2h[\frac{N}{2}-n] b[n]=2h[2N−n] |
3型 | π 2 − N − 1 2 ω \frac{\pi}{2}-\frac{N-1}{2}\omega 2π−2N−1ω |
2
∑
n
=
0
(
N
−
3
)
/
2
h
[
n
]
s
i
n
ω
(
N
−
1
2
−
n
)
=
∑
n
=
1
(
N
−
1
)
/
2
c
[
n
]
s
i
n
ω
n
2\sum_{n=0}^{(N-3)/2}h[n]sin\omega(\frac{N-1}{2}-n)=\sum_{n=1}^{(N-1)/2}c[n]sin \omega n
2∑n=0(N−3)/2h[n]sinω(2N−1−n)=∑n=1(N−1)/2c[n]sinωn c [ n ] = 2 h [ N − 1 2 − n ] c[n]=2h[\frac{N-1}{2}-n] c[n]=2h[2N−1−n] |
4型 | π 2 − N − 1 2 ω \frac{\pi}{2}-\frac{N-1}{2}\omega 2π−2N−1ω |
2
∑
n
=
0
N
/
2
−
1
h
[
n
]
s
i
n
ω
(
N
−
1
2
−
n
)
=
∑
n
=
1
N
/
2
d
[
n
]
s
i
n
ω
(
n
−
1
/
2
)
2\sum_{n=0}^{N/2-1}h[n]sin\omega(\frac{N-1}{2}-n)=\sum_{n=1}^{N/2}d[n]sin \omega (n-1/2)
2∑n=0N/2−1h[n]sinω(2N−1−n)=∑n=1N/2d[n]sinω(n−1/2) d [ n ] = 2 h [ N 2 − n ] d[n]=2h[\frac{N}{2}-n] d[n]=2h[2N−n] |
从上表中可见,3/4型滤波器相比1/2型滤波器需要增加额外的
π
/
2
\pi/2
π/2的相移器,此外,对3/4型滤波器来说,
H
r
(
0
)
=
0
H_r(0)=0
Hr(0)=0,并不适用于低通滤波器的设计。一般情况下,我们设计时更多采用I型滤波器。
基于上表中的表达式,我们期望将
H
r
(
e
j
ω
)
H_r(e^{j\omega})
Hr(ejω)用统一的形式,即
H
r
(
ω
)
=
Q
(
ω
)
P
(
ω
)
H_r(\omega)=Q(\omega)P(\omega)
Hr(ω)=Q(ω)P(ω),下面直接给出推导结果
滤波器类型 | Q ( ω ) Q(\omega) Q(ω) | P ( ω ) P(\omega) P(ω) |
---|---|---|
1型 | 1 1 1 | ∑ n = 0 ( N − 1 ) / 2 a [ n ] c o s ( ω n ) \sum_{n=0}^{(N-1)/2}a[n]cos(\omega n) ∑n=0(N−1)/2a[n]cos(ωn) |
2型 | c o s ( ω 2 ) cos(\frac{\omega}{2}) cos(2ω) |
∑
n
=
0
N
/
2
−
1
b
‾
[
n
]
c
o
s
(
ω
n
)
\sum_{n=0}^{N/2-1}\overline{b}[n]cos(\omega n)
∑n=0N/2−1b[n]cos(ωn) b ‾ [ 0 ] = 1 2 b [ 1 ] \overline{b}[0]=\frac{1}{2}b[1] b[0]=21b[1] b ‾ [ n ] = 2 b [ n ] − b [ n − 1 ] , n = 1 , . . . N / 2 \overline{b}[n]=2b[n]-b[n-1],n=1,...N/2 b[n]=2b[n]−b[n−1],n=1,...N/2 b ‾ [ N 2 − 1 ] = 2 b [ N 2 ] \overline{b}[\frac{N}{2}-1]=2b[\frac{N}{2}] b[2N−1]=2b[2N] |
3型 | s i n ( ω ) sin(\omega) sin(ω) |
∑
n
=
0
(
N
−
3
)
/
2
c
‾
[
n
]
c
o
s
(
ω
n
)
\sum_{n=0}^{(N-3)/2}\overline{c}[n]cos(\omega n)
∑n=0(N−3)/2c[n]cos(ωn) c ‾ [ M − 3 2 ] = c [ M − 1 2 ] \overline{c}[\frac{M-3}{2}]=c[\frac{M-1}{2}] c[2M−3]=c[2M−1] c ‾ [ M − 5 2 ] = 2 c [ M − 3 2 ] \overline{c}[\frac{M-5}{2}]=2c[\frac{M-3}{2}] c[2M−5]=2c[2M−3] c ‾ [ n − 1 ] − c ‾ [ n + 1 ] = 2 c [ n ] , 2 ≤ n ≤ M − 5 2 \overline{c}[n-1]-\overline{c}[n+1]=2c[n],2\leq n\leq \frac{M-5}{2} c[n−1]−c[n+1]=2c[n],2≤n≤2M−5 c ‾ [ 0 ] − 1 2 c ‾ [ 2 ] = c [ 1 ] \overline{c}[0]-\frac{1}{2}\overline{c}[2]=c[1] c[0]−21c[2]=c[1] |
4型 | s i n ( ω 2 ) sin(\frac{\omega}{2}) sin(2ω) |
∑
n
=
0
N
/
2
−
1
d
‾
[
n
]
c
o
s
(
ω
n
)
\sum_{n=0}^{N/2-1}\overline{d}[n]cos(\omega n)
∑n=0N/2−1d[n]cos(ωn) d ‾ [ M 2 − 1 ] = 2 d [ M 2 ] \overline{d}[\frac{M}{2}-1]=2d[\frac{M}{2}] d[2M−1]=2d[2M] d ‾ [ n − 1 ] − d ‾ [ n ] = 2 d [ n ] , 2 ≤ n ≤ M 2 − 1 \overline{d}[n-1]-\overline{d}[n]=2d[n],2\leq n\leq \frac{M}{2}-1 d[n−1]−d[n]=2d[n],2≤n≤2M−1 d ‾ [ 0 ] − 1 2 d ‾ [ 1 ] = d [ 1 ] \overline{d}[0]-\frac{1}{2}\overline{d}[1]=d[1] d[0]−21d[1]=d[1] |
2. 最优滤波器设计
在上一节中,我们知道四类FIR滤波器的频响实部可统一表示成如下形式:
H
r
(
ω
)
=
Q
(
ω
)
P
(
ω
)
=
Q
(
ω
)
(
∑
n
=
0
L
α
(
n
)
c
o
s
(
ω
n
)
)
H_r(\omega)=Q(\omega)P(\omega)=Q(\omega)(\sum_{n=0}^{L}\alpha(n)cos(\omega n))
Hr(ω)=Q(ω)P(ω)=Q(ω)(n=0∑Lα(n)cos(ωn))
基于上述通用表达式,我们进行最优FIR滤波器设计。最优FIR滤波器设计的基本思路是在某种优化目标下,让设计的滤波器频响和理论频响的误差最小。由于上面4中类型的FIR滤波器的相频特性与滤波器系数无关,故优化过程无需考虑相频的作用。设
H
d
r
(
ω
)
H_{dr}(\omega)
Hdr(ω)为期望的频响的实部,通常地,期望频响在通带内的增益为1,在阻带内的增益为0。此外,定义一个加权误差函数
W
(
ω
)
W(\omega)
W(ω)用于控制通带和阻带的波动。具体地
W
(
ω
)
=
{
δ
2
/
δ
1
,
通带内
1
,
阻带内
W(\omega)=\left\{\begin{matrix} \delta_2/\delta_1,通带内 \\ 1, 阻带内 \end{matrix}\right.
W(ω)={δ2/δ1,通带内1,阻带内
其中,
δ
1
\delta_1
δ1,
δ
2
\delta_2
δ2的定义如下图所示
定义误差函数:
E
(
ω
)
=
W
(
ω
)
[
H
d
r
(
ω
)
−
H
r
(
ω
)
]
=
W
(
ω
)
[
H
d
r
(
ω
)
−
Q
(
ω
)
P
(
ω
)
]
=
W
(
ω
)
Q
(
ω
)
[
H
d
r
(
ω
)
Q
(
ω
)
−
P
(
ω
)
]
=
W
^
(
ω
)
[
H
d
r
^
(
ω
)
−
P
(
ω
)
]
E(\omega)=W(\omega)[H_{dr}(\omega)-H_r(\omega)]=W(\omega)[H_{dr}(\omega)-Q(\omega)P(\omega)]\\=W(\omega)Q(\omega)[\frac{H_{dr}(\omega)}{Q(\omega)}-P(\omega)]=\hat{W}(\omega)[\hat{H_{dr}}(\omega)-P(\omega)]
E(ω)=W(ω)[Hdr(ω)−Hr(ω)]=W(ω)[Hdr(ω)−Q(ω)P(ω)]=W(ω)Q(ω)[Q(ω)Hdr(ω)−P(ω)]=W^(ω)[Hdr^(ω)−P(ω)]
我们期望找到一组系数
α
(
n
)
\alpha(n)
α(n)(根据具体滤波器类型,可对应到表中的
a
‾
[
n
]
、
b
‾
[
n
]
、
c
‾
[
n
]
、
d
‾
[
n
]
\overline{a}[n]、\overline{b}[n]、\overline{c}[n]、\overline{d}[n]
a[n]、b[n]、c[n]、d[n],并进一步对应到滤波器的冲激响应
h
[
n
]
h[n]
h[n]),使得在要求的频带内
E
(
ω
)
E(\omega)
E(ω)的最大值可以最小化,所以优化问题可描述为
m
i
n
[
max
ω
∈
S
∣
E
(
ω
)
∣
]
(5)
min[\mathop{\max}\limits_{\omega \in S}|E(\omega)|] \tag{5}
min[ω∈Smax∣E(ω)∣](5)
其中,
S
S
S表示最优化的频带集合,通常由通带和阻带组成。下面考虑对
P
(
ω
)
P(\omega)
P(ω)进行形式转换。对于三角函数有
c
o
s
(
2
ω
)
=
2
c
o
s
ω
−
1
、
c
o
s
(
3
ω
)
=
4
c
o
s
3
ω
−
3
c
o
s
ω
.
.
.
.
cos(2\omega)=2cos\omega-1、cos(3\omega)=4cos^3 \omega-3cos\omega....
cos(2ω)=2cosω−1、cos(3ω)=4cos3ω−3cosω....所以
P
(
ω
)
=
∑
n
=
0
L
α
(
n
)
c
o
s
(
ω
n
)
=
∑
n
=
0
L
β
(
n
)
c
o
s
n
ω
P(\omega)=\sum_{n=0}^{L}\alpha(n)cos(\omega n)=\sum_{n=0}^{L}\beta(n)cos^n\omega
P(ω)=∑n=0Lα(n)cos(ωn)=∑n=0Lβ(n)cosnω,即
P
(
ω
)
P(\omega)
P(ω)为
c
o
s
ω
cos\omega
cosω的
L
L
L阶多项式,由于
c
o
s
ω
cos\omega
cosω在
0
<
ω
<
π
0<\omega<\pi
0<ω<π区间上是单调函数,所以可以将
P
(
ω
)
P(\omega)
P(ω)按常规的
L
L
L阶多项式进行处理。所以,优化问题转化为找到一组多项式系数使其满足(5)式的目标。这个问题即为切比雪夫最佳逼近问题。下面以一阶多项式逼近为例,直观说明该优化问题的基本思路。
当多项式为1阶时,优化问题可写成
m
i
n
[
max
a
,
b
∈
[
M
,
N
]
∣
f
(
x
)
−
a
x
−
b
∣
]
(6)
min[\mathop{\max}\limits_{a,b\in[M,N]}|f(x)-ax-b|]\tag{6}
min[a,b∈[M,N]max∣f(x)−ax−b∣](6)
上面的优化问题可以用下图进行形象说明,我们期望期望找到一条直线,使与图中曲线的最大偏差最小。首先画两条平行线将曲线夹在其中,则期望的最佳直线即为这两条平行线的中线,最小误差即为这两条平行线与y轴截距的中点。
首先我们定义偏差点的概念。偏差值
E
E
E定义为满足(6)式的最大偏差的最小值。若
x
0
∈
[
M
,
N
]
x_0\in[M,N]
x0∈[M,N]使得
∣
f
(
x
0
)
−
a
x
0
−
b
∣
=
E
|f(x_0)-ax_0-b|=E
∣f(x0)−ax0−b∣=E,则称
x
0
x_0
x0为正偏差点,若
x
0
∈
[
M
,
N
]
x_0\in[M,N]
x0∈[M,N]使得
∣
f
(
x
0
)
−
a
x
0
−
b
∣
=
−
E
|f(x_0)-ax_0-b|=-E
∣f(x0)−ax0−b∣=−E,则称
x
0
x_0
x0为负偏差点。由上图我们可以总结存在上述最佳逼近直线的充要条件:
01、 f ( x ) f(x) f(x)为连续函数;
02、至少存在3个偏差点,且正负偏差点依次出现。
不加证明地将上述1次多项式的情况推广到本文的
L
L
L次多项式上,可知若存在
P
(
ω
)
P(\omega)
P(ω)为
H
d
r
^
(
ω
)
\hat{H_{dr}}(\omega)
Hdr^(ω)的唯一逼近,则
E
(
ω
)
E(\omega)
E(ω)在
S
S
S上至少存在
L
+
2
L+2
L+2个偏差点,且这些偏差点正负依次出现。该结论被称为交错定理。
对于函数
P
(
ω
)
=
∑
n
=
0
L
α
(
n
)
c
o
s
(
ω
n
)
=
∑
n
=
0
L
β
(
n
)
c
o
s
n
ω
P(\omega)=\sum_{n=0}^{L}\alpha(n)cos(\omega n)=\sum_{n=0}^{L}\beta(n)cos^n\omega
P(ω)=∑n=0Lα(n)cos(ωn)=∑n=0Lβ(n)cosnω,若
c
o
s
(
ω
)
cos(\omega)
cos(ω)进行求导,显然在
0
<
ω
<
π
0<\omega<\pi
0<ω<π区间内至多有
L
−
1
L-1
L−1个导数值为0,对应
L
−
1
L-1
L−1个偏差点。另外,
ω
=
0
\omega=0
ω=0、
ω
=
π
\omega=\pi
ω=π、
ω
=
ω
p
\omega=\omega_p
ω=ωp和
ω
=
ω
s
\omega=\omega_s
ω=ωs通常也可为偏差点,故至多存在
L
+
3
L+3
L+3个偏差点。设这些偏差点对应的频率为
ω
n
\omega_n
ωn,则根据交错定理,可得
W
^
(
ω
n
)
[
H
d
r
^
(
ω
n
)
−
P
(
ω
n
)
]
=
(
−
1
)
n
δ
,
n
=
0
,
1
,
.
.
.
,
L
+
1
\hat{W}(\omega_n)[\hat{H_{dr}}(\omega_n)-P(\omega_n)]=(-1)^n\delta,n=0,1,...,L+1
W^(ωn)[Hdr^(ωn)−P(ωn)]=(−1)nδ,n=0,1,...,L+1
其中,
δ
\delta
δ表示偏差值。将上面的表达式调整为
P
(
ω
n
)
+
(
−
1
)
n
δ
W
^
(
ω
n
)
=
H
d
r
^
(
ω
n
)
→
∑
n
=
0
L
α
(
n
)
c
o
s
(
ω
n
n
)
+
(
−
1
)
n
δ
W
^
(
ω
n
)
=
H
d
r
^
(
ω
n
)
P(\omega_n)+\frac{(-1)^n\delta}{\hat{W}(\omega_n)}=\hat{H_{dr}}(\omega_n)\rightarrow \sum_{n=0}^L \alpha(n)cos(\omega_n n)+\frac{(-1)^n\delta}{\hat{W}(\omega_n)}=\hat{H_{dr}}(\omega_n)
P(ωn)+W^(ωn)(−1)nδ=Hdr^(ωn)→n=0∑Lα(n)cos(ωnn)+W^(ωn)(−1)nδ=Hdr^(ωn)
可将上式写成矩阵的形式:
[
1
c
o
s
(
ω
0
)
.
.
.
c
o
s
(
L
ω
0
)
1
W
^
(
ω
0
)
1
c
o
s
(
ω
1
)
.
.
.
c
o
s
(
L
ω
1
)
1
W
^
(
ω
1
)
⋮
.
.
.
⋮
⋮
1
c
o
s
(
ω
L
+
1
)
.
.
.
c
o
s
(
L
ω
L
+
1
)
(
−
1
)
L
+
1
W
^
(
ω
L
+
1
)
]
[
α
(
0
)
α
(
1
)
⋮
α
(
L
)
δ
]
=
[
H
d
r
^
(
ω
0
)
H
d
r
^
(
ω
1
)
⋮
H
d
r
^
(
ω
L
)
H
d
r
^
(
ω
L
+
1
)
]
(7)
\left[\begin{matrix} 1 & cos(\omega_0) & ... & cos(L\omega_0) & \frac{1}{\hat{W}(\omega_0)} \\ 1 & cos(\omega_1) & ... & cos(L\omega_1) & \frac{1}{\hat{W}(\omega_1)}\\ \vdots & ... & \vdots & \vdots\\ 1 & cos(\omega_{L+1}) & ... & cos(L\omega_{L+1}) & \frac{(-1)^{L+1}}{\hat{W}(\omega_{L+1})} \end{matrix} \right] \left[\begin{matrix} \alpha(0) \\ \alpha(1) \\ \vdots\\ \alpha(L) \\ \delta \end{matrix} \right] = \left[\begin{matrix} \hat{H_{dr}}(\omega_0) \\ \hat{H_{dr}}(\omega_1) \\ \vdots\\ \hat{H_{dr}}(\omega_L) \\ \hat{H_{dr}}(\omega_{L+1}) \end{matrix} \right] \tag{7}
11⋮1cos(ω0)cos(ω1)...cos(ωL+1)......⋮...cos(Lω0)cos(Lω1)⋮cos(LωL+1)W^(ω0)1W^(ω1)1W^(ωL+1)(−1)L+1
α(0)α(1)⋮α(L)δ
=
Hdr^(ω0)Hdr^(ω1)⋮Hdr^(ωL)Hdr^(ωL+1)
(7)
在上面的矩阵中,
w
n
w_n
wn、
δ
\delta
δ、
α
(
n
)
\alpha(n)
α(n)均为未知数,无法进行直接求解。一个直观的思路是先猜一组
w
n
w_n
wn,然后确定
P
(
ω
)
P(\omega)
P(ω)和
δ
\delta
δ,最后计算
E
(
ω
)
E(\omega)
E(ω),如结果不满足要求就重新选择新的偏差点频率,不断迭代求解,但上述求解过程中涉及矩阵求逆,计算量较大,故考虑用Remez算法进行求解。
首先根据矩阵(7)可以求得
δ
\delta
δ的表达式如下:
δ
=
γ
0
H
d
r
^
(
ω
0
)
+
γ
1
H
d
r
^
(
ω
1
)
+
.
.
.
+
γ
L
+
1
H
d
r
^
(
ω
L
+
1
)
γ
0
W
^
(
ω
0
)
−
γ
1
W
^
(
ω
1
)
+
.
.
.
+
(
−
1
)
L
+
1
γ
L
+
1
W
^
(
ω
L
+
1
)
(8)
\delta = \frac{\gamma_0 \hat{H_{dr}}(\omega_0) +\gamma_1 \hat{H_{dr}}(\omega_1)+...+\gamma_{L+1} \hat{H_{dr}}(\omega_{L+1})}{\frac{\gamma_0}{\hat{W}(\omega_0)}-\frac{\gamma_1}{\hat{W}(\omega_1)}+...+\frac{(-1)^{L+1}\gamma_{L+1}}{\hat{W}(\omega_{L+1})}} \tag{8}
δ=W^(ω0)γ0−W^(ω1)γ1+...+W^(ωL+1)(−1)L+1γL+1γ0Hdr^(ω0)+γ1Hdr^(ω1)+...+γL+1Hdr^(ωL+1)(8)
其中
γ
k
=
∏
n
=
0
,
n
≠
k
L
+
1
1
c
o
s
ω
k
−
c
o
s
ω
n
\gamma_k=\prod_{n=0,n\neq k}^{L+1}\frac{1}{cos\omega_k-cos\omega_n}
γk=n=0,n=k∏L+1cosωk−cosωn1
显然,若已知所有极值点频率,则
δ
\delta
δ是可求的。到这里,我们再来重新梳理一下我们想干什么事,又已知什么东西。
首先,从上面的分析我们知道如果目标函数
m
i
n
[
max
ω
∈
S
∣
E
(
ω
)
∣
]
min[\mathop{\max}\limits_{\omega \in S}|E(\omega)|]
min[ω∈Smax∣E(ω)∣]有唯一解,那么一定至少存在
L
+
2
L+2
L+2个频点满足
W
^
(
ω
n
)
[
H
d
r
^
(
ω
n
)
−
P
(
ω
n
)
]
=
(
−
1
)
n
δ
\hat{W}(\omega_n)[\hat{H_{dr}}(\omega_n)-P(\omega_n)]=(-1)^n\delta
W^(ωn)[Hdr^(ωn)−P(ωn)]=(−1)nδ,其中
δ
\delta
δ表示最大误差。但问题是,我们不知道
ω
n
\omega_n
ωn具体是多少,所以只能随便猜几个。如果运气爆棚,一次就猜准了,那么根据(8)式算出来的
δ
\delta
δ就是真实的我们想要的最大误差,此时对于所有
ω
∈
S
\omega \in S
ω∈S(注意
ω
\omega
ω是连续的)都会满足
∣
E
(
ω
)
∣
≤
δ
|E(\omega)|\leq \delta
∣E(ω)∣≤δ。但是如果极值频率点选择错误,计算出来的
δ
′
\delta'
δ′将比真实
δ
\delta
δ小(因为真实
δ
\delta
δ本身表示的就是最大误差)。所以这个时候误差
∣
E
(
ω
)
∣
|E(\omega)|
∣E(ω)∣曲线上一定存在一些局部峰值点使得
∣
E
(
ω
)
∣
>
δ
|E(\omega)|>\delta
∣E(ω)∣>δ。基于这个分析,我们就可以不断去迭代极值频率点,直到找到所有正确的极值频率点。
将整体操作步骤总结如下:
(1)设置滤波器类型和阶数
N
N
N,并根据表2计算
L
L
L。设置通带/阻带纹波大小
δ
1
\delta_1
δ1/
δ
2
\delta_2
δ2,设置理想频响
H
d
r
(
ω
)
H_{dr}(\omega)
Hdr(ω);
(2)基于
L
+
2
L+2
L+2个极值频率点
ω
0
,
.
.
.
ω
L
+
1
\omega_0,...\omega_{L+1}
ω0,...ωL+1(第一次随机选择),并根据(8)计算
δ
\delta
δ;
(3)对
P
(
ω
n
)
P(\omega_n)
P(ωn)进行插值(算法中用的是拉格朗日插值) 得到
P
(
ω
)
P(\omega)
P(ω),并计算
E
(
ω
)
=
W
^
(
ω
)
[
H
d
r
^
(
ω
)
−
P
(
ω
)
]
E(\omega)=\hat{W}(\omega)[\hat{H_{dr}}(\omega)-P(\omega)]
E(ω)=W^(ω)[Hdr^(ω)−P(ω)],显然插值的目的是用密集采样点来逼近连续的情况,一般插值到
16
N
16N
16N点就足够;
(4)找出
E
(
ω
)
E(\omega)
E(ω)中满足
∣
E
(
ω
)
∣
≥
δ
|E(\omega)|\geq \delta
∣E(ω)∣≥δ的全部局部极值点(特别注意这边找的是极值点,不是一般的点,因为对于所有预设频率均能满足
∣
E
(
ω
n
)
∣
=
δ
|E(\omega_n)|=\delta
∣E(ωn)∣=δ,但这些预设频率点不一定是极值点);
(5)如果找到满足步骤(4)要求的极值点,且找到的极值频率点与预设的极值频率点不一致,则将预设极值频率点替换成找到的极值频率点,若找到的极值频率点数大于
L
+
2
L+2
L+2个,则取前
L
+
2
L+2
L+2个极值频率点替换预设极值频率点。重新跳转到步骤(2)进行迭代处理。
(6)当找到的极值频率点数为
L
+
2
L+2
L+2且与预设极值频率点一致时,跳出迭代过程。
(7)当找到正确的极值频率点后,步骤(3)插值得到的
P
(
ω
)
P(\omega)
P(ω)即为正确的,据此可得所设计的滤波器频响的实部为
H
r
(
ω
)
=
Q
(
ω
)
P
(
ω
)
H_r(\omega)=Q(\omega)P(\omega)
Hr(ω)=Q(ω)P(ω),结合其线性相位的虚部,即可得到完整的滤波器频响。
【附录】
- 拉格朗日插值