数字信号处理实验报告-7

数字信号处理实验报告专业班级电信1101姓名曾文学号20111186020指导老师吴莉华中科技大学武昌分校2014年5月21日实验一信号、系统及系统响应实验目的(1)加深对离散线性移不变(LSI)系统基本理论的理解，明确差分方程与系统函数之间的关系。(2)初步了解用MATLAB语言进行离散时间系统研究的基本方法。(3)掌握求解离散时间系统冲激响应和阶跃响应程序的编写方法，了解常用子函数。(4)通过实验进一步理解卷积定理，了解卷积的过程。(5)了解MATLAB中有关卷积的子函数及其应用方法。实验原理(1)离散LSI系统的响应与激励由离散时间系统的时域和频域分析方法可知，一个线性移不变离散系统可以用线性常系数差分方程表示：系统函数H(z)反映了系统响应与激励的关系。一旦上式中的bm和ak的数据确定了，则系统的性质也就确定了。其中特别注意：a0必须进行归一化处理，即a0＝1。对于复杂信号激励下的线性系统，可以将激励信号在时域中分解为单位脉冲序列或单位阶跃序列，把这些单元激励信号分别加于系统求其响应，然后把这些响应叠加，即可得到复杂信号加于系统的零状态响应。因此，求解系统的冲激响应和阶跃响应尤为重要。由图1－1可以看出一个离散LSI系统响应与激励的关系。同时，图1－1显示了系统时域分析方法和z变换域分析法的关系。如果已知系统的冲激响应h(n)，则对它进行z变换即可求得系统函数H(z)；反之，知道了系统函数H(z)，对其进行z逆变换，即可求得系统的冲激响应h(n)。y(n)=x(n)*h(n)x(n)h(n)y(n)=x(n)*h(n)x(n)h(n)H(z)H(z)Y(z)=X(z)H(z)Y(z)=X(z)H(z)X(z)X(z)图1－1离散LSI系统响应与激励的关系图1－1离散LSI系统响应与激励的关系(2)离散LSI系统的线性卷积由理论学习我们已知，对于线性移不变离散系统，任意的输入信号x(n)可以用及其位移的线性组合来表示，即当输入为 时，系统的输出y(n)=h(n)，由系统的线性移不变性质可以得到系统对x(n)的响应y(n)为称为离散系统的线性卷积，简记为也就是说，如果已知系统的冲激响应，将输入信号与系统的冲激响应进行卷积运算，即可求得系统的响应。实验内容(1)已知一个因果系统的差分方程为满足初始条件y(-1)=0,x(-1)=0,求系统的冲激响应和阶跃响应。编写仿真程序，并调试得到结果，进行分析。(2)已知两个信号序列：f1=0.8n(0<n<20)f2=u(n)(0<n<10)求两个序列的卷积和。编写仿真程序，并调试得到结果，进行分析。程序：①a=[1,0,1/3,0];%分母多项式系数b=[1/6,1/2,1/2,1/6];%分子多项式系数N=32;n=0:N-1;%一维数组，含N个分量，间隔为1hn=impz(b,a,n);%单位冲激响应gn=dstep(b,a,n);%单位阶跃响应subplot(1,2,1),stem(n,hn,'k');title('系统的单位冲激响应');ylabel('h(n)');xlabel('n');subplot(1,2,2),stem(n,gn,'k');title('系统的单位阶跃响应');ylabel('g(n)');xlabel('n');②nf1=0:20;f1=0.8.^nf1;subplot(2,2,1);stem(nf1,f1,'filled');title('f1(n)');nf2=0:10;lf2=length(nf2);f2=ones(1,lf2);subplot(2,2,2);stem(nf2,f2,'filled');title('f2(n)');y=conv(f1,f2);subplot(2,1,2);stem(y,'filled');title('y(n)');。实验二频域采样1.实验目的：(1)掌握频率域采样会引起时域周期化的概念，以及频率域采样定理及其对频域采样点数选择的指导作用。(2)会用MATLAB语言进行频域抽样与恢复时程序的编写方法。实验原理：了解频域采样定理的要点，掌握采样理论的结论：“频域采样时域信号周期延拓”。实验内容：（1）频域采样理论的验证。给定信号如下：（2）编写程序分别对频谱函数在区间上等间隔采样32和16点，得到，再分别对进行32点和16点IFFT，得到。（3）分别画出、的幅度谱，并绘图显示x(n)、的波形，进行对比和分析，验证总结频域采样理论。程序：①M=27;N=32;n=0:M-1;%产生M长三角波序列x(n)xa=1:ceil(M/2);%生成序列n+1xb=floor(M/2):-1:1;%生成序列27-n；floor是向下取整，ceil是向上取整xn=[xa,xb];Xk=fft(xn,1024);%1024点FFT[x(n)],用于近似序列x(n)的FTX32k=fft(xn,32);%32点FFT[x(n)]x32n=ifft(X32k);%32点IFFT[X32(k)]得到x32(n)X16k=X32k(1:2:N);%隔点抽取X32k得到X16(K)x16n=ifft(X16k,N/2);%16点IFFT[X16(k)]得到x16(n)subplot(3,2,2);stem(n,xn,'.');%stem画离散序列图boxon%给图形加边框title('(b)三角波序列x(n)');xlabel('n');ylabel('x(n)');axis([0,32,0,20])k=0:1023;wk=2*k/1024;%产生1024点DFT对应的采样点频率(关于π归一化值)subplot(3,2,1);%画子图，3表示行数，2表示列数，1表示当前子图的序号数(以行元素优先顺序排列)。plot(wk,abs(Xk));%绘制1024点DFT的幅频特性图title('(a)FT[x(n)]');xlabel('\omega/\pi');%'\'是转义符号等价于w/piylabel('|X(e^j^\omega)|');%等价于e的jw次方axis([0,1,0,200])%确定轴的范围，横轴从0~1，纵轴从0~到200k=0:N/2-1;subplot(3,2,3);stem(k,abs(X16k),'.');boxontitle('(c)16点频域采样');xlabel('k');ylabel('|X_1_6(k)|');%下划线表示下标(转义字符‘_’)axis([0,8,0,200])n1=0:N/2-1;subplot(3,2,4);stem(n1,x16n,'.');boxontitle('(d)16点IDFT[X_1_6(k)]');xlabel('n');ylabel('x_1_6(n)');axis([0,32,0,20])k=0:N-1;subplot(3,2,5);stem(k,abs(X32k),'.');boxontitle('(e)32点频域采样');xlabel('k');ylabel('|X_3_2(k)|');axis([0,16,0,200])n1=0:N-1;subplot(3,2,6);stem(n1,x32n,'.');boxontitle('(f)32点IDFT[X_3_2(k)]');xlabel('n');ylabel('x_3_2(n)');axis([0,32,0,20])②wp=0.25*pi;%滤波器的通带截止频率ws=0.4*pi;%滤波器的阻带截止频率Rp=1;As=15;%输入滤波器的通阻带衰减指标ripple=10^(-Rp/20);Attn=10^(-As/20);Fs=100;T=1/Fs;Omgp=(2/T)*tan(wp/2);%原型通带频率预修正Omgs=(2/T)*tan(ws/2);%原型阻带频率预修正[n,Omgc]=buttord(Omgp,Omgs,Rp,As,'s');%计算阶数n和3dB截止频率[z0,p0,k0]=buttap(n);%归一化原型设计ba=k0*real(poly(z0));%求原型滤波器系数baa=real(poly(p0));%求原型滤波器系数a[ba1,aa1]=lp2lp(ba,aa,Omgc);%变换为模拟低通滤波器系数b,a[bd,ad]=bilinear(ba1,aa1,Fs)%用双线性变换法求数字滤波器系数b,a[sos,g]=tf2sos(bd,ad);%由直接型变换为级联型[H,w]=freqz(bd,ad);dbH=20*log10((abs(H)+eps)/max(abs(H)));%化为分贝值subplot(2,2,1),plot(w/pi,abs(H));ylabel('|H|');title('幅度响应');axis([0,1,0,1.1]);set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);gridsubplot(2,2,2),plot(w/pi,angle(H)/pi);ylabel('\phi');title('相位响应');axis([0,1,-1,1]);set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);set(gca,'YTickMode','manual','YTick',[-1,0,1]);gridsubplot(2,2,3),plot(w/pi,dbH);title('幅度响应(dB)');ylabel('dB');xlabel('频率(\pi)');axis([0,1,-40,5]);set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);set(gca,'YTickMode','manual','YTick',[-50,-15,-1,0]);gridsubplot(2,2,4),zplane(bd,ad);axis([-1.1,1.1,-1.1,1.1]);title('零极图');思考题：如果序列x(n)的长度为M，希望得到其频谱在上的N点等间隔采样，当N<M时，如何用一次最少点数的DFT得到该频谱采样？答：先对原序列x(n)以N为周期进行周期延拓后取主值区序列，再计算N点DFT则得到N点频域采样：实验三信号的快速傅里叶变换实验1、实验目的：使学生进一步了解快速付里叶变换的理论；及用Matlab软件设计快速付里叶变换的方法，使学生进一步了解数字信号的频谱概念。2、实验内容与步骤：1）根据数字信号序列的长度决定2的整数次幂N，用Matlab软件设计出快速付里叶变换及反变换的程序，或按照范例程序进行修改，输入信号序列，运行程序，检查程序是否有错。2）添加绘图语句，画出数字信号变换前后的波形。运行程序，观察变换前后的波形有什么不同。3）添加绘图语句，画出变换后的频谱图，运行程序，仔细观察频谱图，了解变换的意义。3.实验程序及现象m=10;fori=1:m+1//系统函数中B为分子多项式的系数，用矩阵的形式表示ifi==1B(i)=1;elseifi==m+1B(i)=-1;elseB(i)=0;endendA=[1,-1];//系统函数中A为分母多项式的系数N=8192;[H,f]=freqz(B,A,N);plot(f*25/pi,abs(H));grid;figure,plot(f*25/pi,angle(H));grid;figure,zplane(B,A);k=0:N-1;f=2*k/N;load('RawData.mat');x=rawdata(1,1:N);w=filter(B,A,x);X=abs(fft(x,N));W=abs(fft(w,N));figure;plot(x);title('输入信号');figure;plot(w);title('输出信号');figure;plot(f,abs(X));title('输入信号的幅频响应');figure;plot(f,abs(W));title('输出信号的幅频响应');figure;plot(f,abs(X),'b',f,abs(W),'r');实验现象如下图图形如下：（3）极点分布图，一共有9个极点（4）输入信号与（5）输出信号的不同，是由于输入信号经过了低通滤波器的滤波，滤除了不必要的波形。（6）输入频响与（7）输出频响的区别，也与低通滤波器有关快速傅里叶变换的目的：快速傅里叶变化是为了为应用各种信号的实时处理提供条件的方法，大大的提高了傅里叶变化的计算量。实验四用双线性变换法设计IIR数字滤波器1.实验目的(1)熟悉用双线性变换法设计IIR数字滤波器的原理与方法。(2)掌握用双线性变换法设计数字滤波器的计算机仿真方法。(3)了解MATLAB有关双线性变换法的子函数。2.实验涉及的MATLAB子函数bilinear功能：双线性变换——将s域映射到z域的标准方法，将模拟滤波器变换成离散等效滤波器。调用格式：[numd,dend]=bilinear(num,den,Fs)；将模拟域传递函数变换为数字域传递函数，Fs为取样频率。3.实验原理(1)双线性变换法是将整个s平面映射到整个z平面，其映射关系为双线性变换法克服了脉冲响应不变法从s平面到z平面的多值映射的缺点，消除了频谱混叠现象。但其在变换过程中产生了非线性的畸变，在设计IIR数字滤波器的过程中需要进行一定的预修正。(2)双线性变换法设计IIR数字滤波器的步骤：①输入给定的数字滤波器设计指标；②根据公式 进行预修正，将数字滤波器指标转换成模拟滤波器设计指标；③确定模拟滤波器的最小阶数和截止频率；④计算模拟低通原型滤波器的系统传递函数；⑤利用模拟域频率变换法，求解实际模拟滤波器的系统传递函数；⑥用双线性变换法将模拟滤波器转换为数字滤波器。 4.实验内容采用双线性变换法设计一个巴特沃斯数字低通滤波器，要求：滤波器采样频率思考题用双线性变换法设计数字滤波器过程中，变换公式中T的取值，对设计结果有无影响?为什么?答：没有，一般取2/T=1，方便计算实验程序现象程序如下：wp=0.25*pi;%滤波器的通带截止频率ws=0.4*pi;%滤波器的阻带截止频率Rp=1;As=15;%输入滤波器的通阻带衰减指标ripple=10^(-Rp/20);Attn=10^(-As/20);Fs=100;T=1/Fs;Omgp=(2/T)*tan(wp/2);%原型通带频率预修正Omgs=(2/T)*tan(ws/2);%原型阻带频率预修正[n,Omgc]=butter(Omgp,Omgs,Rp,As,'s');%计算阶数n和3dB截止频率[z0,p0,k0]=buttap(n);%归一化原型设计ba=k0*real(poly(z0));%求原型滤波器系数baa=real(poly(p0));%求原型滤波器系数a[ba1,aa1]=lp2lp(ba,aa,Omgc);%变换为模拟低通滤波器系数b,a[bd,ad]=bilinear(ba1,aa1,Fs)%用双线性变换法求数字滤波器系数b,a[sos,g]=tf2sos(bd,ad);%由直接型变换为级联型[H,w]=freqz(bd,ad);dbH=20*log10((abs(H)+eps)/max(abs(H)));%化为分贝值subplot(2,2,1),plot(w/pi,abs(H));ylabel('|H|');title('幅度响应');axis([0,1,0,1.1]);set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);set(gca,'YTickMode','manual','YTick',[0,Attn,ripple,1]);gridsubplot(2,2,2),plot(w/pi,angle(H)/pi);ylabel('\phi');title('相位响应');axis([0,1,-1,1]);set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);set(gca,'YTickMode','manual','YTick',[-1,0,1]);gridsubplot(2,2,3),plot(w/pi,dbH);title('幅度响应(dB)');ylabel('dB');xlabel('频率(\pi)');axis([0,1,-40,5]);set(gca,'XTickMode','manual','XTick',[0,0.25,0.4,1]);set(gca,'YTickMode','manual','YTick',[-50,-15,-1,0]);gridsubplot(2,2,4),zplane(bd,ad);axis([-1.1,1.1,-1.1,