基于MATLAB下数字滤波器设计程序.doc

附录二1.窗函数法的Matlab实现(1) 利用fir1函数及hamming窗设计低通滤波器f1=100;f2=200;%待滤波正弦信号频率fs=2000;%采样频率m=(0.3*f1)/(fs/2);%定义过度带宽M=round(8/m);%定义窗函数的长度N=M-1;%定义滤波器的阶数b=fir1(N,0.5*f2/(fs/2);%使用fir1函数设计滤波器%输入的参数分别是滤波器的阶数和截止频率figure(1)h,f=freqz(b,1,512);%滤波器的幅频特性图%H,W=freqz(B,A,N)当N是一个整数时函数返回N点的频率向量和幅频响应向量plot(f*fs/(2*pi),20*log10(abs(h)%参数分别是频率与幅值xlabel(频率/赫兹);ylabel(增益/分贝);title(滤波器的增益响应);figure(2)subplot(211)t=0:1/fs:0.5;%定义时间范围和步长s=sin(2*pi*f1*t)+sin(2*pi*f2*t);%滤波前信号plot(t,s);%滤波前的信号图像xlabel(时间/秒);ylabel(幅度);title(信号滤波前时域图);subplot(212)Fs=fft(s,512);%将信号变换到频域AFs=abs(Fs);%信号频域图的幅值f=(0:255)*fs/512;%频率采样plot(f,AFs(1:256);%滤波前的信号频域图xlabel(频率/赫兹);ylabel(幅度);title(信号滤波前频域图);figure(3)sf=filter(b,1,s);%使用filter函数对信号进行滤波%参数分别为滤波器系统函数的分子和分母多项式系数向量和待滤波信号输入subplot(211)plot(t,sf)%滤波后的信号图像xlabel(时间/秒);ylabel(幅度);title(信号滤波后时域图);axis(0.2 0.5 -2 2);%限定图像坐标范围subplot(212)Fsf=fft(sf,512);%滤波后的信号频域图AFsf=abs(Fsf);%信号频域图的幅值f=(0:255)*fs/512;%频率采样plot(f,AFsf(1:256)%滤波后的信号频域图xlabel(频率/赫兹);ylabel(幅度);title(信号滤波后频域图);（2）利用fir1函数及Kaiser窗设计带通滤波器Rs=0.01;fs=8000;%采样频率fcuts=1000 1300 2210 2410;a=0 1 0;dev=Rs*ones(1,length(a);M,Wc,beta,ftype=kaiserord(fcuts,a,dev,fs);%M为能够满足设计要求的滤波器的最小阶数，Wc为滤波器的截止频率点%第一个元素f为待设计滤波器的过渡带的起始点和结束点%第二个元素a指定第一个元素频率段的理想幅度值%第三个元素dev中的元素为各通带和阻带内允许的幅度最大误差M=mod(M,2)+M;window=Kaiser(M+1,beta);b=fir1(M,Wc,ftype,window);%输入的第一个参数是滤波器的阶数%第二个参数是滤波器的截止频率%第三个参数是滤波器的类型，stop为带阻滤波器%第四个参数是采用的窗函数h,f=freqz(b,1,512);%滤波器的幅频特性图%H,W=freqz(B,A,N)当N是一个整数时函数返回N点频率向量和幅频响应向量figure(1)plot(f*fs/(2*pi),20*log10(abs(h)%参数分别是频率与幅值xlabel(频率/赫兹);ylabel(增益/分贝);title(滤波器的增益响应);f1=500;f2=1500;f3=2000;f4=3000;%待滤波正弦信号频率t=(0:200)/fs;%定义时间的步长s=sin(2*f1*pi*t)+sin(2*f2*pi*t)+sin(2*f3*pi*t)+sin(2*f4*pi*t);sf=filter(b,1,s);%使用filter函数对信号进行滤波figure(2)subplot(211)plot(t,s);%滤波前的信号图像xlabel(时间/秒);ylabel(幅度);title(信号滤波前时域图);subplot(212)Fs=fft(s,512);AFs=abs(Fs);f=fs/512*(0:255);plot(f,AFs(1:256);%滤波前的信号频域图xlabel(频率/赫兹);ylabel(幅度);title(信号滤波前频域图);figure(3)subplot(211)plot(t,sf)%滤波后的信号图像xlabel(时间/秒);ylabel(幅度);title(信号滤波后时域图);axis(0.005 0.025 -4 4);subplot(212)Fsf=fft(sf,512);%滤波后的信号频域图AFsf=abs(Fsf);%信号频域图的幅值f=(0:255)*fs/512;%频率采样plot(f,AFsf(1:256)%滤波后的信号频域图xlabel(频率/赫兹);ylabel(幅度);title(信号滤波后频域图);（3）利用fir1函数及Kaiser窗设计多通带滤波器Rs=0.01;fs=200;%采样频率fcuts=10 20 40 50 60 70 80 90;a=0,1,0,1,0;dev=Rs*ones(1,length(a);M,Wc,beta,ftype=kaiserord(fcuts,a,dev,fs);M=mod(M,2)+M;window=Kaiser(M+1,beta);b=fir1(M,Wc,ftype,window);h,f=freqz(b,1,512);%滤波器的幅频特性图figure(1)plot(f*fs/(2*pi),20*log10(abs(h)%参数分别是频率与幅值xlabel(频率/赫兹);ylabel(增益/分贝);title(滤波器的增益响应);f1=5;f2=20;f3=30;f4=55;f5=75;f6=95;%待滤波正弦信号频率t=(0:200)/fs;%定义时间的步长s1=sin(2*f1*pi*t)+sin(2*f2*pi*t)+sin(2*f3*pi*t);s=s1+sin(2*f4*pi*t)+sin(2*f5*pi*t)+sin(2*f6*pi*t);%滤波前信号sf=filter(b,1,s);%使用filter函数对信号进行滤波figure(2)subplot(211)plot(t,s);%滤波前的信号图像xlabel(时间/秒);ylabel(幅度);title(信号滤波前时域图);subplot(212)Fs=fft(s,512);AFs=abs(Fs);f=fs/512*(0:255);plot(f,AFs(1:256);%滤波前的信号频域图xlabel(频率/赫兹);ylabel(幅度);title(信号滤波前频域图);figure(3)subplot(211)plot(t,sf)%滤波后的信号图像xlabel(时间/秒);ylabel(幅度);title(信号滤波后时域图);subplot(212)Fsf=fft(sf,512); AFsf=abs(Fsf);%滤波后的信号频域图及信号频域图的幅值f=(0:255)*fs/512;%频率采样plot(f,AFsf(1:256)%滤波后的信号频域图xlabel(频率/赫兹);ylabel(幅度);title(信号滤波后频域图)2.频率取样法的Matlab实现(1) 利用频率取样法设计低通滤波器M=63; Wp=0.5*pi;%所需频率采样点个数及通带截止频率m=0:(M+1)/2; Wm=2*pi*m./(M+1);%通频带上的采样点及阻带截止频率mtr=floor(Wp*(M+1)/(2*pi)+2;%向负方向入floor(3.5)=3;floor(-3.2)=-4Ad=Wm=Wp;Ad(mtr)=0.28;Hd=Ad.*exp(-j*0.5*M*Wm);%构造频域采样向量H(k)Hd=Hd conj(fliplr(Hd(2:M/2+1);%fliplr函数实现矩阵的左右翻转conj是求复数的共轭h=real(ifft(Hd);%h（n）=IDFTH(k)w=linspace(0,pi,1000);%用于产生0,pi之间的1000点行矢量H=freqz(h,1,w);%滤波器的幅频特性图figure(1)plot(w/pi,20*log10(abs(H);%参数分别是归一化频率与幅值xlabel(归一化频率);ylabel(增益/分贝);title(滤波器的增益响应);axis(0 1 -50 0);f1=200;f2=700;f3=800%待滤波正弦信号频率fs=2000;%采样频率figure(2)subplot(211)t=0:1/fs:0.25;%定义时间范围和步长s=sin(2*pi*f1*t)+sin(2*pi*f2*t)+sin(2*pi*f3*t);%滤波前信号plot(t,s);%滤波前的信号图像xlabel(时间/秒);ylabel(幅度);title(信号滤波前时域图);subplot(212)Fs=fft(s,512);%将信号变换到频域AFs=abs(Fs);%信号频域图的幅值f=(0:255)*fs/512;%频率采样plot(f,AFs(1:256);%滤波前的信号频域图xlabel(频率/赫兹);ylabel(幅度);title(信号滤波前频域图);figure(3)sf=filter(h,1,s);%使用filter函数对信号进行滤波subplot(211)plot(t,sf)%滤波后的信号图像xlabel(时间/秒);ylabel(幅度);title(信号滤波后时域图);axis(0.2 0.25 -2 2);%限定图像坐标范围subplot(212)Fsf=fft(sf,512);%滤波后的信号频域图AFsf=abs(Fsf);%信号频域图的幅值f=(0:255)*fs/512;%频率采样plot(f,AFsf(1:256)%滤波后的信号频域图xlabel(频率/赫兹);ylabel(幅度);title(信号滤波后频域图);3.优化设计的Matlab实现(1) 利用Remez函数设计等波纹低通滤波器fs=2000;%设定采样频率rp=3;%通带波纹rs=40;%阻带波纹f=500 600;%截止频率a=1 0;%期望幅度dev=(10(rp/20)-1)/(10(rp/20)+1) 10(-rs/20);n,fo,ao,w=remezord(f,a,dev,fs);%函数remezord返回参数n表示滤波器的阶数%FIR滤波器有B个频带时，f,a,dev分别为2B-2，B,B个元素的向量b=remez(n,fo,ao,w);%函数remez的返回值为n阶FIR滤波器的系数%fo，ao是2B个元素的向量，分别表示B个频带的2B边界频率及幅度值%w是B个元素的向量，表示各频带的加权值figure(1)freqz(b,1,1024,fs);%滤波器的特性图f1=400;f2=700;%待滤波正弦信号频率t=0:1/fs:0.1;%定义时间范围和步长s=sin(2*pi*f1*t)+sin(2*pi*f2*t);%滤波前信号figure(2)subplot(211)plot(t,s);%滤波前的信号图像xlabel(时间/秒);ylabel(幅度);title(信号滤波前时域图);subplot(212)Fs=fft(s,512);%将信号变换到频域AFs=abs(Fs);%信号频域图的幅值f=(0:255)*fs/512;%频率采样plot(f,AFs(1:256);%滤波前的信号频域图xlabel(频率/赫兹);ylabel(幅度);title(信号滤波前频域图);figure(3)sf=filter(b,1,s);%使用filter函数对信号进行滤波subplot(211)plot(t,sf)%滤波后的信号图像xlabel(时间/秒);ylabel(幅度);title(信号滤波后时域图);subplot(212)Fsf=fft(sf,512);%滤波后的信号频域图AFsf=abs(Fsf);%信号频域图的幅值f=(0:255)*fs/512;%频率采样plot(f,AFsf(1:256)%滤波后的信号频域图xlabel(频率/赫兹);ylabel(幅度);title(信号滤波后频域图);(2) 利用Remez函数设计等波纹带通滤波器fs=2000;%设定采样频率rp=3;%通带波纹rs=40;%阻带波纹f=200 300 600 700;%截止频率a=0 1 0;%期望幅度dev=10(-rs/20) (10(rp/20)-1)/(10(rp/20)+1) 10(-rs/20);n,fo,ao,w=remezord(f,a,dev,fs); b=remez(n,fo,ao,w);figure(1)freqz(b,1,1024,fs);%滤波器的特性图f1=100;f2=400;f3=500;f4=800;%待滤波正弦信号频率t=0:1/fs:0.1;%定义时间范围和步长s=sin(2*pi*f1*t)+sin(2*pi*f2*t)+sin(2*pi*f3*t)+sin(2*pi*f4*t);%滤波前信号figure(2)subplot(211)plot(t,s);%滤波前的信号图像xlabel(时间/秒);ylabel(幅度);title(信号滤波前时域图);subplot(212)Fs=fft(s,512);%将信号变换到频域AFs=abs(Fs);%信号频域图的幅值f=(0:255)*fs/512;%频率采样plot(f,AFs(1:256);%滤波前的信号频域图xlabel(频率/赫兹);ylabel(幅度);title(信号滤波前频域图);figure(3)sf=filter(b,1,s);%使用filter函数对信号进行滤波subplot(211)plot(t,sf)%滤波后的信号图像xlabel(时间/秒);ylabel(幅度);title(信号滤波后时域图);subplot(212)Fsf=fft(sf,512);%滤波后的信号频域图AFsf=abs(Fsf);%信号频域图的幅值f=(0:255)*fs/512;%频率采样plot(f,AFsf(1:256)%滤波后的信号频域图xlabel(频率/赫兹);ylabel(幅度);title(信号滤波后频域图);(3) 利用Remez函数设计等波纹带阻滤波器fs=2000;%设定采样频率rp=3;%通带波纹rs=40;%阻带波纹f=200 300 600 700;%截止频率a=1 0 1;%期望幅度dev=10(-rs/20) (10(rp/20)-1)/(10(rp/20)+1) 10(-rs/20);n,fo,ao,w=remezo