基于matlab的数字滤波器的设计及仿真

数字滤波器的Matlab设计、仿真研究主要内容前言Matlab下IIR滤波器的设计Matlab下FIR滤波器的设计两个有趣的实验1.声音滤波2.图像滤波数字滤波器的分类按功能分：低通、高通、带通、带阻滤波器按滤波器的网络结构分：IIR滤波器和FIR滤波器它们的函数分别为：N阶IIR滤波器：N-1阶FIR滤波器：滤波器的指标滤波器的指标常在频域给出数字滤波器的频率响应特性曲线：H(ejw)=|H(ejw)|ejө(w)|H(ejw)|为幅频特性函数，表示信号通过该滤波器后各频率的衰减情况。ө(w)为相频特性函数，表示信号通过滤波器后在时间上的延时情况。二、MATLAB环境下IIR数字滤波器的设计典型设计先按照一定规则将给出的数字滤波器的技术指标转换成模拟低通滤波器的技术指标，据此产生模拟滤波器原型，然后把模拟低通滤波器原型转换成数字滤波器。完全设计利用函数直接设计出低通、高通、带通和带阻滤波器，用巴特沃斯滤波器来实现。IIR滤波器设计函数任务目标设计一个滤波器滤除指定正弦信号中的噪声，还原正弦信号。该正弦信号为sin(2*pi*100*t)，噪声信号是sin(2*pi*50*t)和sin(2*pi*150*t)。信号采样频率为2KHz。滤波器的性能指标1.巴特沃斯IIR带通滤波器2.Fs=2000Hz,Wp1=90Hz,Wp2=110Hz,Ws1=80Hz,Ws2=120Hz3.Rp=1，Rs=20用MATLAB进行典型的数字滤波器的设计步骤：按一定规则将给出的数字滤波器的技术指标转换成模拟低通滤波器的技术指标；根据转换后的技术指标使用滤波器阶数选择函数，确定最小阶数N和固有频率Wn。运用最小阶数N产生模拟滤波器原型。运用固有频率Wn把模拟低通滤波器原型转换成模拟低通、高通、带通、带阻滤波器。运用冲击响应不变法或双线性变换法来把模拟滤波器转换成数字滤波器。1.典型法设计clearWp1=90;Wp2=110;Ws1=80;Ws2=120;Fs=2000;Rp=1;Rs=20;Wp=[Wp1,Wp2]*2*pi;Ws=[Ws1,Ws2]*2*pi;Bw=Wp2*2*pi-Wp1*2*pi;Wo=2*pi*sqrt(Wp2*Wp1);[n,Wn]=buttord(Wp,Ws,Rp,Rs,'s');[z,p,k]=buttap(n);[b,a]=zp2tf(z,p,k);[Bb,Ab]=lp2bp(b,a,Wo,Bw);[Bbz,Abz]=impinvar(Bb,Ab,Fs);[h,w]=freqz(Bbz,Abz,512);subplot(2,3,1);plot((w/pi)*Fs/2,20*log10(abs(h)));grid;xlabel('w/pi');ylabel('数字带通幅度(dB)')t1=0:1/2000:0.5;A=0.5*sin(2*pi*50*t1);B=sin(2*pi*100*t1);C=2*sin(2*pi*150*t1);subplot(2,3,2);plot(t1,A);subplot(2,3,3);plot(t1,B);subplot(2,3,4);plot(t1,C);D=A+B+C;subplot(2,3,5);plot(t1,D);E=filter(Bbz,Abz,D);subplot(2,3,6);plot(t1,E);滤波器的幅频和相频曲线2.完全滤波器设计clearcloseallt1=0:1/2000:0.25;A=2*sin(2*pi*100*t1);B=sin(2*pi*50*t1);C=3*sin(2*pi*150*t1);Wp1=90;Wp2=110;Ws1=80;Ws2=120;Fs=2000;Rp=1;Rs=20;Wp=[Wp1,Wp2];Ws=[Ws1,Ws2];[n,Wn]=buttord(Wp/(Fs/2),Ws/(Fs/2),Rp,Rs);[b,a]=butter(n,Wn);[h,w]=freqz(b,a,512);subplot(2,2,1);plot((w/pi)*Fs/2,20*log10(abs(h)));grid;xlabel('Hz');ylabel('数字带通幅度(dB)')D=A+B+C;E=filtfilt(b,a,D);subplot(2,2,2);plot(t1,A);subplot(2,2,3);plot(t1,D);subplot(2,2,4);plot(t1,E);滤波器的幅频和相频曲线3.FDATOOL设计在Simulink下进行设计MATLAB命令行键入simulink，进入Simulink界面；选择SignalProcessingBlockset模块下DSPsources中SineWave来产生信号；选择Sum来叠加三个正弦信号；选择DSPsinks下的VectorScope进行示波；选择Filtering下的DigitalFilterDesign模块进行滤波器的导入。滤波器的仿真三个信号叠加：滤波后的结果：可见FDAtool设计的滤波器滤波的效果很好！二、MATLAB环境下FIR数字滤波器的设计窗函数法等波纹逼近法频率采样法窗函数法用窗函数法设计FIR滤波器时，先根据wc和N求出相应的理想滤波器单位脉冲响应hd(n)；选择合适的窗函数w(n)来截取hd(n)的适当长度(即阶数)，以保证实现要求的阻带衰减；最后得到FIR滤波器单位脉冲响应:h(n)=hd(n)*w(n).窗函数法设计滤波器时阶数选择名称近似过度带宽精确过度带宽最小阻带衰减矩形4π/M1.8π/M21dB巴特利特8π/M6.1π/M25dB汉宁8π/M6.2π/M44dB哈明8π/M6.6π/M51dB布莱克曼12π/M11π/M74dB取kaiser窗时用MATLAB中的kaiserord函数来得到长度M任务设计用窗函数法设计一个FIR带通滤波器低端阻带截止频率：wls=0.2π；低端通带截止频率：wlp=0.35π；高端通带截止频率：whp=0.65π；高端阻带截止频率：whs=0.8π；通带最大衰减：Rp=1dB；阻带最小衰减：Rs=60dB.建模Wc=[wlp/π,whp/π];根据阻带衰减选blackman窗(最小阻带衰减为74dB);其窗口长度由过度带宽B=0.15π决定，blackman窗设计的滤波器过渡带宽为12π/M，故M取80，因为M=N+1，所以N=79.使用工具箱函数fir1的格式b=fir1(N,wc,window)编程。程序编写clearcloseallwls=0.2*pi;wlp=0.35*pi;whp=0.65*pi;B=wlp-wls;N=ceil(12/0.15);wc=[wlp/pi-6/N,whp/pi+6/N];hn=fir1(N-1,wc,blackman(N));n=0:79;subplot(2,1,1);stem(n,hn,'.')[h,w]=freqz(hn,1,256);subplot(2,1,2);plot(w/pi,20*log10(abs(h)));grid;四、两个有趣的实验1、利用MATLAB进行声音滤波：利用MATLAB读取wav格式的声音文件，对信号进行离散傅里叶变换，得到其频谱图，给信号加杂音，然后设计一个数字滤波器将杂音滤除，得到原音，最后比较滤波器的滤波效果。语音读取和加噪clear;closeall;[Y,Fs,NBITS]=WAVREAD('WindowsXP.wav');n=length(Y);f=0:Fs/n:Fs*(n-1)/n;K=0:1/(n-1):1;subplot(2,3,1);plot(K,Y);title('语音信号的时域波形');Y1=fft(Y,n);subplot(2,3,3);plot(f,abs(Y1));title('语音信号的频谱');noise=0.02*sin(2*pi*4000*K);subplot(2,3,2);plot(K(1:50),noise(1:50));title('噪声的时域波形');s=length(noise);noise1=fft(noise,s);subplot(2,3,5);plot(f,abs(noise1));title('噪声的频谱');A=Y'+noise;subplot(2,3,4);plot(K,A);title('含噪语音的时域波形');A1=fft(A,s);subplot(2,3,6);plot(f,abs(A1));title('含噪语音的频谱');使用IIR带阻滤波器滤波fs=4000;fw=200;f0=fs/(Fs/2);fc=fw/(Fs/2);wp=[f0-fcf0+fc];ws=[f0-0.5*fcf0+0.5*fc];Ap=1;Rs=30;[n,wc]=buttord(wp,ws,Ap,Rs);[b,a]=butter(n,wc,'stop');figure,freqz(b,a,512);y=filter(b,a,A);y1=fft(y,s);figure,plot(f,abs(y1));sound(y);滤波后的波形图从上图可以看出，滤波器成功的把噪音滤除，用耳机可以听到清晰的声音。2、利用MATLAB进行图像滤波：clear;closeall;I1=imread('xmu.jpg');I2=rgb2gray(I1);I3=imnoise(I2,'salt');[MN]=size(I3);F=fft2(double(I3));F1=fftshift(F);m1=fix(M/2);n1=fix(N/2);foru=1:Mforv=1:ND=sqrt((u-m1)^2+(v-n1)^2);H(u,v)=1/(1+0.414*(D/50)^8);endendF2=H.*F1;F3=ifftshift(F2);I4=abs(ifft2(F3));subplot(211),imshow(I3,[]);title('含有椒盐噪音的校徽')subplot(212),imshow(I4,[]);title('消除椒盐噪音的校徽')高通滤波clear;closeall;I1=imread('xmu.jpg');I2=rgb2gray(I1);[MN]=size(I2);F=fft2(double(I2));F1=fftshift(F);m1=fix(M/2);n1=fix(N/2);foru=1:Mforv=1:ND=sqrt((u-m1)^2+(v-n1)^2);ifD==0H(u,v)=0;elseH(u,v)=1/(1+0.414*(5/D)^4);endendendF2=H.*F1;F3=ifftshift(F2);I3=abs(ifft2(F3));subplot(211),imshow(I2,[]);title('正常显示的校徽')subplot(212),imshow(I3,[]);title('边缘加强的校徽')THEEND！