滤波器设计实例.doc

_%凯泽窗设计示例1：n,Wn,beta,ftype=kaiserord(7 13 17 23,0 1 0,0.01 0.01 0.01,100); %得出滤波器的阶数n=38，beta=3.4fc1=10;fc2=20;fs=100;w1=2*fc1/fs; w2=2*fc2/fs;%将模拟滤波器的技术指标转换为数字滤波器的技术指标window=kaiser(n+1,beta);%使用kaiser窗函数b=fir1(n,w1 w2,window);%使用标准频率响应的加窗设计函数fir1freqz(b,1,512);%数字滤波器频率响应t = (0:100)/fs;s = sin(2*pi*t*5)+sin(2*pi*t*15)+sin(2*pi*t*30);%混和正弦波信号sf = filter(b,1,s);%对信号s进行滤波%凯泽窗设计FIR滤波器高通，带通，低通滤波器设计示例2：clear all;wb=0.317*pi;ws=0.272*pi;wa=0.136*pi;As=65;tr_width=wb-ws; %计算过渡带宽N=(ceil(As-7.95)*2*pi/(14.36*tr_width)+1)/2; %计算凯泽窗的长度disp(窗函数的长度为,num2str(N); %在Matlab命令窗口显示N的数值beta=0.1102*(As-8.7); %计算凯泽窗的线性相位斜率值，即凯泽窗系数n=0:1:N-1;disp(线性相位斜率为,num2str(beta); %在Matlab命令窗口显示线性相位斜率值w_kai=(kaiser(N,beta); %凯泽窗函数hd_L=ideal_lp(wa,N); %调用理想脉冲函数程序,Lowpass，这里调用函数ideal_lph_L=hd_L.*w_kai; %实际脉冲响应=理想脉冲响应*凯泽窗函数db_L,mag,pha,grd,w=freqz_m(h_L,1); %求幅频响应的各个参数，db为幅度响应，w为对应的频率值（弧度表示）%绘图程序figure(4)plot(w/pi,db_L);hold onhd_B=ideal_lp(wb,N)-ideal_lp(wa,N); %调用理想脉冲函数程序,Bandpassh_B=hd_B.*w_kai; db_B,mag,pha,grd,w=freqz_m(h_B,1); plot(w/pi,db_B,m);hold onwc_H=(wb+ws)/2;hd_H=ideal_lp(pi,N)-ideal_lp(wc_H,N); %调用理想脉冲函数程序,Highpassh_H=hd_H.*w_kai; db_H,mag,pha,grd,w=freqz_m(h_H,1); plot(w/pi,db_H,:r);legend(lowpass,bandpass,highpass,1);title(幅度响应/dB);axis(0 1 -120 10);xlabel(以为单位的频率);ylabel(分贝数/dB);%相关函数1functiondb,mag,pha,grd,w=freqz_m(b,a);%db是0到内的相对振幅，mag是0到内的绝对振幅，pha是0到内的相位响应%grd是0到内的群延迟，w是0到内的501个频率样本向量%b是直接型滤波器分子多项式系数，a是直接型滤波器分母多项式系数H,w=freqz(b,a,1000,whole);H=(H(1:1:501);w=(w(1:1:501);mag=abs(H);db=20*log10(mag+eps)/max(mag); %num2str(eps) 的数值是2.22e-16pha=angle(H);grd=grpdelay(b,a,w);%相关函数2function hd=ideal_lp(wc,N);alpha=(N-1)/2;%计算冲击脉冲响应的相位n=0:1:(N-1);m=n-alpha+eps;%加入eps，确保分母是非零数hd=sin(wc*m)./(pi*m);%直接传递函数设计低通滤波器示例3：Fs=50;N=25;tao=(N-1)/2;wc=2*15/50;%n是奇数H(N)为偶数对称，采用N为奇数的I型，取样化为归一化频率m=fix(wc/(2/N)+1);%fix向下取整omega=0:N-1*2*pi/N;ABS_H=ones(1,m),0.4,zeros(1,N-2*m-1),0.4,ones(1,m-1);H=ABS_H.*exp(-j*tao.*omega);h=ifft(H);figure(1);subplot(2,1,1);stem(omega,abs(H);title(理想低通频谱取样点);subplot(2,1,2);stem(0:N-1,real(h);title(频率取样法设计的低通脉冲响应);Hh=fft(h,512);Hh_db=20*log10(abs(Hh);figure(2);omega=0:255*pi/256;plot(omega,Hh_db(1:256);title(滤波器的低通衰耗特性);set(gca,ytick,(-120:20:-20),-3,0);%gca是指最近绘制的那个图像的句柄%ytick表示对Y坐标按照后面的形式显示t=(0:99)/Fs;s=sin(2*pi*t*10)+sin(2*pi*t*20);sf=fft(s,512);y=conv(s,h);%conv多项式的乘法，s*hyf=fft(y,512);figure(3);subplot(3,1,1);plot(0.5*(1:length(y)-1),y); %subplot(3,1,2);w1=(0:255)/256*(50/2);plot(w1,abs(sf(1:256);subplot(3,1,3);plot(w1,abs(yf(1:256);%窗函数设计滤波器示例4：%低通滤波器设计fp=1000;fc=1200;as=100;ap=1;fs=22000;wp=2*fp/fs;wc=2*fc/fs;N=ceil(as-7.95)/(14.36*(wc-wp)/2)+1;beta=0.1102*(as-8.7);window=Kaiser(N+1,beta);b=fir1(N,wc,window);freqz(b,1,512,fs); %高通滤波器：fs=22000;Wp=2*5000/fs;Ws=2*4800/fs;Ap=1;As=100;N=ceil(8*pi/(Wp-Ws)+1;N=N+mod(N+1,2)+1;Wc=(Wp+Ws)/2/pi;h=fir1(N,Wc,high);omega=linspace(0,pi,512);freqz(h,1,omega); 带通滤波器：fs=22000;Wp1=2*1200/fs;Wp2=2*3000/fs;Wc1=2*1000/fs;Wc2=2*3200/fs;Ap=1;As=100;W1=(Wp1+Wc1)/2;W2=(Wp2+Wc2)/2;wdth=min(Wp1-Wc1),(Wc2-Wp2);N=ceil(11*pi/wdth)+1;b = fir1(N,W1 W2);freqz(b,1,512,fs)双线性法设计低通滤波器：fp=1000;fc=1200;as=100;ap=1;fs=22000;wp=2*fp/fs;wc=2*fc/fs; %归一化截止频率n,wn=ellipord(wp,wc,ap,as);%求数字滤波器的最小阶数和归一化截止频率b,a=ellip(n,ap,as,wn);%求传递函数的分子分母系数freqz(b,a,512,fs);%汉明窗设计FIR低通滤波器，处理音频信号示例5：Fp=500;Fs=900;fedge=Fp Fs;Fc=(Fp+Fs)/2;mval=1 0;dev=0.0114469 0.01778279;Ft=11025;N,ftps,mag,wt=remezord(fedge,mval,dev,Ft);fprintf(filter oder is %dn,N)b=remez(N,ftps,mag,wt);figure(1)freqz(b)title(滤波器频率响应)figure(2)x,fs=wavread(北风.wav);y=awgn(x,25);y1=filter(b,1,y);subplot(3,1,1);plot(x);title(原始语音信号)subplot(3,1,2);plot(y);title(加噪语音信号)subplot(3,1,3);plot(y1);title(滤波后的语音信号)wavwrite(y,fs,北风0.wav)wavwrite(y1,fs,北风1.wav)soundsc(y1,fs)figure(3)n=4092;i=0:n-1;x1=fft(x,n);x2=ff