MATLAB实现数字FIR的高通_和带通等滤波器的源程序

1、利用汉宁窗设计I型数字高通滤波器clear all;Wp=0.6*pi;Ws=0.4*pi; tr_width=Wp-Ws;N=ceil(6.2*pi/tr_width) n=0:1:N-1;Wc=(Ws+Wp)/2; hd=ideal_hp1(Wc,N); w_han=(hanning(N)' h=hd.*w_han;% 过渡带宽度% 滤波器长度% 理想低通滤波器的截止频率% 理想低通滤波器的单位冲激响应% 汉宁窗% 截取得到实际的单位脉冲响应db,mag,pha,w=freqz_m2(h,1);%计算实际滤波器的幅度响应delta_w=2*pi/1000;Ap=-(min(db(W

2、p/delta_w+1:1:501)%实际通带纹波As=-round(max(db(1:1:Ws/delta_w+1) % 实际阻带纹波 subplot(221)stem(n,hd)title(' 理想单位脉冲响应subplot(222)hd(n)')stem(n,w_han) title(' 汉宁窗 w(n)') subplot(223)stem(n,h)title(' 实际单位脉冲响应subplot(224)hd(n)')plot(w/pi,db) title(' 幅度响应 (dB)') axis(0,1,-100,10) c

3、lear all;Wp=0.6*pi;Ws=0.4*pi; tr_width=Wp-Ws;N=ceil(6.2*pi/tr_width) n=0:1:N-1;Wc=(Ws+Wp)/2; hd=ideal_hp1(Wc,N); w_han=(hanning(N)' h=hd.*w_han;% 过渡带宽度% 滤波器长度% 理想低通滤波器的截止频率% 理想低通滤波器的单位冲激响应% 汉宁窗% 截取得到实际的单位脉冲响应db,mag,pha,w=freqz_m2(h,1);%计算实际滤波器的幅度响应delta_w=2*pi/1000;Ap=-(min(db(Wp/delta_w+1:1:501

4、)%实际通带纹波As=-round(max(db(1:1:Ws/delta_w+1) % 实际阻带纹波subplot(221)stem(n,hd)title(' 理想单位脉冲响应 hd(n)') subplot(222)stem(n,w_han)title(' 汉宁窗 w(n)') subplot(223) stem(n,h)title(' 实际单位脉冲响应 hd(n)')subplot(224)plot(w/pi,db)title(' 幅度响应 (dB)')axis(0,1,-100,10)基于切比雪夫一致逼近法设计 FIR 数

5、字低通滤波器给定频率轴分点给定在这些频率分点上理想的幅频响应 给定在这些频率分点上的加权 设计出切比雪夫最佳一致逼近滤波器clear all;f=0 0.6 0.7 1;%A=1 1 0 0;%weigh=1 10;%b=remez(32,f,A,weigh); % h,w=freqz(b,1,256,1); h=abs(h);h=20*log10(h); subplot(211) stem(b,'.'); grid;title(' 切比雪夫逼近滤波器的抽样值 ') subplot(212) plot(w,h);grid;title(' 滤波器幅频特性

6、(dB)')利用汉宁窗设计I型数字带阻滤波器clear all;Wpl=0.2*pi;Wph=0.8*pi;Wsl=0.4*pi;Wsh=0.6*pi;tr_width=min(Wsl-Wpl),(Wph-Wsh);N=ceil(6.2*pi/tr_width)n=0:1:N-1;Wcl=(Wsl+Wpl)/2;Wch=(Wsh+Wph)/2;hd=ideal_bs(Wcl,Wch,N);% 过渡带宽度% 滤波器长度%理想低通滤波器的截止频率% 理想低通滤波器的单位冲激响应w_hann=(hanning(N)'h=hd.*w_hann;db,mag,pha,w=freqz_m2

7、(h,1);过渡带宽度滤波器长度理想低通滤波器的截止理想低通滤波器三角窗截取得到实际的单计算实际滤波器实际通带纹波实际阻带纹波% 汉宁窗%截取得到实际的单位脉冲响应% 计算实际滤波器的幅度响应delta_w=2*pi/1000;Ap=-(min(db(1:1:Wpl/delta_w+1) %实际通带纹波As=-round(max(db(Wsl/delta_w+1:1:Wsh/delta_w+1)%实际阻带纹波 subplot(221)stem(n,hd)title(' 理想单位脉冲响应 hd(n)') subplot(222)stem(n,w_hann) title('

8、 汉宁窗 w(n)')subplot(223) stem(n,h) title(' 实际单位脉冲响应 hd(n)') subplot(224) plot(w/pi,db) title(' 幅度响应 (dB)') axis(0,1,-100,10)利用三角窗设计川型数字带通滤波器clear all;Wpl=0.4*pi;Wph=0.6*pi;Wsl=0.2*pi;Wsh=0.8*pi;tr_width=min(Wpl-Wsl),(Wsh-Wph); %N=ceil(6.1*pi/tr_width) %n=0:1:N-1;Wcl=(Wsl+Wpl)/2; %

9、 频率Wch=(Wsh+Wph)/2;hd=ideal_bp2(Wcl,Wch,N); % 的单位冲激响应w_tri=(triang(N)' %h=hd.*w_tri; % 位脉冲响应db,mag,pha,w=freqz_m2(h,1); % 的幅度响应delta_w=2*pi/1000;Ap=-(min(db(Wpl/delta_w+1:1:Wph/delta_w+1) %As=-round(max(db(Wsh/delta_w+1:1:501) % subplot(221)stem(n,hd)title(' 理想单位脉冲响应 hd(n)') subplot(222)

10、 stem(n,w_tri)title(' 三角窗 w(n)') subplot(223) stem(n,h)title(' 实际单位脉冲响应 hd(n)') subplot(224) plot(w/pi,db)title(' 幅度响应 (dB)')axis(0,1,-100,10)利用布拉克曼窗设计口型数字带通滤波器clear all;Wpl=0.4*pi;Wph=0.6*pi;Wsl=0.2*pi;Wsh=0.8*pi;tr_width=min(Wpl-Wsl),(Wsh-Wph);%过渡带宽度N=ceil(11*pi/tr_width)+1

11、%滤波器长度n=0:1:N-1;Wcl=(Wsl+Wpl)/2; %理想低通滤波器的截止频率Wch=(Wsh+Wph)/2;hd=ideal_bp1(Wcl,Wch,N);%理想低通滤波器的单位冲激响应w_bman=(blackman(N)'%布拉克曼窗h=hd.*w_bman; %截取得到实际的单位脉冲响应db,mag,pha,w=freqz_m2(h,1);%计算实际滤波器的幅度响应delta_w=2*pi/1000;Ap=-(min(db(Wpl/delta_w+1:1:Wph/delta_w+1) % 实际通带纹波 As=-round(max(db(Wsh/delta_w+1:

12、1:501) % 实际阻带纹波 subplot(221)stem(n,hd)title(' 理想单位脉冲响应 hd(n)') subplot(222)stem(n,w_bman)title(' 布拉克曼窗 w(n)') subplot(223) stem(n,h)title(' 实际单位脉冲响应 hd(n)') subplot(224) plot(w/pi,db)title(' 幅度响应 (dB)')axis(0,1,-100,10)利用海明窗设计口型数字低通滤波器clear all;Wp=0.2*pi;Ws=0.4*pi; tr

13、_width=Ws-Wp;N=ceil(6.6*pi/tr_width)+1 n=0:1:N-1;Wc=(Ws+Wp)/2;hd=ideal_lp1(Wc,N); w_ham=(hamming(N)' h=hd.*w_ham;%过渡带宽度%滤波器长度%理想低通滤波器的截止频率% 理想低通滤波器的单位冲激响应 % 海明窗db,mag,pha,w=freqz_m2(h,1);%计算实际滤波器的幅度响应%截取得到实际的单位脉冲响应delta_w=2*pi/1000;Ap=-(min(db(1:1:Wp/delta_w+1)%实际通带纹波As=-round(max(db(Ws/delta_w+

14、1:1:501)%实际阻带纹波subplot(221)stem(n,hd)title(' 理想单位脉冲响应 hd(n)')subplot(222)stem(n,w_ham)title(' 海明窗 w(n)')subplot(223) stem(n,h)title(' 实际单位脉冲响应 hd(n)') subplot(224) plot(w/pi,db) title(' 幅度响应 (dB)') axis(0,1,-100,10)% functiondb,mag,pha,w=freqz_m2(b,a) %滤波器的幅值响应 (相对、绝对

15、 )、相位响应%db:相对幅值响应%mag:绝对幅值响应%pha: 相位响应%w 采样频率；%b系统函数H(z)的分子项(对FIR，b=h)%a系统函数H(z)的分母项(对FIR，a=1)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); % 相对幅值响应 pha=angle(H);%相位响应利用模拟 Butterworth 滤波器设计数字低通滤波器% exa4-8_pulseDF for example4

16、-8%using Butterworth analog lowpass filter to design digital lowpass filter%利用模拟 Butterworth 滤波器设计数字低通滤波器%脉冲响应不变法wp=0.2*pi;ws=0.3*pi;Rp=1;As=15;T=1;%性能指标Rip=10(Rp/20);At n=10(As/20);OmgP=wp*T;OmgS=ws*T;N,OmgC=buttord(OmgP,OmgS,Rp,As,'s');% 选取模拟滤波器的阶数 cs,ds=butter(N,OmgC,'s'); %设计出所需的

17、模拟低通滤波器b,a=impinvar(cs,ds,T); %应用脉冲响应不变法进行转换%求得相对、绝对频响及相位、群迟延响应db,mag,pha,grd,w=freqz_m(b,a);%下面绘出各条曲线subplot(2,2,1);plot(w/pi,mag);title(' 幅频特性 ');xlabel('w(/pi)');ylabel('|H(jw)|');axis(0,1,0,1.1);set(gca,'XTickMode','manual','XTick',0 0.2 0.3 0.5 1)

18、; set(gca,'YTickMode','manual','YTick',0 Atn Rip 1);gridsubplot(2,2,2);plot(w/pi,db);title(' 幅频特性 (dB)');xlabel('w(/pi)');ylabel('dB');axis(0,1,-40,5); set(gca,'XTickMode','manual','XTick',0 0.2 0.3 0.5 1); set(gca,'YTickMod

19、e','manual','YTick',-40 -As -Rp 0); grid subplot(2,2,3);plot(w/pi,pha/pi); title(' 相频特性 ');xlabel('w(/pi)'); ylabel('pha(/pi)');axis(0,1,-1,1); set(gca,'XTickMode','manual','XTick',0 0.2 0.3 0.5 1); gridsubplot(2,2,4); plot(w/pi,grd

20、); title('群延迟 ');xlabel('w(/pi)'); ylabel('Sample'); axis(0,1,0,12);set(gca,'XTickMode','manual','XTick',0 0.2 0.3 0.5 1); gridfunctiondb,mag,pha,grd,w=freqz_m(b,a) %滤波器幅值响应 (绝对、相对 ) 、相位响应及群延迟 %Usage: db,mag,pha,grd,w=freqz_m(b,a) %500 点对应 0,pi%db 相对幅值

21、响应； mag 绝对幅值响应； pha 相位响应； grd 群延迟响应 %w采样频率；b系统函数H(z)的分子项(对FIR, b= h)%a系统函数H(z)的分母项(对FIR, a= 1)H,w=freqz(b,a,500);%500 点的复频响应 mag=abs(H);db=20*log10(mag+eps)/max(mag); pha=angle(H);grd=grpdelay(b,a,w);基于频域抽样法的 FIR 数字带阻滤波器设计clear all;N=41;T1=0.598;alpha= (N-1)/2;l=0:N-1;wl= (2*pi/N)*l; Hrs=ones(1,6),T

22、1,zeros(1,7),T1,ones(1,11),T1,zeros(1,7),T1,ones(1,6);%理想振幅采样响应 Hdr=1,1,0,0,1,1;wdl=0,0.3,0.3,0.7,0.7,1;k1=0:floor(N-1)/2); k2=floor(N-1)/2)+1:N-1; angH=pi/2-alpha*(2*pi)/N*(k1+0.5),-pi/2+alpha*(2*pi)/N*(N-k2-0.5);%相位约束条件 subplot(221) plot(wl/pi+1/N,Hrs,'.',wdl,Hdr) title(' 频率样本 Hd(k) :N

23、=41') axis(0 1 -0.1 1.2) subplot(222) stem(l,h)Hdk=Hrs.*exp(j*angH);h1=ifft(Hdk,N);n=0:1:N-1;h=real(h1.*exp(j*pi*n/N);db,mag,pha,w=freqz_m2(h,1);Hr,ww,a,L=hr_type3(h);%构成 Hd(k)%实际单位冲激响应%实际振幅响应title(' 实际单位脉冲响应 h(n)') subplot(223)plot(ww/pi,Hr,wl/pi+1/N,Hrs,'.') title(' 实际振幅响应

24、H(w)') axis(0 1 -0.1 1.2) subplot(224) plot(w/pi,db) title(' 幅度响应 (dB)') axis(0 1 -80 10) function db,mag,pha,w = freqz_m(b,a); %滤波器的幅值响应 (相对、绝对 )、相位响应 %db:相对幅值响应%mag:绝对幅值响应 %pha: 相位响应 %w 采样频率；%b系统函数H(z)的分子项(对FIR , b=h) %a系统函数H(z)的分母项(对FIR , a=1) H,w = freqz(b,a,1000,'whole');H =

25、 (H(1:1:501)' w = (w(1:1:501)' mag = abs(H);db = 20*log10(mag+eps)/max(mag); pha = angle(H);% pha = unwrap(angle(H); function Hr,w,c,L=hr_type3(h);%计算所设计的 3 型滤波器的振幅响应%Hr=振幅响应%b=3 型滤波器的系数%L=Hr 的阶次%h=3 型滤波器的单位冲击响应 M=length(h);L=(M-1)/2;c= 2*h(L+1:-1:1); n=0:1:L; w=0:1:500'*2*pi/500;Hr=sin(

26、w*n)*c'基于频域抽样法的 FIR 数字带通滤波器设计wsl=0.12*pi;% 低阻带边缘 wsh=0.82*pi;% 高阻带边缘 wpl=0.32*pi;% 低通带边缘 wph=0.62*pi;% 高通带边缘 delta=(wpl-wsl);% 过度带 M=ceil(2*pi*3/delta);% 抽样点数 al=(M-1)/2;wl=(2*pi/M); % 抽样间隔k=0:M-1;T1=0.12; T2=0.6;% 过渡带样本点Hrs=zeros(1,ceil(0.12*pi/wl)+1),T2,T1,ones(1,ceil(0.3*pi/wl),T1,T2,zeros(1,

27、ceil(0.3734*pi/wl) ,T2,T1,ones(1,ceil(0.3*pi/wl),T1,T2,zeros(1,ceil(0.12*pi/wl)+1);wdl=0 0.12 0.32 0.62 0.82 1; k1=0:floor(M-1)/2);k2=floor(M-1)/2)+1:M-1; angH=-al*(2*pi)/M*k1,al*(2*pi)/M*(M-k2);H=Hrs.*exp(j*angH);h=real(ifft(H);% 傅立叶反变换 figure(1);% 冲击响应图 stem(k,h);title('impulse response');

28、 xlabel('n');ylabel('h(n)'); grid;figure(2);% 幅频曲线图 Hf=abs(H);w=k*wl/pi;plot(w,Hf,'*b-')axis(0 1 -0.1 1.1);title('amplitude response');xlabel('frequency in pi units');ylabel('Hr(w)'); set(gca,'xtickmode','manual','xtick',wdl); s

29、et(gca,'ytickmode','manual','ytick',0 0.12 0.6 1); grid;figure(3);fs=15000;c,f3=freqz(h,1);f3=f3/pi*fs/2; plot(f3,20*log10(abs(c); title(' 频谱特性 ');xlabel('频率 /HZ');ylabel('衰减 /dB');grid; t=(0:100)/fs; x=sin(2*pi*t*700)+sin(2*pi*t*3200)+sin(2*pi*t*6200)

30、; q=filter(h,1,x);a,f1=freqz(x); f1=f1/pi*fs/2;b,f2=freqz(q); f2=f2/pi*fs/2;figure(4); subplot(2,1,1);plot(f1,abs(a);title(' 输入波形频谱图 ');xlabel('频率');ylabel('幅度')subplot(2,1,2);plot(f2,abs(b);title(' 输出波形频谱图 ');xlabel('频率');ylabel('幅度')基于汉宁窗的 FIR 数字高通滤波

31、器设计 function s2Fs=15000;t=(0:100)/Fs;x=sin(2*pi*500*t)+sin(2*pi*3000*t)subplot(245);stem(x);title(' 原始信号 ');axis(0,100,-2,2);Ws=7*pi/30;Wp=13*pi/30;tr_wid=Wp-Ws;%过渡带宽度N=ceil(11*pi/tr_wid)%滤波器长度n=0:1:N-1;Wc=(Ws+Wp)/2;%理想高通滤波器的截止频率hd=ideal_hp1(Wc,N);%理想高通滤波器的单位冲激响应w_bla=(blackman(N)'%布拉克曼h

32、=hd.*w_bla;%截取得到实际的单位脉冲响应db,mag,pha,grd,w=freqz_m(h,1); %计算实际滤波器的幅度响应delta_w=2*pi/1000;As=-round(max(db(1:1:Ws/delta_w+1) % 数y=filter(h,1,x) subplot(246) plot(y) title(' 滤波后的信号 '); axis(0,100,-1,1) subplot(241) stem(n,hd) title(' 理想单位脉冲响应 hd(n)') subplot(242) stem(n,w_bla)title('

33、 布拉克满窗 w(n)') subplot(243) stem(n,h)title(' 实际单位脉冲响应 hd(n)') subplot(244) plot(w/pi,db)title(' 幅度响应 (dB)') axis(0,1,-100,10) function db,mag,pha,grd,w = freqz_m(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); pha = angle(H);grd = grpdelay(b,a,w); subplot(247);plot(pha) title(' 相频响应 ') function hd=ideal_hp1(Wc,N) alp=(N-1)/2;n=0:1:N-1; m=n-alp+