高级数字滤波器设计及Matlab实现

高级数字滤波器设计及Matlab 实现
利用Parks-McClellan 算法设计线性相位FIR 滤波器
一、
算法原理
长度为2n+1的线性相位数字滤波器的传输函数为：20
()n
k
k
k G z h
Z
-==
∑。

当
Z=exp(j T ω)=exp(2j F π)时，可得到频率响应：
()exp(2)cos 2n
k k G F j nF d k F ππ==-∑
exp(2)()j F H F π=-
其中2n k k d h -=，k=0，……，n-1，0n d h =
max ()()()W F D F H F -
对于一个理想的低通滤波器上式中的H(F)可以表示为：
误差加权函数：
允许设计者自己给定通带和阻带内的误差范围。

令p s B B A = ，设计长度为2n+1的线性相位低通滤波器只要找到k d 并使得m a x ()()()W F D F H F -最小。

设空间A 是[0,1/2]的封闭子空间，为了使0
()cos 2n
k
k H F d
kF
π==
∑是D(F)
在A 上唯一的最佳逼近，加权误差方程()()[()()]E F W F D F H F =-在A 上至少要有n+2个交错点。

因此1()()i i E F E F E -=-=±，011,n i F F F F A +<∈ ，
max ()E E F =。

算法的流程如图1所示。

对于给定的n+2个频率点，需要计算n+2个方程：
()(()())(1)k k
k k W F H F D F ρ-=--
写成矩阵的形式就是：
图 1.
0000001
1
1
1
1011
1
1
1
11
1cos 2cos 4cos 2()()11cos 2cos 4cos 2()()()(1)1cos 2cos 4cos 2()n n n n n n n F F nF W F d D F d F F nF D F W F d D F F F nF W F ππππππρπππ++++++⎡⎤
⎢⎥⎡⎤
⎢⎥
⎡⎤
⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥-⎣⎦
⎢⎥⎢⎥⎣
⎦
通过该方程组可得：
其中：
利用拉格朗日插值公式可得：
这里
利用求得的H(F)求出误差函数E(F)。

如果对所有的频率都有()E F ρ≤，说明ρ是纹波极值，交错频率点121,n F F F + 是交错频率点。

若存在某些频率使得()E F ρ>，说明初始交错点组中的某些点需要交换。

对于上次确定的121,n F F F + 中每一点，都检查其附近是否存在某一频率
()E F ρ>，如果有再在该点附近找出局部极值点，并用该，点代替原来的
点。

待所有点都检查过后，得到新的交错点组。

然后重新求出ρ，H(F)，E(F)再做一次比较。

重复上述步骤，ρ最后收敛到自己的上限此时H(F)就是D(F)的最佳逼近。

由最后一组121,n F F F + 算H(F)。

二、
编程仿真
设计一个低通滤波器，Fp 为500Hz ，Fs 为600Hz ，通带波纹为3dB,阻带波纹为40dB 。

rp = 3; rs = 40; fs = 2000; f = [500 600]; a = [1 0];
dev = [(10^(rp/20)-1)/(10^(rp/20)+1) 10^(-rs/20)]; [n,fo,ao,w] = firpmord(f,a,dev,fs); b = firpm(n,fo,ao,w); freqz(b,1,1024,fs);
title('Lowpass Filter Designed to Specifications');
010*******
4005006007008009001000
-1500
-1000
-500
Frequency (Hz)
P h a s e (d e g r e e s )
010*******
4005006007008009001000
-100
-50
50Frequency (Hz)
M a g n i t u d e (d B )
Lowpass Filter Designed to Specifications
三、 参考文献
[1] T.W. Parks and J.H. McClellan, “Chebyshev approximation for nonrecursive digital filters with linear phase,’’ IEEE Trans. Circuit Theory , vol. CT-19, no. 2, pp. 189–194, 1972.
[2]丁玉美 高西全 《数字信号处理》（第二版） 西安电子科技大学出版社。