数字信号处理数字滤波器的设计

1、 数字信号处理1 fir（有限脉冲响应）数字滤波器设计1.1设计题目：用汉宁窗设计线性相位fir低通数字滤波器。1.2设计原理：窗函数设计的基本思想是用firdf（有限长单位脉冲响应数字滤波器）逼近希望滤波特性。设希望逼近的滤波器的频率响应函数为，其脉冲响应用表示。为了计算简便方便，通常选择为具有片段常数特性的理想滤波器。因此是无限长非因果序列，不能直接作为firdf的单位脉冲响应。窗函数设计法就是截取为有限长的一段因果序列，并且合适的窗函数进行加权作为firdf的单位脉冲响应。后面看到，截取的长度和加权窗函数的类型都直接影响逼近精度（滤波器指标）。具体的设计步骤如下： 1） 构造希望逼近的频

2、率响应函数。以低通线性相位firdf设计为例，一般选择为线性相位理想低通滤波器，即， ， 2 ）求出。对 进行ift得到3 ）加窗得到firdf的单位脉冲响应 式中，成为窗函数，其长度为n。1.3窗函数：1.3.1 汉宁窗（hanning） w=hanning（n）1.3.2汉宁窗表达式， k=1,2.n1.3.3汉宁窗w的绘图程序:n=21;n=0:n-1; %数据总数和序列序号figure(1) w=hanning(n); %汉宁窗 posplot=2,2, int2str(i); %指定绘制窗函数的图形位置 stem(n,w); %绘出窗函数 hold on xlabel(n); yla

3、bel(w(n); title(汉宁窗); hold off; grid on；1.3.4汉宁窗的仿真波形： 1.4设计指标：通带截止频 率=0.4 rad，阻带截止频率=0.6 rad，阻带最小衰减=40db，1.5设计过程：1） 选择窗函数，计算窗函数长度n。由阻带衰减40db知道，hanning和hamming窗 函数都满足要求，这两种窗函数过渡带宽度相同，选择hanning窗。hanning窗函数的过渡带宽度为8/n,所以应当满足： 8/nb；解得n8/b=402）确定希望逼近的理想低通滤波器截止频率（该点为过渡中心）处滤波器频响幅度为0.5（等价于衰减6db），所以 =0.53）计算

4、4)加窗得到设计结果 1.6 加窗后h（n）及损耗的绘图源程序：wp=0.4*pi;ws=0.6*pi;b=ws-wp;n=ceil(8*pi/b)wc=(wp+ws)/2/pi;hn=fir1(n-1,wc) m=1024;hk=fft(hn,m);n=0:n-1;subplot(1,2,1);stem(n,hn,.);xlabel(n);ylabel(h(n);grid onk=1:m/2+1;w=2*(0:m/2)/m;subplot(1,2,2);plot(w,20*log10(abs(hk(k);axis(0,1,-100,5);xlabel(w/pi);ylabel(20lg|hg

5、(w)|);1.7 运行结果：n =40hn = columns 1 through 12 -0.0009 -0.0010 0.0013 0.0018 -0.0025 -0.0034 0.0047 0.0062 -0.0081 -0.0104 0.0133 0.0168 -0.0211 -0.0267 0.0341 0.0443 -0.0598 -0.0869 0.1483 0.4503 0.4503 0.1483 -0.0869 -0.0598 0.0443 0.0341 -0.0267 -0.0211 0.0168 0.0133 -0.0104 -0.0081 0.0062 0.0047

6、-0.0034 -0.0025 0.0018 0.0013 -0.0010 -0.00091.8 加窗后h（n）及损耗的仿真波形由图显示：1)由于加窗使得，所以，存在误差。 2）用汉宁窗进行数字滤波器的设计简单，实用，但滤波器的边界频率不易精确控制。所设计的滤波器在阻带边界频率附近的衰减最小，距阻带边界频率越远，衰减越大。因此所设计的滤波器性价比较低。2.iir（无限脉冲响应）模拟滤波器设计2.1 设计题目：设计一个切比雪夫i型模拟滤波器2.2设计原理：模拟滤波器的设计是通过查表或计算确定滤波器的系统函数，再进行安装调试，一般的设计步骤如下：1） 根据信号处理要求确定设计指标2） 选择滤波器的

7、类型3） 计算滤波器的阶数n4） 通过查表或计算确定滤波器的系统函数5） 综合实现并安装调试2.3技术指标如下： 通带边界频率wp=0.2*pi，通带最小衰减rp=1阻带截止频率wr=0.4*pi，阻带最小衰减rs=50db2.4 设计源程序：% 主程序：wp=0.2*pi;wr=0.4*pi;rp=1;rs=50; b,a=afd_chb1(wp,wr,rp,rs) c,b,a=sdir2cas(b,a)db,mag,pha,w=freqs_m(b,a,pi);ha,x,t=impulse(b,a);subplot(2,2,1);plot(w/pi,mag);title(幅度响应|ha(jo

8、mega)|);subplot(2,2,2);plot(w/pi,db);title(幅度响应(db);subplot(2,2,3);plot(w/pi,pha/pi);title(相位响应);axis(0,1,-1,1);subplot(2,2,4);plot(t,ha);title(单位脉冲响应ha(t);axis(0,max(t),min(ha),max(ha);% 切比雪夫i型模拟滤波器的设计主程序functionb,a=afd_chb1(wp,wr,rp,rs);ep=sqrt(10(rp/10)-1);a=10(rs/20);omegac=rp;omegar=wr/wp;g=sqr

9、t(a*a-1)/ep;n=ceil(log10(g+sqrt(g*g-1)/log10(omegar+sqrt(omegar*omegar-1);fprintf(n=%2.0fn,n);b,a=u_chblap(n,rp,omegac);% 设计归一化切比雪夫i型模拟滤波器原型子程序functionb,a=u_chblap(n,rp,omegac);z,p,k=cheb1ap(n,rp);a=real(poly(p);ann=a(n+1);p=p*omegac;a=real(poly(p);anu=a(n+1);k=k*anu/ann;b0=k;b=real(poly(z);b=k*b;%系

10、统函数的幅度响应和相位响应子程序function db,mag,pha,w=freqs_m(b,a,wmax);w=0:1:500*wmax/500;h=freqs(b,a,w);mag=abs(h);db=20*log10(mag+eps)/max(mag);pha=angle(h);%直接形式转化成级联形式子程序 function c,b,a=sdir2cas(b,a);na=length(a)-1;nb=length(b)-1;b0=b(1);b=b/b0;a0=a(1);a=a/a0;c=b0/a0;p=cplxpair(roots(a);k=floor(na/2);if k*2=na

11、 a=zeros(k,3); for n=1:2:na arow=p(n:1:n+1,:);arow=poly(arow); a(fix(n+1)/2),:)=real(arow); endelseif na=1 a=0 real(poly(p);else a=zeros(k+1,3); for n=1:2:2*k arow=p(n:1:n+1,:);arow=ploy(arow); a(fix(n+1)/2),:)=real(arow); end a(k+1,:)=0 real(poly(p(na);endz=cplxpair(roots(b);k=floor(nb/2);if nb=0 b

12、=0 0 poly(z);elseif k*2=nb b=zeros(k,3); for n=1:2:nb brow=z(n:1:n+1,:);brow=poly(brow); b(fix(n+1)/2),:)=real(brow); endelseif nb=1 b=0 real(poly(z); for n=1:2:2*k brow=z(n:1:n+1,:);brow=poly(brow); b(fix(n+1)/2),:)=real(brow); end b=0 real(poly(z(nb);end2.5 运行结果:n= 6b = 0.0614a = 1.0000 0.9283 1.9

13、308 1.2021 0.9393 0.3071 0.0689c = 0.0614b = 0 0 1a = 1.0000 0.4641 0.1247 1.0000 0.3398 0.5577 1.0000 0.1244 0.99072.6仿真波形：3 iir（无限脉冲响应）数字滤波器设计3.1 设计题目：用脉冲响应不变法设计低通滤波器 3.2设计指标：通带截止频率wp=0.1*pi，通带衰减ap=1db，阻带截止频率wr=0.6*pi，阻带衰减ar=15。 3.3设计原理：先设计与该数字滤波器的指标相同的模拟滤波器，再通过s平面映射到z平面。一般的设计步骤如下：1）确定数字低通滤波器的技术指标

14、：通带截止频率wp，通带衰减ap，阻带截止频率wf，阻带衰减ar=15db2）将数字低通滤波器的技术指标换成模拟信号的技术指标。3）按照模拟低通滤波器的指标设计数字低通滤波器的指标。4）将模拟低通滤波器的，从s平面映射到z平面，得到低通数字滤波器的系统函数h(z)。3.4 设计源程序：%主程序wp=0.1*pi;wr=0.6*pi;ap=1;ar=15;t=1;omegap=wp/t;omegar=wr/t;cs,ds=afd_chb1(omegap,omegar,ap,ar)c,b,a=sdir2cas(cs,ds);db,mag,pha,omega=freqs_m(cs,ds,pi);su

15、bplot(3,2,3);plot(omega/pi,mag);title(模拟滤波器幅度响应|ha(jomega)|);b,a=imp_invr(cs,ds,t);h,n=impz(b,a);c,b,a=dir2par(b,a)db,mag,pha,grd,w=freqz_m(b,a);subplot(3,2,1);plot(w/pi,mag);title(数字滤波器幅度响应|ha(jomega)|);subplot(3,2,4);plot(w/pi,db);title(数字滤波器幅度响应（db）);subplot(3,2,2);plot(w/pi,pha/pi);title(数字滤波器相位

16、响应);%脉冲响应不变法程序function b,a=imp_invr(c,d,t)r,p,k=residuez(c,d);p=exp(p*t);b,a=residuez(r,p,k);b=real(b).*t;a=real(a);%数字滤波器响应子程序function db,mag,pha,grd,w=freqz_m(b,a);h,w=freqz(b,a,1000,whole);h=(h(1:501);w=(w(1:501);mag=abs(h);db=20*log10(mag+eps)/max(mag);pha=angle(h);grd=grpdelay(b,a,w);%直接型转换成并联型

17、子程序function c,b,a=dir2par(b,a);m=length(b);n=length(a);r1,p1,c=residuez(b,a);p=cplxpair(p1,10000000*eps);i=cplxcomp(p1,p);r=r1(i);k=floor(n/2);b=zeros(k,2);a=zeros(k,3);if k*2=n; for i=1:2:(n-2) brow=r(i:1:(i+1),:); arow=p(i:1:(i+1),:); brow,arow=residuez(brow,arow,); b(fix(i+1)/2),:)=real(brow); a(

18、fix(i+1)/2),:)=real(arow); end brow,arow=residuez(r(n-1),p(n-1),); b(k,:)=real(brow) 0;a(k,:)=real(arow) 0;else for i=1:2:(n-1) brow=r(i:1:(i+1),:); arow=p(i:1:(i+1),:); brow,arow=residuez(brow,arow,); b(fix(i+1)/2),:)=real(brow); a(fix(i+1)/2),:)=real(arow); endend%比较两个含同样标量元素但（可能）有不同下标的复数对及其相应留数向量

19、子程序function i=cplxcomp(p1,p2);i=;for j=1:length(p2) for i=1:length(p1) if(abs(p1(i)-p2(j)0.0001) i=i,i; end endend3.5 运行结果：n= 2cs = 0.0970ds =1.0000 0.3449 0.1088c =b = 0.0970 -0.0923a = 1.0000 -1.6171 0.70833.6 仿真波形：4利用fft对一信号进行频谱分析 4.1设计指标：x(n)=0.1*sin(2*pi*f1*n/fs)+0.5*cos(2*pi*f2*n/fs)+0.5*sin(2

20、*pi*f3*n/fs)中含有三种频率成分，f1=2hz，f2=2.05hz，f3=1.9hz，采样频率为fs=10hz，序列：分别取n1=64,n2=128点有效数据作频谱特性分析，分别在四个图形窗口绘出x(n)， x(k)64点dft，x(k)补零到128点dft，x(k)128点dft。比较得出在哪种情况下可以清楚地分辨出信号的频谱成分。4.2 设计源程序如下： f1=2;f2=2.05;f3=1.9;fs=10;n=64;n=0:n-1; xn=0.1*sin(2*pi*f1*n/fs)+0.5*cos(2*pi*f2*n/fs)+0.5*sin(2*pi*f3*n/fs);k=0:n

21、-1;wn=exp(-j*2*pi/n); %旋转因子nk=n*k; %产生一个台nk值的n乘n维矩阵wnnk=wn.nk; %dft矩阵xk=xn*wnnk; %dft系数的行subplot(2,2,1); %将多个图形画到一个平面plot(n,xn); % 绘制函数曲线 title(a)输入序列);subplot(2,2,2);plot(k,abs(xk); %绘制幅频特性图title(b)x(k)64点dft);n=128;n=0:n-1;xn=0.1*sin(2*pi*f1*n/fs)+0.5*cos(2*pi*f2*n/fs)+0.5*sin(2*pi*f3*n/fs);k=0:n-1;wn=exp(-j*2*pi/n);nk=n*k;wnnk=wn.nk;xk=xn*wnnk;subplot(2,2,3);plot(1:128,abs(xk);title(c)x(k)补零到128点dft);subplot(2,2,4);plot(k,abs(xk);title(d)x(k)128点dft);4.3运行结果：4.4仿真