数字信号处理matlab实验

数字信号处理实验报告题目冲激响应不变法设计IIR数字滤波器姓名班级学号物理与电子学院2020年6月一．技术指标f1kHz通带内，幅度特性下降小于1dB,频率大于fst=1.5kHz的阻带内，衰减大于15dB,抽样频率为fs=10kHz,尝试采用冲激响应不变法设计一个数字巴特沃兹低通滤波器。二．设计结果1.过程数据：N=6OmegaC=7.0865e+03bz=-0.00000.00070.01050.01670.00420.0001az=1.0000-3.34435.0183-4.21902.0725-0.56000.0647C=[]B=1.8701-0.6294-2.15941.14750.2893-0.4503A=1.0000-0.99180.25441.0000-1.06280.36711.0000-1.28980.6929dbHx=0.920215.00032.幅频特性曲线3.结果分析,达到了设计要求，结合dbHx和幅频特性曲线可知，f=1kHz通带时，幅度特性下降为0.9202dB小于1dB,f=1.5kHz的阻带内，衰减为15.003dB大于15dB.曲线走向可证明符合设计要求。三．源程序主函数：%采用冲激响应不变法设计一个数字巴特沃兹低通滤波器%当fp1000Hz，Rp1dB,当fs1500Hz,AS15dB,fs=10KHzclc;clearallOmegaP=2*pi*1000;OmegaS=2*pi*1500;%巴特沃兹低通滤波器技术指标Rp=1;As=15;Fs=10*10^3;%抽样频率wp=OmegaP/Fs;ws=OmegaS/Fs;%数字频率[N,OmegaC]=buttord(OmegaP,OmegaS,Rp,As,'s')%确定巴特沃兹模拟滤波器阶数和OmegaC[b,a]=butter(N,OmegaC,'s');%确定系统函数分子分母系数[bz,az]=impinvar(b,a,Fs)%冲激响应不变法从AF到DF的变换[C,B,A]=tf2par(bz,az)%直接形式转化为并联形式w0=[wp,ws];%数字临界频率Hx=freqz(bz,az,w0);%检验衰减指标[H,w]=freqz(bz,az);dbHx=-20*log10(abs(Hx)/max(abs(H)))[db,mag,pha,grd,w]=freqz_m(bz,az);plot(w/pi,db);%画图xlabel('\omega/\pi');ylabel('dB');axis([0,0.5,-20,5]);set(gca,'xtickmode','manual','xtick',[0,0.1,0.2,0.3,0.4,0.5]);set(gca,'ytickmode','manual','ytick',[-20,-15,-10,-5,-1]);grid;impinvar函数:function[bz,az]=impinvar(b,a,Fs)[R,P,k]=residue(b,a);%将模拟滤波器系数向量转换成极点，留数和增益T=1/Fs;pz=exp(P*T);%将模拟滤波器极点转换为数字滤波器极点pz[bz,az]=residuez(T*R,pz,k);%求数字滤波器的分子，分母系数向量bz=real(bz');az=real(az');%去除运算产生的微小虚数误差freqz_m函数：function[db,mag,pha,grd,w]=freqz_m(b,a);%freqz子程序的改进[H,w]=freqz(b,a,1000,'whole');H=(H(1:1:501))';w=(w(1:1:501))';%区间内501个频率样本向量mag=abs(H);%绝对振幅db=20*log10((mag+eps)/max(mag));%相对振幅pha=angle(H);%相位响应grd=grpdelay(b,a,w);%群延迟tf2par函数:function[C,B,A]=tf2par(b,a);M=length(b);N=length(a);[r1,p1,C]=residuez(b,a);p=cplxpair(p1,1e-9);I=cplxcomp(p1,p);r=r1(I);K=floor(N/2);B=zeros(K,2);A=zeros(K,3);ifK*2==N;fori=1:2:N-2pi=p(i:i+1,:);ri=r(i:i+1,:);[Bi,Ai]=residuez(ri,pi,[]);B(fix((i+1)/2),:)=real(Bi);A(fix((i+1)/2),:)=real(Ai);end[Bi,Ai]=residuez(r(N-1),p(N-1),[]);B(K,:)=[real(Bi)0];A(K,:)=[real(Ai)0];elsefori=1:2:N-1pi=p(i:i+1,:);ri=r(i:i+1,:);[Bi,Ai]=residuez(ri,pi,[]);B(fix((i+1)/2),:)=real(Bi);A(fix((i+1)/2),:)=real(Ai);endendcplxcomp函数:functionI=cplxcomp(p1,p2)I=[];forj=1:1:length(p2)fori=1:1:length(p1)if(abs(p1(i)-p2(j))0.0001)I=[I,i];endendendI=I';