发布时间 : 星期一 文章数字信号处理上机指导 - 1dl更新完毕开始阅读a620bd4b852458fb770b56fc
function [bz,az]=fun8 %求模拟滤波器参数
wp=100*2*pi; %通带截止频（率弧度/秒） ws=200*2*pi; %阻带截止频率（率弧度/秒） Rp=2; Rs=15;
Fs=500; %采样频率 Ts=1/Fs;
[N,Wn]=buttord(wp,ws,Rp,Rs,'s'); %选择低通滤波器的最小阶数N
[Z,P,K]=buttap(N); %创建零极点模型的巴特沃思型模拟低通滤波器 [Bap,Aap]=zp2tf(Z,P,K); %把滤波器零极点模型转化为传递函数模型
[b,a]=lp2lp(Bap,Aap,Wn); %把模拟滤波器原型转换成截止频率为的低通滤波器 [bz,az]=bilinear(b,a,Fs); %用双线性变换法实现模拟滤波器到数字滤波器的转换 [H,W]=freqz(bz,az);
plot(W*Fs/(2*pi),20*log10(abs(H))); %绘制频率响应曲线 grid
xlabel('频率/Hz');
ylabel('幅频响应(dB)') 运行结果为右上图。
§6-5 FIR滤波器的设计与实现
相对于无限冲激响应数字滤波器，有限冲激响应数字滤波器可实现具有精确线性相位、总是稳定的、过渡过程有限、硬件易实现等优点，但所需滤波器的阶数较高，其延迟也较大。
N?1有限长单位冲击响应（FIR）数字滤波器系统的传递函数模型为：H(z)??h(n)zn?0?n，
其设计方法有窗函数法和频率采样法等，在MATLAB的数字信号处理工具箱中提供了fir1等设计函数。fir1是采用经典窗函数法设计线性相位FIR数字滤波器的函数，且具有标准低通、带通、高通和带阻等类型。
语法格式： B = fir1 (n , Wn , ’ ftype ’ , window)
其中n为FIR滤波器的阶数即窗口长度为n+1，对于高通、带阻滤波器 n只取偶数。Wn为对Nyquist频率（即大于抽样频率的二分之一）归一化的滤波器截止频率，取值范围0~1，即当Wn=1时，滤波器实际截止频率等于Nyquist频率。对于带通、带阻滤波器，Wn=[W1,W2]，且W1 例6-5 用窗函数法设计一个线性相位FIR低通滤波器，性能指标：通带截止频率Wp=0.2π，阻带截止频率Ws=0.3π，阻带衰减不小于40dB，通带衰减不大于3dB。设抽样频率Fs=1 Hz. 由阻带衰减不小于40dB的情况下，选择汉宁窗将比较经济。用汉宁窗生成线性相位FIR低通滤波器的实现程序如下： function [b,k]=fun7 wp=0.2*pi; ws=0.3*pi; wdelta=ws-wp; %计算过渡带宽 N=ceil(8*pi/wdelta); ?il：沿正无穷方向取整，估算窗口长度 Wn=(0.2+0.3)*pi/2; %估算窗口截止频率 - 20 - Wn=Wn/pi; %对Nyquist频率(即πrad/s)归一化 b=fir1(N,Wn,hanning(N+1)); %用汉宁窗生成线性相位FIR低通滤波器 freqz(b,1,512); %绘出此系统的频响特性 100n=0:511; W1=exp(-j*wp*n); 0W2=exp(-j*ws*n); -100HW1=20*log10(abs(sum(b.*W1))) -200HW2=20*log10(abs(sum(b.*W2))) 00.10.20.30.40.50.60.70.8Normalized Frequency (?? rad/sample)if(abs(HW1)<3)&(abs(HW2)>40) 0k=1; %满足技术指标 -1000else k=-1 %不满足技术指标 -2000end Magnitude (dB)Phase (degrees)-30000.91常用窗函数的MATLAB命令如下： 窗函数 命令 矩形窗 boxcar(n) 三角形窗 triang(n) 汉宁窗 00.10.20.30.40.50.60.70.8Normalized Frequency (?? rad/sample)0.91海明窗 布拉克曼窗 凯泽窗 hanning(n) hamming(n) blackman(n) kaiser(n,β) 注：窗函数的使用方法相同。n均为窗的长度，β为凯泽窗的形状参数，它们的取值参考教材相应的表格。 例6-6 设抽样频率为Fs=1000Hz，已知原信号为x=sin(2π×80t)+2sin(2π×140t)，由于某种原因，信号被白噪声污染，实际获得的信号为xn=x+rand(size(t))，要求设计一个FIR滤波器恢复出原始信号。 滤波器设计要求； 频带(Hz) [0,65] [75,85] [95,125] [135,145] [155,500] 幅度 0 1 0 1 0 function lvbo (见Fig.1和Fig.3) %生成实际信号xn function lvbo_cls (见Fig2和Fig.3) Fs=1000; %生成实际信号xn Fs=1000; t=0:1/Fs:2; t=0:1/Fs:2; x=sin(2*pi*80*t)+sin(2*pi*140*t); x=sin(2*pi*80*t)+sin(2*pi*140*t); xn=x+randn(size(t)); xn=x+randn(size(t)); %生成一个阶数为100的多带滤波器 %生成一个阶数为100的多带滤波器 n=100; n=100; w1=[75,145]/(Fs/2);w2=[85,135]/(Fs/2); f=[0 0.13 0.15 0.17 0.19 0.25 0.27 0.29 0.31 1]; %生成以w1为通带的带通滤波器 m=[0 0 1 1 0 0 1 1 0 0] b1=fir1(n,w1,hamming(n+1)); b=firls(n,f,m) % 生成以w2为阻带的带阻滤波器 [H,W]=freqz(b,1,512,2); b2=fir1(n,w2,'stop',hamming(n+1)); %绘系统的频响特性 %计算两个系统的频率响应 figure(1) plot(W,abs(H)); [H1,W]=freqz(b1,1,512,2); grid [H2,W]=freqz(b2,1,512,2); xlabel('归一化频率(Nyquist频率=1)') %绘系统的频响特性 ylabel('幅度') figure(1) %对实际信号xn进行滤波 plot(W,abs(H1),'r', W,abs(H2)); xo=filter(b,1,xn); grid %检查滤波效果 xlabel('归一化频率(Nyquist频率=1)') figure(2); ylabel('幅度') nn=500:750; %对实际信号xn进行滤波 tt=nn/Fs; y1=filter(b1,1,xn); subplot(2,1,1); y=1.5*filter(b2,1,y1); plot(tt,x(nn),tt,xo(nn),'r'); legend('原始信号','滤波信号'); %检查滤波效果 grid figure(2); subplot(2,1,2); nn=500:750; plot(tt,xn(nn)); legend('污染信号');grid; xlabel('时间(秒)'); - 21 - 1.4tt=nn/Fs; subplot(2,1,1); plot(tt,x(nn),tt,y(nn),'r'); legend('原始信号','滤波信号'); grid subplot(2,1,2); plot(tt,xn(nn)); legend('污染信号');grid; xlabel('时间(秒)'); 1.210.8幅度0.60.40.200 1.40.10.20.30.40.50.60.7归一化频率(Nyquist频率=1)0.80.91Fig.1 FIR滤波器的幅频特性 01.210.8幅度0.60.40.20 0.1 0.2 0.3 0.40.50.60.7归一化频率(Nyquist频率=1)0.80.91Fig.2 FIR滤波器的幅频特性 Fig.3 滤波效果对比 实验一 序列线性卷积、圆周卷积的计算及其关系 一、实验目的 1．实现序列线性卷积、圆周卷积的计算； 2．研究序列线性卷积、圆周卷积的关系。 二、实验内容 设有一离散因果系统，其输入输出关系由以下差分方程确定 y(n)-0.5y(n-1)=x(n)+0.5x(n-1) 1．求该系统的单位抽样响应； 2．利用卷积和求系统在x(n)=e j3n激励下的零状态响应y(n)。 三、实验原理 1．当一个离散时间序列信号x(n)通过一个线性离散时间系统（其单位冲击响应为h(n)）后的输出y(n)，用数学语言表达为x(n)与h(n)的线性卷积，即 ? y(n)= x(n)* h(n)=?h(k)x(n?k) k?02．两个有限长的序列x1（n）、N1和x2（n）、N2的N点圆周卷积z（n） - 22 - N?1Nx2（n）=z（n）= x1（n）○?x1(k)x2(n?k)RN(n) k?03． 若x(n)的长度为N1，h(n)的长度为N2，在N＞N1+ N2－1的条件下，令 ?0x1(n)???x(n)N1?n?N?10?n?N1?1，x2(n)???0?h(n)N2?n?N?10?n?N2?1 Nx2（n） 则 x1（n）*x 2（n）=x1（n）○ 即可用序列的圆周卷积计算序列的线性卷积。 四、实验步骤 编写MATLAB应用程序，以实现如下功能 1．描述系统； 2．计算系统的单位抽样响应，并绘出频率特性； 3．直接计算x(n)与h(n)的线性卷积，并绘出频率特性； 4．在N＞N1+ N2－1、N=N1+ N2－1、N＜N1+ N2－1三种情况下，计算x(n)与h(n)的圆周卷积，并绘出频率特性。 五、实验要求 1．了解MATLAB软件，认真学习并练习每一个例题。 2．针对实验内容，首先进行理论分析，给出预测结果；再编写程序，上机运行，给出仿真结果；进行结果验证，如有不同，查找原因，解决问题。 六、程序示例 function [hn,yn]=shiyan1(M,N) M=160; N=180; b=[1,0.5];a=[1,-0.5]; %线性离散系统的传递函数法的描述 n=0:M-1; %令M= N2（=N1） xh=[1,zeros(1,M-1)]; %创建单位抽样激励 hn=filter(b,a,xh); %求该系统的单位抽样响应hn figure(1); stem(n,hn); %绘制系统的单位抽样响应hn的离散序列图 axis([0,M,0,1.1]); xn=exp(j*3*n); %创建激励x(n)=e j3n yn=conv(xn,hn); %直接计算线性卷积y(n) ynn=circonv(xn,hn,N); %调用函数circonv和cirshiftd程序在（18页例6-1）计算xn与hn %N点的圆周卷积 nh=0:1:length(yn)-1; ny=0:1:length(ynn)-1; figure(2) subplot(3,1,1); stem(n,abs(xn)) %绘制激励x(n)=e j3n幅度的离散序列图 axis([0,M,0,1.1]); title('输入激励幅度的离散序列图') subplot(3,1,2); stem(nh,abs(yn)); %绘制零状态响应yn幅度的离散序列图 axis([0,N,0,1.1]); title('线性卷积幅度的离散序列图') subplot(3,1,3); stem(ny,abs(ynn),'r'); %绘制零状态响应ynn幅度的离散序列图 axis([0,N,0,1.1]); - 23 -