数字信号处理matlab仿真

数字信号处理作业设计报告一、目的1.增进对Matlab的认识，加深对数字信号处理理论方面的理解。2.掌握数字信号处理中IIR和FIR滤波器的设计。3.了解和掌握用Matlab实现IIR和FIR滤波器的设计方法、课程，为以后的设计打下良好根底。二、设计内容1.IIR〔无限脉冲响应〕模拟滤波器设计〔1〕设计题目：椭圆型模拟带通IIR滤波器技术指标：通带下截止频率fpl=2kHz，上截止频率fph=5kHz，通带内最大衰减ap=1dB；阻带下截止频率fsl=1.5kHz上截止频率fsh=5.5kHz,阻带最小衰减as=40dB.设计原理：①确定模拟带通滤波器的技术指标，并对边界频率做归一化处理；②确定归一化低通技术要求；③设计归一化低通G(p)；④将低通G(p)转换成带通H(s)。Matlab原程序如下：clearall;fp=[2000,5000];ap=1;fs=[1500,5500];as=40;wp=2*pi*fp;ws=2*pi*fs;%归一化的截止频率[N,wn]=ellipord(wp,ws,ap,as,'s');%求椭圆形滤波器的最小阶数和归一化截止频率[B,A]=ellip(N,ap,as,wn,'s');%求传递函数的分子分母系数[H,w]=freqs(B,A);%频率响应函数f=0:8000;w=2*pi*f;H=freqs(B,A,w);%求系统在指定频率点w上的频响Hplot(f,20*log10(abs(H)));%绘图显示axis([07000-800])仿真波形图如下：〔2〕设计题目：巴特沃斯低通模拟滤波器技术指标：通带截止频率fp=5kHz，通带内最大衰减ap=2dB；阻带截止频率fs=12kHz，阻带最小衰减as=30dB。设计原理：①确定模拟带通滤波器的技术指标，并对边界频率做归一化处理；②确定归一化低通技术要求并求出归一化低通原型系统函数Ga(p);③将Ga(p)去归一化。Matlab原程序如下：clearall;wp=2*pi*5000;ws=2*pi*12000;Rp=2;As=30;%设置滤波器参数[N,wc]=buttord(wp,ws,Rp,As,'s');%计算滤波器阶数N和3dB截止频率wc[B,A]=butter(N,wc,'s');%计算滤波器系统函数分子分母多项式系数k=0:511;fk=0:14000/512:14000;wk=2*pi*fk;Hk=freqs(B,A,wk);subplot(2,2,1);plot(fk/1000,20*log10(abs(Hk)));gridonxlabel('频率(kHz)');ylabel('幅度(dB)')axis([0,14,-40,5])仿真波形图如下：2.IIR〔无限脉冲响应〕数字滤波器设计〔1〕设计题目：IIR数字带通滤波器设计技术指标：采用巴特沃斯滤波器,一个数字系统的抽样频率Fs=2000Hz，设计一个为此系统使用的带通数字滤波器Hdbp(z)。要求通带范围为300~400Hz，在带边频率处的衰减不大于3dB；在200Hz以下和500Hz以上衰减不小于18dB。设计原理：对于数字带通滤波器的设计，通常所用方法为双线性变换。可借助于模拟滤波器的频率变换设计一个所需类型的过渡模拟滤波器，再通过双线性变换将其转换成所需类型的数字滤波器。Matlab原程序如下：clearall;fp=[300400];fs=[200500];rp=3;rs=18;Fs=2000;wp=fp*2*pi/Fs;ws=fs*2*pi/Fs;%%Firstlytofinishfrequencyprewarping;wap=2*Fs*tan(wp./2);was=2*Fs*tan(ws./2);[n,wn]=buttord(wap,was,rp,rs,'s');%Note:'s'![z,p,k]=buttap(n);[bp,ap]=zp2tf(z,p,k);%bw=wap(2)-wap(1);w0=sqrt(wap(1)*wap(2));[bs,as]=lp2bp(bp,ap,w0,bw);%[h1,w1]=freqs(bp,ap);figure(1);plot(w1,abs(h1));grid;ylabel('lowpassG(p)')%w2=[0:Fs/2-1]*2*pi;h2=freqs(bs,as,w2);%Note:z=(2/Ts)(z-1)/(z+1);[bz1,az1]=bilinear(bs,as,Fs);[h3,w3]=freqz(bz1,az1,1000,Fs);figure(2);plot(w2/2/pi,20*log10(abs(h2)),w3,20*log10(abs(h3)));grid;ylabel('BandpassAFandDF')xlabel('Hz');仿真波形图如下：图①模拟低通原型图图②模拟和数字带通滤波器的幅频响应〔2〕设计题目：基于Matlab函数直接设计切比雪夫Ⅰ型数字低通IIR滤波器技术指标：通带截止频率fp为100HZ，阻带截止频率fs为200HZ，通带最大衰减ap为3dB，阻带最小衰减as为30dB,采样频率Fs为1000HZ.设计原理：Matlab工具箱中提供了直接设计数字滤波器的函数使设计更加简单。切比雪夫Ⅰ型滤波器在通带是等波纹的，而在阻带是单调下降的，利用cheby1函数MATLAB语法可以直接设计；而切比雪夫Ⅰ型滤波器得阶次可以通过MATLAB中的cheb1ord函数确定。Matlab原程序如下：clearall;fp=100; ap=3;%通带衰减fs=200;as=30;%阻带衰减Fs=1000;[N,Wn]=cheb1ord(fp/(Fs/2),fs/(Fs/2),ap,as);%计算ChebyshevⅠ型滤波器的最小阶数和截止频率[b,a]=cheby1(N,ap,Wn);freqz(b,a,512,1000);%采用1000HZ的采样频率绘出显示滤波器的幅频和相频响应仿真波形图如下：FIR〔有限脉冲响应〕数字滤波器设计〔1〕设计题目：基于矩形窗和汉明窗的FIR滤波器设计技术指标：分别用矩形窗和汉明窗，设计一低通滤波器，令N=10，频率响应截止频率Hd(ejω)在0≤ω≤ππ≤ω≤π之间为0。设计原理：要设计FIR滤波器，h(n)必然是有限长的，所以要用有限长的h(n)来逼近无限长的hd(n)，最有效的方法是截断hd(n)，即用一个有限长度的窗函数ω(n)来截取hd(n)，并将截短后的hd(n)移位，得由于ω(n)的长度为N，所以h(n)是因果的。因此窗函数序列ω(n)的形状及长度的选择就很关键。Matlab原程序如下：clearall;N=10;b1=fir1(N,0.25,boxcar(N+1));%用矩形窗作为冲激响应的窗函数b2=fir1(N,0.25,hamming(N+1));%用汉明窗作为冲激响应的窗函数%M=128;h1=freqz(b1,1,M);h2=freqz(b2,1,M);%分别求两个滤波器的频率响应t=0:10subplot(221)stem(t,b2,'.');holdon;plot(t,zeros(1,11));grid;f=0:0.5/M:0.5-0.5/M;M1=M/4;fork=1:M1hd(k)=1;hd(k+M1)=0;hd(k+2*M1)=0;hd(k+3*M1)=0;endsubplot(222)plot(f,abs(h1),'b-',f,abs(h2),'g-',f,hd,'-');grid;仿真波形图如下：〔2〕设计题目：海明窗函数法FIR数字滤波器技术指标：ππ，阻带最小衰减为As＝50dB。设计原理：π/N,因此具有较小的阶数。在MATLAB中利用ideal_lp函数，计算理想低通滤波器的单位脉冲响应hd〔n〕。functionhd=ideal_lp(Wc,N);alpha=(N-1)/2;n=0:1:N-1;m=n-alpha+eps;hd=sin(Wc*m)./(pi*m);利用freqz_m函数计算绝对和相对的幅度响应。function[db,mag,pha,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);以上两个函数应先分别保存在两个新建m文件中存放在Matlab软件工作路径下。Matlab原程序如下：clearall;Wp=0.2*pi;Ws=0.4*pi;%滤波器的边界频率tr_width=Ws-Wp;%过渡带宽N=ceil(6.6*pi/tr_width)+1%根据过度带宽等于海明窗函数主瓣宽度求得滤波器所用海明窗的最小长度n=0:1:N-1;wc=(Ws+Wp)/2;%截止频率在通带和阻带边界频率的中点hd=ideal_lp(wc,N);%调用理想低通滤波器脉冲响应函数w_ham=(hamming(N))';%采用海明窗h=hd.*w_ham;%在时间域乘积对应于频率域的卷积，截取得到实际单位脉冲响应[db,mag,pha,w]=freqz_m(h,[1]);%调用freqz_m函数计算实际滤波器绝对和相对的幅度响应delta_w=2*pi/1000;Ap=-(min(db(1:1:Wp/delta_w+1)));%实际通带最大衰减As=-round(max(db(Ws/delta_w+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('实际单脉冲响应h(n)')subplot(224);plot(w/pi,db);title('幅度响应(db)')axis([0,1,-100,10]);运行结果：N=34仿真波形图如下：4、利用FFT进行频谱分析设一序列中含有三种频率成分，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。比拟得出在哪种情况下可以清楚地分辨出信号的频谱成分设计要求：因为f1=2Hz,f2=2.05Hz,f3=1.9Hz采样频率为fs=10Hz,所以构造正弦序列x(n)=sin(2nπf1/fs)+sin(2nπf2/fs)+sin(2nπf3/fs).设计原理：信号采样的频谱分辨率为F=0.05Hz,所以采样记录最小点数为N=2fh/F=2f2/F=82点。利用Matlab离散傅里叶变换函数对序列进行频谱变换，通过Matlab频谱波形图分析信号频谱成分Matlab原程序如下：clearall;N1=64;fs=10;n=0:1:N1-1;f1=2;f2=2.05;f3=1.9;xn=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs)+sin(2*pi*f3*n/fs);%满足条件的正弦序列x(n)Xk=fft(xn);AXk=abs(Xk(1:N1/2));%对x(n)进行64点DFT运算figure(1)subplot(2,1,1)plot(n,xn);%绘出64点采样正弦序列x(n)figure(2)k=(0:N1/2-1)*fs/N1;subplot(2,1,2)plot(k,AXk);%绘出X(k)64点DFT频谱幅度%X(k)补零到128点DFTN2=128;xn=[xnzeros(1,N2-N1)];%对64点采样序列x(n)进行补零到128点Xk=fft(xn);AXk=abs(Xk(1:N2/2));%对补零到128点的序列进行DFT运算m=0:1:N2-1;figure(3)k=(0:N2/2-1)*fs/N2;subplot(2,1,2)plot(k,AXk)%绘出X(k)补零到128点DFT频谱幅度%X(k)128点DFTn=0:1:N2-1;xn=sin(2*pi*f1*n/fs)+sin(2*pi*f2*n/fs)+sin(2*pi*f3*n/fs);%满足条件的正弦序列Xk=fft(xn);AXk=abs(Xk(1:N2/2));%对x(n)进行128点DFT运算figure(4)k=(0:N2/2-1)*fs/N2;subplot(2,1,2)plot(k,AXk);%绘出X(k)的128点DFT频谱幅度结果与仿真波形图如下：图①64点采样正弦序列x(n)图②X(k)64点DFT图③X(k)补零到128点DFT图④X(k)128点DFT因为采样最小点数为82点，所以当N1=64点采样时进过D