基于MATLAB的FIR低通滤波器文档格式.docx
《基于MATLAB的FIR低通滤波器文档格式.docx》由会员分享，可在线阅读，更多相关《基于MATLAB的FIR低通滤波器文档格式.docx（25页珍藏版）》请在冰点文库上搜索。
∙矩形窗w（n）=RN（n）；
∙Hanning窗
；
∙Hamming窗
；
∙Blackmen窗
窗函数法设计的步骤
∙确定数字滤波器的性能要求：
临界频率wp和ws，滤波器单位脉冲响应长度N；
∙根据性能要求，合理选择单位脉冲响应h（n）的奇偶对称性，从而确定理想频率响应Hd（ejω）的幅频特性和相频特性；
∙选择适当的窗函数w（n），根据h（n）=hd（n）w（n）求所需设计的FIR滤波器单位脉冲响应；
∙求H（ejω），分析其幅频特性，若不满足要求，可适当改变窗函数形式或长度N，重复上述设计过程，以得到满意的结果。
2）频率采样法设计的原理
频率采样法是从频域出发，将给定的理想频率响应Hd（ejω）加以等间隔采样
然后以此Hd（k）作为实际FIR数字滤波器的频率特性的采样值H（k）,即令
由H（k）通过IDFT可得有限长序列h（n）
将上式代入到Z变换中去可得
其中Φ（ω）是内插函数
3）等波纹最佳逼近法设计FIR滤波器的原理
FIR滤波器的优化设计是按照最大误差最小化准则，使所设计的频响与理想频响之间的最大误差，在通带和阻带范围均为最小，而且是等波动逼近的。
为了简化起见，在优化设计中一般将线性相位FIR滤波器的单位脉冲响应h（n）的对称中心置于n=0处，此时，线性相位因子α=0。
当N为奇数，且N=2M+1，则
如逼近一个低通滤波器，这里M，ωp和ωs固定为某个值。
在这种情况下有
定义一逼近误差函数：
E（ω）为在希望的滤波器通带和阻带内算出的误差值，W（ω）为加权函数。
k应当等于比值δ1/δ2，δ1为通带波动，δ2为阻带波动。
在这种情况下，设计过程要求|E（ω）|在区间
的最大值为最小，它等效于求最小δ2。
根据数学上多项式逼近连续函数的理论，用三角多项式逼近连续函数，在一定条件下存在最佳逼近的三角多项式，而且可以证明这个多项式是唯一的。
这一最佳逼近定理通常称作交替定理。
在逼近过程中，可以固定k，M，ωp，ωs而允许改变δ2，按照交替定理，首先估计出（M+2）个误差函数的极值频率点{ωi}，i=0,1,...,M+1，共计可以写出（M+2）个方程
式中ρ表示峰值误差。
一般仅需求解出ρ，接着便可用三角多项式找到一组新的极值频率点，并求出新的峰值误差ρ。
依此反复进行，直到前、后两次ρ值不变化为止，最小的ρ即为所求的δ2。
（3）三种方法设计方法的MATLAB程序和仿真图
【1】窗函数法设计FIR低通滤波器的MATLAB程序
clearall
%用窗函数法设计FIR低通滤波器
Wc=0.4*pi;
N=61;
w1=boxcar（N）;
w2=bartlett（N）;
w3=hamming（N）;
w4=hanning（N）;
w5=blackman（N）;
n=1:
1:
61;
hd=sin（Wc*（n-31））./（pi*（n-31））;
%理想低通滤波器单位脉冲响应函数
hd（31）=Wc/pi;
h1=hd.*rot90（w1）;
h2=hd.*rot90（w2）;
h3=hd.*rot90（w3）;
h4=hd.*rot90（w4）;
%rot90（w）是将窗函数的列向量转换成行向量
h5=hd.*rot90（w5）;
%hd.*rot90（w）是窗函数截取的一段单位脉冲响应函数
[H1,W]=freqz（h1）;
[H2,W]=freqz（h2）;
[H3,W]=freqz（h3）;
[H4,W]=freqz（h4）;
[H5,W]=freqz（h5）;
%截取的一段单位脉冲响应函数的幅频特性
figure
（1）;
subplot（2,2,1）;
stem（1:
N,w1）;
title（'
矩形窗的时域波形和幅频特性'
）;
subplot（2,2,2）;
plot（W/pi,20*log10（abs（H1）/max（H1）））;
xlabel（'
频率/Hz'
ylabel（'
振幅/dB'
gridon;
subplot（2,2,3）;
plot（W,20*log10（abs（H1）））;
矩形窗设计的低通滤波器的幅频特性和相频特性'
[h1,w1]=freqz（h1,1,129,2）;
subplot（2,2,4）;
plot（w1,unwrap（angle（h1）））;
w/pi'
角弧度/rad'
figure
（2）;
N,w2）;
三角窗的时域波形和幅频特性'
plot（W/pi,20*log10（abs（H2）/abs（max（H2））））;
plot（W,20*log10（abs（H2）））;
三角窗设计的低通滤波器的幅频特性和相频特性'
[h2,w2]=freqz（h2,1,129,2）;
plot（w2,unwrap（angle（h2）））;
figure（3）;
N,w3）;
海明窗的时域波形和幅频特性'
plot（W/pi,20*log10（abs（H3）/abs（max（H3））））;
plot（W,20*log10（abs（H3）））;
海明窗设计的低通滤波器的幅频特性和相频特性'
[h3,w3]=freqz（h2,1,129,2）;
plot（w3,unwrap（angle（h3）））;
figure（4）;
N,w4）;
汉宁窗的时域波形和幅频特性'
plot（W/pi,20*log10（abs（H4）/abs（max（H4））））;
plot（W,20*log10（abs（H4）））;
汉宁窗设计的低通滤波器的幅频特性和相频特性'
[h4,w4]=freqz（h4,1,129,2）;
plot（w4,unwrap（angle（h4）））;
figure（5）;
N,w5）;
布莱克曼窗的时域波形和幅频特性'
plot（W/pi,20*log10（abs（H5）/abs（max（H5））））;
plot（W,20*log10（abs（H5）））;
布莱克曼窗设计的低通滤波器幅频特性和相频特性'
[h5,w5]=freqz（h5,1,129,2）;
plot（w5,unwrap（angle（h5）））;
其执行的图形：
【2】频率采样法设计FIR低通滤波器
（1）已知FIR低通滤波器通带截止频率wp，阻带截止频率ws，通带波纹最大为0.04，阻带波纹最大为0.02,滤波器的阶数N通过remezord函数估算,过渡带的采样取值T时的MATLAB程序：
%用频率采样法设计FIR低通滤波器
wp=0.2*pi;
ws=0.3*pi;
%理想低通滤波器的截止频率
T=0.5;
%过渡带采样取值
Bt=[wp/piws/pi];
a=[10];
dev=[0.040.02];
%给出滤波器的参数
[nf0a0w]=remezord（Bt,a,dev）;
N=n;
%N为奇数
alpha=（N-1）/2;
k=0:
N-1;
wc=（wp+ws）/2;
m=fix（wc*N/（2*pi）+1）;
%通带采样点数
Hk=[ones（1,m）,T,zeros（1,N-2*m-1）,T,ones（1,m-1）];
%幅度采样向量
k1=0:
floor（alpha）;
k2=floor（alpha+1）:
phai=[-alpha*（2*pi）/N*k1,alpha*（2*pi）/N*（N-k2）];
%相位采样向量
Hdk=Hk.*exp（j*phai）;
%构造频域采样向量
h=ifft（Hdk,N）;
%进行反离散傅里叶变化
[h1,w1]=freqz（h,1,256,1）;
h1=20*log10（abs（h1））;
%画出FIRDF的单位脉冲响应
plot（n,h1）;
stem（k,h,'
k.'
axis（[0,N-1,1.1*min（real（h））,1.1*max（real（h））]）;
%使横纵坐标成比例
n'
ylabel（'
h（n）'
grid;
频率采样法设计FIR低通滤波器的单位脉冲响应h（n）'
%画出FIRDF的低通的幅频特性
plot（w1,h1）;
频率（×
πrad/sample）'
幅度（dB）'
gridon;
频率采样法设计FIR低通滤波器的幅频特性'
（2）已知通带截止频率wp，过渡带宽度Bt，过渡带的采样取值T时，由滤波器阻带最小衰减Rs决定采样点的个数m值的MATLAB程序：
%用频率采样法设计FIR低通滤波器
%通带截止频率
Bt=0.1*pi;
%过渡带宽度
rs=input（'
rs='
%阻带衰减
ifrs<
=54
m=1;
T=input（'
T='
N=ceil（（m+1）*2*pi/Bt）;
N=N+mod（N+1,2）;
%使N为奇数
Np=fix（wp/（2*pi/N））;
%Np+1为通带采样点数
Ns=N-2*Np-1;
%阻带采样点数
Hk=[ones（1,Np+1）,zeros（1,Ns）,ones（1,Np）];
Hk（Np+2）=T;
Ak（N-Np）=T;
%过度采样
elseifrs<
=75
m=2;
T1=input（'
T1='
T2=input（'
T2='
Hk（Np+2）=T1;
Hk（Np+3）=T2;
Ak（N-Np）=T1;
Ak（N-Np+1）=T2;
else
m=3;
T3=input（'
T3='
Hk（Np+2）=T1;
Hk（Np+4）=T3;
Ak（N-Np）=T1;
Ak（N-Np+2）=T3;
end
thetak=-pi*（N-1）*（0:
N-1）/N;
Hdk=Hk.*exp（j*thetak）;
%频域采样向量
hn=real（ifft（Hdk））;
%时域取实部，忽略计算误差引起的虚部
Hw=fft（hn,1024）;
%计算频率响应函数
wk=2*pi*[0:
1023]/1024;
Hgw=Hw.*exp（j*wk*（N-1）/2）;
%计算幅度响应函数
figure
（1）;
stem（k,hn,'
axis（[0,N-1,1.1*min（real（hn））,1.1*max（real（hn））]）;
%使横纵成比例
xlabel（'
title（'
figure
（2）;
plot（wk,20*log（abs（Hgw）））;
grid;
%画出FIRDF的幅频响应
【3】等波纹最佳逼近法设计FIR低通滤波器
（1）已知通带截止频率wp，阻带截止频率ws，模拟信号的采样频率fs，通带最大衰减rp，阻带最小衰减rs时的MATLAB程序：
的MATLAB程序：
%用等波纹最佳逼近法设计FIR低通滤波器
fs=input（'
fs='
%对模拟信号的采样频率
rp=input（'
rp='
%通带最大衰减
%阻带最小衰减
f=[wp*fs/（2*pi）,ws*fs/（2*pi）];
m=[1,0];
dat1=（10^（rp/20）-1）/（10^（rp/20）+1）;
dat2=10^（-rs/20）;
%求出通带波纹幅度dat1和阻带波纹幅度dat2
rip=[dat1,dat2];
[M,fo,mo,w]=remezord（f,m,rip,fs）;
hn=remez（M,fo,mo,w）;
hw=fft（hn,512）;
w=0:
511;
w=2*pi*w/512;
plot（w,20*log（abs（hw）））;
gridon;
%画出FIRDF的幅频特性
等波纹最佳逼近法设计FIR低通滤波器的幅频特性'
（2）已知通带截止频率wp，阻带截止频率ws，通带波纹最大为0.04，阻带波纹最大为0.02时的MATLAB程序：
（要求给出滤波器的阶数n值）
%用等波纹最佳逼近法设计FIR低通滤波器
n=input（'
n='
f=[0wp/piws/pi1];
a=[1100];
w=[0.040.02];
hn=remez（n,f,a,w）;
[H,W]=freqz（hn）;
plot（W,20*log10（abs（H）））;
[h,w]=freqz（hn,1,129,2）;
plot（w,unwrap（angle（h）））;
%画出FIRDF的相频特性
等波纹最佳逼近法设计FIR低通滤波器的相频特性'
（4）功能说明
稳定和线性相位特性是FIR滤波器最突出的优点。
FIR数字滤波器的优点在于它可以做成具有严格线性相位,而同时可以具有任意的幅度特性;
它的传递函数没有极点;
这保证了设计出的FIR数字滤波器一定是平稳的。
因此，FIR滤波器的功能是：
很容易实现具有严格线性相位的系统,使信号经过处理后不产生相位失真,舍入误差小,而且稳定。
（5）实现步骤
1．对五种窗函数设计的低通滤波器进行比较，滤波器的好差的指标是：
Rn——窗函数的幅频函数的最大旁瓣的最大值相对主瓣最大值的衰减值（越小越好）
Bg——低通滤波器的过渡带宽度=窗函数的幅度特性的主瓣宽度（越窄越好）
Rs——低通滤波器的阻带最小衰减（越小越好）
当窗口长度N值相等时：
矩形窗：
Rn，Bg，Rs三个参数都较大，形成的低通滤波器的截断效应影响较大，从而影响滤波器的通带的平稳性和阻带的衰减。
三角窗：
相对矩形窗各参数有所改善，但仍不是最佳的方法。
汉宁窗：
旁瓣互相抵消，能量集中在主瓣。
这样使RnBgRs各参数值相对较好，能较好的实现低通滤波器的设计。
海明窗：
能量集中于主瓣，主瓣能量约占99.96%。
这样使RnBgRs更接近于理想低通滤波器的设计。
Rn以相当的小，截断效应影响小。
布莱克曼窗：
旁瓣进一步互相抵消，能量更集中与主瓣。
比海明窗稍好一些。
2.比较窗函数法的程序、频率采样法和等波纹逼近法中的程序，同样是设计一个N阶的FIR低通数字滤波器，可以看出：
（1）窗函数法在阶数N值较低时，阻带特性不满足设计要求，只有滤波器阶数较高时，使用海明窗和凯塞窗可以达到阻带衰耗要求；
（2）频率采样法偏离设计指标最明显，阻带衰减最小，而且设计的程序比采用窗函数法复杂。
只有适当选取过渡带样点值，才会取得较好的衰耗特性；
（3）利用等波纹逼近法则的设计可以获得最佳的频率特性和衰耗特性，具有通带和阻带平坦，过渡带窄等优点。
四、参考文献
1.数字信号处理【西安电子科技大学出版社（第三版）高西全丁玉美编著】
2.MATLAB及其在理工课程中的应用指南【西安电子科技大学出版社（第三版）陈怀深编著】
3.仿真软件教程【清华大学出版社常华袁刚常敏嘉编著】
五、实习体会
FIR滤波器设计时选择有限长度的h（n），使频率响应函数H（ejω）满足技术指标要求。
我本次实习采用三种设计方法：
窗函数法、频率采样法和切比雪夫等波纹逼近法。
1.在设计完低通滤波器后，从中可以看出，理想低通滤波器的单位脉冲响应加窗处理后，频率响应函数与原理想低通频率响应函数的差别有以下两点：
（1）在理想特性不连续点ω=ωc附近形成过渡带。
过渡带的宽度近似等于WRg（ω）主瓣宽度4π/N。
（2）通带内产生了波纹，最大的峰值在ωc－2π/N处。
阻带内产生了余振，最大的负峰在ωc+2π/N处。
通带与阻带中波纹的情况与窗函数的幅度谱有关，WRg（ω）旁瓣幅度的大小直接影响Hg（ω）波纹幅度的大小。
2.加窗后，由于窗函数的影响产生截断效应，下面是对截断效应的介绍和采用怎样的方法减小截断效应的影响。
对理想低通滤波器的单位脉冲响应hd（n）用窗函数截断后，在频域通带和阻带都产生的波纹，称为截断效应。
这种效应直接影响滤波器的性能。
通带内的波纹影响滤波器通带的平稳性，阻带内的波纹影响阻带内的衰减，这样可能使最小衰减不满足技术指标要求。
一般滤波器要求过渡带愈窄愈好。
为减少截断效应的影响，可以加大窗函数的长度（加大N）。
但N加大，幅度加高，旁瓣也加高，主瓣和旁瓣幅度的相对值不变；
另一方面，N加大，WRg（ω）的主瓣和旁瓣宽度变窄，波动的频率加快。
因此加大N只能使Hg（ω）过渡带变窄，并不是减小截断效应的有效方法。
减小截断效应最有效的方法是采用合适的窗函数。
3.通过用这三种方法去实现低通滤波器的设计，可以看出对于同样的技术指标，用等波纹最佳逼近法设计的滤波器阶数低，其原因是：
①用窗函数法设计的滤波器，如果在阻带截止频率附近刚好满足，
升级会员 