数字信号处理课程设计.docx

数字信号处理课程设计一、课程设计目的1、掌握matlab用于数字滤波器的设计技术和方法；2、掌握fir滤波器的设计方法和步骤；3、掌握iir滤波器的设计方法和步骤。二、课程设计任务及步骤题目1、双线性变换法设计数字滤波器 （1） 采用双线性变换法设计一个巴特沃斯数字低通滤波器，要求：wp=0.25，rp=1db，ws=0.4，as=15db，滤波器的采样频率为fs=100hz。设计步骤：1、用matlab信号处理工具箱函数buttord，butter是巴特沃斯滤波器设计函数，该程序中调用以下格式： n,wc=buttord(wp,ws,rp,as)该格式用于计算巴特沃斯数字滤波器的阶数n和3db截止频率wc，调用参数wp和ws分别为巴特沃斯滤波器的通带边界频率和阻带边界频率的归一化值，0wp1,0ws1, rp，as分别为通带最大衰减和阻带最小衰减（db）。 b,a=butter(wp,ws,rp,as,s)该格式用于计算巴特沃斯数字滤波器系统函数分子和分母多项式的字数向量b,a。调用参数n,和wc分别为巴特沃斯数字滤波器的阶数和3db截止频率的归一化值。开始2、滤波器设计流程图读入要求设计指标的wp,ws,as,rp,fs用双极性变换法转换wp,ws为wp1,ws1调用工具箱函数buttord, butter求n,wc,b,a调用bilinear将模拟转换为数字滤波器参数调用子程序(函数)计算h，w，pha结束调用子程序(函数)绘出h，pha，损耗曲线 3、程序实现代码如下：%双极性变换法设计iir切巴特沃斯数字低通滤波器close all;clear all;%=要求设计指标=fs=100;t=1/fs; rp=1; as=15; wp=0.25; ws=0.4; %=wp1=2/t*tan(wp*pi/2); ws1=2/t*tan(ws*pi/2); %双极性变换法边界频率转换指标n,wc=buttord(wp1,ws1,rp,as,s) %计算滤波器阶数和3db截止频率wcb,a =butter(n,wc,s); %计算滤波器系统函数分子、分母多项式系数bz,az=bilinear(b,a,fs); %用双极性变换法将模拟滤波器转换成数字滤波器h,w =freqz(bz,az); %计算h(z)的幅频响应,freqz(bz,az,计算点数,采样频率)subplot(2,2,1);plot(w/pi,abs(h); xlabel(w/pi); ylabel(幅度(db) h1=10(-rp/20);%通带波纹幅度 rp=-20lg(h1)h2=10(-as/20);%阻带波纹幅度 as=-20lg(h2)set(gca,xtick,0,0.25,0.4,1); set(gca,ytick,0,h2,h1,1); text(0.25,h1,-通带最大衰减,fontsize,8);%文本注释text(0.4,h2,-阻带最小衰减,fontsize,8); %文本注释axis(0,1,0,1.1);title(巴特沃斯低通滤波器);grid onsubplot(2,2,2);pha=angle(h)/pi;plot(w/pi,pha);grid ontitle(相位响应曲线);xlabel(w/pi);ylabel(角度(db)set(gca,xtick,0,0.25,0.4,1); set(gca,ytick,-1,0,1);grid on;subplot(2,2,3);myplot(bz,az,as);%调用自定义绘图函数set(gca,xtick,0,0.25,0.4,1);set(gca,ytick,-50,-15,-1,0); axis(0,1,-40,5);grid on4、运行结果及相应波形n = 5wc = 103.2016（2） 采用双极性变换法设计一个椭圆高通滤波器，要求：通带fp=250hz，rp=1db；阻带fs=150hz，as=20db，滤波器采样频率为fs=100hz。设计步骤：1、用matlab信号处理工具箱函数ellipord，ellip是椭圆滤波器设计函数，该程序中调用以下格式： n,wpo=ellipord(wp,ws,rp,as,s)该格式用于计算椭圆数字滤波器的最低阶数n和通带边界频率wpo，调用参数wp和ws分别为椭圆滤波器的通带边界频率和阻带边界频率的归一化值，rp，as分别为通带最大衰减和阻带最小衰减（db）。 b,a=ellip(n,rp,as,wpo,high,s)该格式用于计算椭圆数字滤波器系统函数分子和分母多项式的字数向量b,a。调用参数n,和wpo分别为椭圆数字滤波器的阶数和通带边界频率wpo，high表示为设计高通滤波器。开始2、滤波器设计流程图读入要求设计指标的fp,fs,as,rp,fs用双极性变换法转换wp,ws为wp1,ws1调用工具箱函数ellipord, ellip求n,wpo,b,a调用bilinear将模拟转为数字滤波器参数sh调用子程序(函数)计算h，w，pha调用子程序(函数)绘出h，pha，损耗曲线结束3、程序实现代码如下：%双极性变换法设计iir椭圆数字高通滤波器close all;clear all;%=rp=1; as=20; fp=250; fs=150; fs=1000;t=1/fs; %=wp=2*pi*fp/fs; ws=2*pi*fs/fs; %设置指标参数wp1=(2/t)*tan(wp/2); ws1=(2/t)*tan(ws/2);%双极性变换法边界频率转换指标h1=10(-rp/20); %通带波纹幅度 rp=-20lg(h1)h2=10(-as/20); %阻带波纹幅度 as=-20lg(h2)n,wpo=ellipord(wp1,ws1,rp,as,s);%计算椭圆高通滤波器最低阶数和通带边界频率b,a =ellip(n,rp,as,wpo,high,s);%计算滤波器系统函数的系数bd,ad=bilinear(b,a,fs);%用双极性变换法转换成数字滤波器h,w =freqz(bd,ad);%计算h(z)的幅频响应subplot(2,2,1);plot(w/(2*pi)*fs,abs(h); xlabel(频率(hz); ylabel(幅度(db) set(gca,xtick,0,fs,fp,fs/2); set(gca,ytick,0,h2,h1,1); text(fp,h1,-通带最大衰减,fontsize,8);%文本注释text(fs,h2,-阻带最小衰减,fontsize,8);%文本注释title(椭圆高通滤波器);axis(0,fs/2,0,1.1);grid on;subplot(2,2,2);pha=angle(h)*180/pi;plot(w/2/pi*fs,pha);grid onset(gca,xtick,0,fs,fp,fs/2); set(gca,ytick,-180,0,180);title(相位响应曲线)subplot(2,2,3);myplot(bd,ad,as);grid onset(gca,xtick,0,0.3,0.5,1);4、运行结果及相应波形n = 3wpo = 2.0000e+003 （3）采用双极性变换法设计一个切比雪夫i型数字高通滤波器，要求：通带wp1=0.3，wp2=0.7，rp=1db，阻带ws1=0.2，ws2=0.8，as=20db，滤波器采样频率fs=1000hz。设计步骤：1、用matlab信号处理工具箱函数cheb1ord，cheby1是切比雪夫i型滤波器设计函数，该程序中调用以下格式： n,wc=cheb1ord(wp,ws,rp,as,s)该格式用于计算切比雪夫i型数字滤波器的最低阶数n和通带边界频率wc， b,a=cheby1(n,rp,wc,s)该格式用于计算切比雪夫i型数字滤波器系统函数分子和分母多项式的字数向量b,a。调用参数n,和wc分别为切比雪夫i型数字滤波器的阶数和通带边界频率wc。开始2、滤波器设计流程图读入要求设计指标的wp1,ws1,wp2,ws2,as,rp,fs用双极性变换法转换wp1,ws1,wp2,ws2为wp1,ws1,wp2,ws2调用工具箱函数cheb1ord, cheby1求n,wc,b,a调用bilinear将模拟转换为数字滤波器参数调用子程序(函数)计算h，w，pha调用子程序(函数)绘出h，pha，损耗曲线结束3、程序实现代码如下：%双极性变换法设计iir切比雪夫i型数字带通滤波器close all;clear all;%=要求设计指标=wp1=0.3*pi; wp2=0.7*pi; ws1=0.2*pi; ws2=0.8*pi; rp=1; as=20; fs=1000;t=1/fs; %=wp1=(2/t)*tan(wp1/2); wp2=(2/t)*tan(wp2/2); ws1=(2/t)*tan(ws1/2); ws2=(2/t)*tan(ws2/2); ws=ws1,ws2;wp=wp1,wp2;%通带宽度h1=10(-rp/20);%通带波纹幅度 rp=-20lg(h1)h2=10(-as/20);%阻带波纹幅度 as=-20lg(h2)n,wc=cheb1ord(wp,ws,rp,as,s) %计算滤波器阶数和3db截止频率wcb,a =cheby1(n,rp,wc,s); %计算滤波器系统函数的系数bz,az=bilinear(b,a,fs); %用双极性变换法转换成数字滤波器h,w=freqz(bz,az); %计算h(z)的幅频响应subplot(2,2,1);plot(w/2/pi*fs,abs(h),k);xlabel(频率(hz); ylabel(幅度(db) ;axis(0,500,0,1.1);set(gca,xtick,0,ws1*fs/(2*pi),wp1*fs/(2*pi),wp2*fs/(2*pi),ws2*fs/(2*pi),fs/2);set(gca,ytick,0,h2,h1,1); title(切比雪夫i型数字带通滤波器);grid onsubplot(2,2,2),plot(w/2/pi*fs,angle(h)/pi*180,k);ylabel(phi);title(相位函数曲线);axis(0,fs/2,-180,180);set(gca,xtick,0,ws1*fs/(2*pi),wp1*fs/(2*pi),wp2*fs/(2*pi),ws2*fs/(2*pi),fs/2);set(gca,ytick,-180,0,180);grid onsubplot(2,2,3);myplot(bz,az,as);grid on;set(gca,xtick,0,ws1/pi,wp1/pi,wp2/pi,ws2/pi);4、运行结果及相应波形n = 3wc = 1.0e+003 * 1.0191 3.9252（4） 采用双极性变换法设计切比雪夫i型数字带阻滤波器，要求：下通带wp1=0.2，上通带wp2=0.8，rp=1db，阻带下限ws1=0.3，阻带上限ws2=0.7，as=20db，滤波器的采样频率fs=1000hz。设计步骤：1、用matlab信号处理工具箱函数cheb1ord，cheby1是切比雪夫i型滤波器设计函数，该程序中调用以下格式： n,wc=cheb1ord(wp,ws,rp,as,s)该格式用于计算切比雪夫i型数字滤波器的最低阶数n和通带边界频率wc， b,a=cheby1(n,rp,wc,stop,s)该格式用于计算切比雪夫i型数字滤波器系统函数分子和分母多项式的字数向量b,a。调用参数n,和wc分别为切比雪夫i型数字滤波器的阶数和通带边界频率wc。stop表示为设计带阻滤波器。2、程序实现代码如下：%双极性变换法设计iir切比雪夫i型数字带阻滤波器close all;clear all;%=要求设计指标=wp1=0.2*pi; wp2=0.8*pi; ws1=0.3*pi; ws2=0.7*pi;rp=1; as=20; t=0.001;fs=1/t; %=wp1=(2/t)*tan(wp1/2); wp2=(2/t)*tan(wp2/2);ws1=(2/t)*tan(ws1/2); ws2=(2/t)*tan(ws2/2); ws=ws1,ws2;wp=wp1,wp2;%通带宽度h1=10(-rp/20);%通带波纹幅度 rp=-20lg(h1)h2=10(-as/20); %阻带波纹幅度 as=-20lg(h2)n,wc=cheb1ord(wp,ws,rp,as,s);%计算滤波器阶数和3db截止频率wcb,a =cheby1(n,rp,wc,stop,s);%计算滤波器系统函数的系数bz,az=bilinear(b,a,fs);%用双极性变换法转换成数字滤波器h,w=freqz(bz,az);%计算h(z)的幅频响应subplot(2,2,1),plot(w/2/pi*fs,abs(h),k);xlabel(频率(hz); ylabel(幅度(db) ;axis(0,fs/2,0,1.1);set(gca,xtick,0,100,150,350,400,500); set(gca,ytick,0,h2,h1,1); title(切比雪夫i型数字带通滤波器);grid onsubplot(2,2,2);plot(w/2/pi*fs,angle(h)/pi*180,k);ylabel(phi);title(相位函数曲线);axis(0,fs/2,-180,180);set(gca,xtick,0,100,150,350,400,500);%获取当前x轴句柄值set(gca,ytick,-180,0,180);grid onsubplot(2,2,3);myplot(bz,az,as);grid on;set(gca,xtick,0,0.2,0.3,0.7,0.8);3、运行结果及相应波形n = 3wc = 1.0e+003 * 0.6498 6.1554题目2、窗函数法设计数字fir滤波器 （1） 用矩形窗设计一个fir数字低通滤波器，要求：n=64,截止频率为wp=0.45，描绘理想和实际滤波器的脉冲相应，窗函数及滤波器的幅频特性响应曲线。设计步骤：1、用matlab信号处理工具箱函数fir1是用于设计线性相位fir数字滤波器的工具箱函数，实现线性相位fir数字滤波器的标准窗函数法设计，调用格式如下： hn=fir1(m,wc,low,boxcar(n)用于计算6db截止频率为wc的m阶单位脉冲响应hn。 理想的单位脉冲响应函数：hd=sin(wc*(n-a)/(pi*(n-a); b=boxcar(n);%返回长度为n的矩形窗函数b开始2、滤波器设计流程图读入窗口的长度计算hd(n)调用窗函数子程序求b计算h(n)=hd(n)*b调用子程序(函数)计算理想和实际的h，w，pha调用子程序(函数)绘出理想和实际的h，pha，损耗曲线结束3、程序实现代码如下：%用矩形窗设计fir数字低通滤波器close all;clear all;%=要求设计指标=n=64;wc=0.4*pi;%=n=0:n-1;a=(n-1)/2;hd=sin(wc*(n-a)/(pi*(n-a);b=boxcar(n);%返回长度为n的矩形窗函数bhn=hd.*(b);%计算理想滤波器的单位脉冲响应h,w=freqz(hn);%计算理想h(z)的幅频响应pha=angle(h);%相频特性subplot(3,2,1);stem(n,hn,.);xlabel(n);ylabel(h(n);title(理想低通滤波器脉冲响应);subplot(3,2,2);hn1=fir1(n-1,wc/pi,low,boxcar(n);%计算hnstem(n,hn1,.);axis(0,70,-0.1,0.4);grid ontitle(实际低通滤波器脉冲响应);xlabel(n);ylabel(hn1);subplot(3,2,3);plot(w,abs(h);grid on ;set(gca,xtick,0,0.4);title(理想幅频响应曲线);xlabel(omega/pi);ylabel(幅值/db);subplot(3,2,4);h1,w1=freqz(hn1);%计算实际h(z)的幅频响应plot(w/pi,abs(h1);axis(0,1,0,1.1);grid on;title(实际幅频响应曲线);xlabel(omega/pi);ylabel(幅值/db);set(gca,xtick,0,0.4); subplot(3,2,5);stem(n,b,.);axis(0,70,0,1.2);title(矩形窗函数曲线);xlabel(n);ylabel(幅度);subplot(3,2,6)a=1;myplot(hn1,a);set(gca,xtick,0,0.4);4、运行结果及相应波形 （2） 用合适的窗函数设计一个fir数字高通滤波器，要求：通带截止频率为wp=0.45，rp=0.5db，阻带截止频率为ws=0.3，as=20db，描绘该滤波器的脉冲响应，窗函数机滤波器的幅频响应曲线和相频响应曲线。设计步骤：1、用matlab信号处理工具箱函数fir1是用于设计线性相位fir数字滤波器的工具箱函数，实现线性相位fir数字滤波器的标准窗函数法设计，调用格式如下： hn=fir1(m,wc,high, bartlett (n)用于计算6db截止频率为wc的m阶单位脉冲响应hn。 b= bartlett (n);%返回长度为n的三角形窗函数b n1=ceil(6.1*pi/bt);%计算三角形窗hn所需长度2、程序实现代码如下：%用三角形窗设计fir数字高通滤波器close all;clear all;%=要求设计指标=wp=0.45*pi;ws=0.3*pi;as=20;%=%由as=20db，查表得知使用三角形窗函数符合设计要求wc=(wp+ws)/2; %计算理想高通滤波器的通带截止频率bt=wp-ws; %过渡带宽度n1=ceil(6.1*pi/bt);%计算三角形窗hn所需长度n=n1+mod(n1+1,2)%如果n为偶数加1，保证n为奇数hn=fir1(n-1,wc/pi,high,bartlett(n);%计算hnsubplot(3,2,1);yn=hn;tstem(hn,yn);%绘制频率响应hn波形title(三角形高通滤波器频率响应函数);grid on;subplot(3,2,2);a=1;myplot(hn,a)%绘制损耗函数波形set(gca,xtick,0,ws/pi,wp/pi); %获取当前x轴句柄值subplot(3,2,3);h,w=freqz(hn);plot(w/pi,abs(h);set(gca,xtick,0,ws/pi,wp/pi);grid on; ylabel(幅值/db);axis(0.25,0.8,0,1.2);title(幅频响应曲线);xlabel(omega/pi);subplot(3,2,4);plot(w/pi,angle(h);set(gca,xtick,0,ws/pi,wp/pi); axis(0,1,-3,3);title(相频响应曲线);grid on;xlabel(omega/pi);ylabel(相位/rad);subplot(3,2,5);wn=bartlett(n);%返回长度为n的三角形窗函数wnn=0:n-1;stem(n,wn,.);axis(0,45,0,1.2);title(三角形窗函数曲线);xlabel(n);ylabel(幅度);grid on;3、运行结果及相应波形n = 41 （3） 选择合适的窗函数设计一个fir数字带通滤波器，要求：下通带截止频率ws1=0.2，as=65db，通带低端截止频率wp1=0.3，rp=0.05db；通带高端截止频率wp2=0.8，as=65db，描绘实际滤波器的脉冲响应，窗函数机滤波器的幅频响应曲线机相频响应曲线。设计步骤：1、用matlab信号处理工具箱函数fir1是用于设计线性相位fir数字滤波器的工具箱函数，实现线性相位fir数字滤波器的标准窗函数法设计，调用格式如下： hn=fir1(m,wc,blackman(n); 用于计算6db截止频率为wc的m阶单位脉冲响应hn。 b=blackman(n);%返回长度为n的布莱克曼窗函数b n1=ceil(11*pi/bt);%计算布莱克曼窗hn所需长度2、程序实现代码如下：%用布莱克曼窗设计fir数字带通滤波器close all;clear all;%=要求设计指标=wp1=0.3*pi;wp2=0.7*pi;ws1=0.2*pi;ws2=0.8*pi;as=65;%=%由as=65db，查表得知使用布莱克曼窗函数符合设计要求wc=(wp1+ws1)/2,(wp2+ws2)/2;bt=wp1-ws1;%过渡带宽度n=ceil(11*pi/bt)%计算布莱克曼窗hn所需长度hn=fir1(n-1,wc/pi,blackman(n);%计算hnsubplot(3,2,1);yn=hn;tstem(hn,yn);%绘制hn波形title(布莱克曼带通滤波器频率响应函数);grid on;subplot(3,2,2);a=1;myplot(hn,a);set(gca,xtick,0,ws1/pi,wp1/pi,wp2/pi,ws2/pi); subplot(3,2,3);h,w=freqz(hn);plot(w/pi,abs(h);axis(0,1,0,1.1);grid onset(gca,xtick,0,ws1/pi,wp1/pi,wp2/pi,ws2/pi); title(幅频响应曲线);xlabel(w);ylabel(幅值/db);subplot(3,2,4);plot(w/pi,angle(h);axis(0,1,-3,3);grid on;set(gca,xtick,0,ws1/pi,wp1/pi,wp2/pi,ws2/pi); xlabel(w);ylabel(相位/rad);title(相频响应曲线);subplot(3,2,5);wn=blackman(n);%返回长度为n的布莱克曼窗函数wnn=0:n-1;stem(n,wn,.);axis(0,120,0,1.2);title(布莱克曼窗函数曲线);xlabel(n);ylabel(幅度);grid on;3、运行结果及相应波形 n = 110 （4） 选择合适的窗函数设计一个fir数字带阻滤波器，要求：下通带截止频率wp1=0.2，rp=0.1db；阻带低端截止频率ws1=0.3，as=40db，阻带高端截止频率ws2=0.7，as=40db，上通带截止频率wp2=0.8，rp=0.1db；描绘实际滤波器的脉冲响应，窗函数机滤波