信号与系统实验指导书(带源程序)

1、实验一 离散时间系统与MATLAB一. 实验目的1. 进一步加深对离散时间系统的理解。2. 学习在MATLAB中怎样表示离散时间信号。3. 熟悉离散时间信号的作图。二. 实验步骤1. 复习离散时间系统的有关内容。2. 复习MATLAB的基本语法。3. 按实验内容熟悉 stem。4. 编写程序。5. 输出结果，总结结论，按要求写出实验报告。三. 实验内容1掌握stem函数STEM(Y) plots the data sequence Y as stems from the x axis terminated with circles for the data value.STEM(X,Y) pl

2、ots the data sequence Y at the values specified in X.例：t=0:0.1:2; x=cos(pi*t+0.6); stem(t,x); xn=4,2,2,3,6,7; stem(xn);思考：STEM(Y)与STEM(X,Y)有什么不同？STEM与PLOT函数有什么不同？ 2掌握subplot函数H = SUBPLOT(m,n,p), or SUBPLOT(mnp), breaks the Figure window into an m-by-n matrix of small axes, selects the p-th axes for

3、the current plot, and returns the axis handle. The axes are counted along the top row of the Figure window, then the second row, etc. 例：n1=0:3;x1=1,1,1,1;subplot(221);stem(n1,x1);title('x1序列');n2=0:7;x2=1,2,3,4,4,3,2,1;subplot(222);stem(n2,x2);title('x2序列');n3=0:7;x3=4,3,2,1,1,2,3,4;

4、subplot(223);stem(n3,x3);title('x3序列');n4=0:7;x41=cos(pi/4)*n4);subplot(224);stem(n4,x41);title('x4序列');思考：subplot是怎样分配各个作图分区的顺序号的？3信号的运算，请作出，的图形。思考：假如与长度不同，序列的求和和乘积运算能否执行，结果怎样？程序：x1=1,0.7,0.4,0.1,0x2=0.1,0.3,0.5,0.7,0.9subplot(121);stem(x1+x2);title('x1+x2');subplot(122);ste

5、m(x1.*x2);title('x1*x2');4掌握freqs与freqz函数freqs用于s域的频率响应的计算，f reqz用于z域的频率响应的计算。请参看两者的help文件。例：建立2个函数文件，分别输入为系统函数的分子，分母，和角频率的最大值，输出为系统频率响应的db，幅度，相位，和角频率。参看下面的程序。functiondb,mag,pha,Omega=freqs_m(b,a,Omega_Max)% s域频率响应的计算Omega=0:1:500*Omega_Max/500;H=freqs(b,a,Omega);mag=abs(H);db=20*log10(mag+e

6、ps)/max(mag);pha=angle(H);functiondb,mag,pha,w=freqz_m(b,a,w_Max);% z域频率响应的计算w=0:1:499*w_Max/500;H=freqz(b,a,w);mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H);（1）若，请画出该系统幅频和相频特性。程序：b=1a=1 5 6Omega_Max=2*pidb,mag,pha,Omega=freqs_m(b,a,Omega_Max)subplot(221);plot(Omega,mag);title('AF的幅度响应&#

7、39;);subplot(222);plot(Omega,db);title('AF的幅度响应db');subplot(223);plot(Omega,pha);title('AF的相位响应');（2）若，请画出该的幅频和相频特性。程序：b=1 1a=1 0.8 0.64w_Max=2*pidb,mag,pha,w=freqz_m(b,a,w_Max);subplot(221);plot(w,mag);title('AF的幅度响应');subplot(222);plot(w,db);title('AF的幅度响应db');subpl

8、ot(223);plot(w,pha);title('AF的相位响应');思考：freqs，freqz的输入的各个参数与系统函数怎么联系起来，各种调用形式之间有什么不同点？观察上述系统的幅频和相频曲线，能得出什么结论？四. 实验报告要求1. 简述实验目的及原理。2. 概括各函数的常用调用方式，记录程序运行的各种结果。3. 回答思考题。实验二 用FFT作谱分析一. 实验目的1. 进一步加深DFT算法原理和基本性质的理解（FFT是DFT的一种快速算法，FFT的运算结果必然满足DFT的基本性质）。2. 熟悉FFT算法原理和FFT 函数的应用。3. 学习用FFT对连续信号和时域离散信号

9、进行谱分析的方法，了解可能出现的分析误差及其原因。二. 实验步骤1. 复习DFT的定义、性质和用 DFT作谱分析的有关内容。2. 复习FFT算法原理与编程思想。熟悉FFT算法的MATLAB实现。注：MATLAB提供fft函数来计算的DFT，fft函数是用机器语言，而不是以MATLAB指令写成的，因此执行速度快。格式（1）： y=fft(x) 计算x的FFT变换y。当x为矩阵，计算x中每一列信号的离散傅氏变换。当x的长度为2的幂时，采用基2算法，否则采用分裂基算法。格式（2）： y=fft(x,n) 计算x的n点FFT，当x长度大于n时，截断x，否则补零。Plot线性绘图函数。stem：绘制离散

10、序列图。subplot：多坐标设置与定位当前坐标系。figure：创建新的图形窗口（用于输出图形的窗口）。3. 产生下列信号并进行谱分析。（1） （2） （3） （4） （5） （6） 4. 编写程序5. 输出结果，总结结论，按要求写出实验报告。三. 实验内容1. 对上述6个信号，逐个进行谱分析对于， ：。对于 ：在一个周期内取点抽样，做频谱分析。（参考）每个信号谱分析的流程设计为：开始读入长度产生实验信号绘制时间波形图用fft函数进行谱分析曲线绘制结束程序：% fftfigure(1);n1=0:3;x1=1,1,1,1;subplot(221);stem(n1,x1);title('

11、;x1序列');k1=0:7;y11=fft(x1,8);magy11=abs(y11);subplot(222);stem(k1,magy11);title('x1的8点FFT');k2=0:15;y12=fft(x1,16);magy12=abs(y12);subplot(224);stem(k2,magy12);title('x1的16点FFT');%x2figure(2);n2=0:7;x2=1,2,3,4,4,3,2,1;subplot(221);stem(n2,x2);title('x2序列');k1=0:7;y21=fft(

12、x2,8);magy21=abs(y21);subplot(222);stem(k1,magy21);title('x2的8点FFT');k2=0:15;y22=fft(x2,16);magy22=abs(y22);subplot(224);stem(k2,magy22);title('x2的16点FFT');%x3figure(3);n3=0:7;x3=4,3,2,1,1,2,3,4;subplot(221);stem(n3,x3);title('x3序列');k1=0:7;y31=fft(x3,8);magy31=abs(y31);subpl

13、ot(222);stem(k1,magy31);title('x3的8点FFT');k2=0:15;y32=fft(x3,16);magy32=abs(y32);subplot(224);stem(k2,magy32);title('x3的16点FFT');%x4figure(4);n41=0:7;x41=cos(pi/4)*n41);subplot(221);stem(n41,x41);title('x4的8点序列');k1=0:7;y41=fft(x41,8);magy41=abs(y41);subplot(222);stem(k1,magy

14、41);title('x4的8点FFT');n42=0:15;x42=cos(pi/4)*n42);subplot(223);stem(n42,x42);title('x4的16点序列');k2=0:15;y42=fft(x42,16);magy42=abs(y42);subplot(224);stem(k2,magy42);title('x4的16点FFT');%x5figure(5);n51=0:7;x51=sin(pi/8)*n51);subplot(221);stem(n51,x51);title('x5的8点序列');k

15、1=0:7;y51=fft(x51,8);magy51=abs(y51);subplot(222);stem(k1,magy51);title('x5的8点FFT');n52=0:15;x52=sin(pi/8)*n52);subplot(223);stem(n52,x52);title('x5的16点序列');k2=0:15;y52=fft(x52,16);magy52=abs(y52);subplot(224);stem(k2,magy52);title('x5的16点FFT');%x6figure(6);f1=4;f2=8;f3=10;t=

16、0:0.01:1;x6=cos(f1*2*pi*t)+cos(f2*2*pi*t)+cos(f3*2*pi*t);subplot(221);plot(t,x6);title('x6模拟信号');%x61 8T=0.5/8;t=0:T:0.5;x6=cos(f1*2*pi*t)+cos(f2*2*pi*t)+cos(f3*2*pi*t);k61=0:7;y61=fft(x6,8);magy61=abs(y61);subplot(222);stem(k61,magy61);title('x6的8点FFT');%x62 16T=0.5/16;t=0:T:0.5;x6=

17、cos(f1*2*pi*t)+cos(f2*2*pi*t)+cos(f3*2*pi*t);k62=0:15;y62=fft(x6,16);magy62=abs(y62);subplot(224);stem(k62,magy62);title('x6的16点FFT');2. 令，求出其的FFT。当已求出，根据DFT的对称性，求出和，与起第1步求出的、比较。程序：%x7figure(7);n4=0:15;X4=cos(pi/4)*n4);n5=0:15;X5=sin(pi/8)*n5);n7=0:15;x7=x42+x52;subplot(221);stem(n7,x7);titl

18、e('x7的16点序列');k7=n7;y7=fft(x7,16);magy7=abs(y7);subplot(222);stem(k7,magy7);title('x7的16点FFT');X4=real(y7);X5=imag(y7)*i;subplot(223);stem(k7,abs(X4);title('x4的16点序列');subplot(224);stem(k7,abs(X5);title('x5的16点序列');3. 令，求出其的FFT。当已求出，根据DFT的对称性，求出和，与起第1步求出的、比较。程序：%x8fig

19、ure(8);n8=0:15;x8=x42+x52*j;subplot(221);stem(n8,abs(x8);title('x8的16点序列');k8=n8;y8=fft(x8,16);magy8=abs(y8);subplot(222);stem(k8,magy8);title('x8的16点FFT');tempk=mod(16-k8),16);tempy8=y8(tempk+1);convy8=real(tempy8)-imag(tempy8)*i;X4=0.5*(y8+convy8);subplot(223);stem(k8,abs(X4);title

20、('x4的16点序列');X5=1/(2*j)*(y8-convy8);subplot(224);stem(k8,abs(X5);title('x5的16点序列');四. 思考题1. 在时，和的幅频特性会相同吗？为什么？呢？2. 如果周期信号的周期预先不知道，如何使用FFT进行分析。五. 实验报告要求1. 简述实验的原理和目的。2. 输出实验中所用的程序以及所得的时域频域性曲线，与理论结果比较，分析误差产生的原因，以及FFT做谱分析时有关参数的选择方法。3. 总结实验所得的主要结论。4. 简要回答思考题。实验三 用双线性变换法设计IIR 数字滤波器一. 实验目的

21、1. 熟悉用双线性变换法设计IIR数字滤波器的原理与方法。2. 掌握数字滤波器的计算机仿真方法。3. 通过对实际心电图信号的滤波作用，获得数字滤波器的感性知识。二. 实验内容人体心电图信号在测量过程中往往受到工业高频干扰，所以必须经过低通滤波处理后，才能作为判断心脏功能的有用信息。下面给出一实际心电图信号采样序列样本，其中存在高频干扰。在实验中，以作为输入序列，滤除其中的干扰成分。要求：程序设计应包含以下几个部分：（1）绘制原始数据图形；（2）设计巴特沃思低通滤波器并绘制其幅频相应曲线；（3）用设计的滤波器对原始数据进行滤波；（4）绘制滤波后的数据图。三. 实验步骤1. 用双线性变换法设计一个

22、巴特沃思低通IIR数字滤波器。设计指标参数为：在通带内频率低于时，最大衰减小于1dB，在阻带内频率区间上，最小衰减大于15dB。以为间隔，输出模拟滤波器原型和数字滤波器在频率区间上的幅频响应特性曲线。 设。2. 用所设计的数字滤波器对实际心电图信号采样序列进行仿真滤波处理，输出比较滤波前后的心电图信号波形图，观察总结滤波作用。3. （参考）可用到的matlab命令： （1）文件I/O：fopen（打开文件） fscanf（读入文件） fclose（关闭文件） （2）信号处理内部函数：buttap（模拟巴特沃思滤波器） poly（多项式系数） freqs（模拟滤波器频响） f reqz（数字滤波

23、器频响） filter（滤波） （3）作图函数：figure（新开图形窗口） subplot（图形定位） title（图形标题） plot（线性图形） stem（序列图形）函数的具体应用请参见 help文件。4. （参考）程序流程图：开始数字指标转化为模拟指标计算滤波器阶数N 得到模拟滤波器原型的H(s)输出频响曲线用双线性变换法得到H(z)结束输出频响曲线对心电图进行滤波读入心电图数据输出滤波后的图形输出滤波前的图形wp=0.2*pi;wr=0.3*pi;Ap=1;Ar=15; %确定数字滤波器指标参数T=1;Omegap=(2/T)*tan(wp/2);Omegar=(2/T)*tan(w

24、r/2); %模拟滤波器指标cs,ds=afd_butt(Omegap,Omegar, db,mag,pha,Omega=freqs_m(cs,ds,pi); %模拟滤波器的频率Ap,Ar); %得到H(s),其中cs为分子多项式系数,ds为分母多项式系数响应figure;subplot(231);plot(Omega,mag);title('AF的幅度响应'); %模拟滤波器响应subplot(232);plot(Omega,db);title('AF的幅度响应db');subplot(233);plot(Omega,pha);title('AF的相位

25、响应');b,a=bilinear(cs,ds,1/T); %双线性变换法,从模拟滤波器到数字滤波器db,mag,pha,w=freqz_m(b,a,pi); %数字滤波器的频率响应subplot(234);plot(w/T,mag);title('DF的幅度响应'); %数字滤波器的频率响应subplot(235);plot(w/T,db);title('DF的幅度响应db');subplot(236);plot(w/T,pha);title('DF的相位响应');%开始滤波fid=fopen('text1.txt',&

26、#39;r'); if fid=0 error('cannot open the input file!')endinputx,N=fscanf(fid,'%d');figure;n=0:1:(N-1);subplot(211);stem(n,inputx);title('未滤波前的心电图');y1=filter(b,a,inputx); %y1=myfilter(inputx);subplot(212);stem(n,y1);title('滤波后的心电图');fclose(fid);functiondb,mag,pha,

27、Omega=freqs_m(b,a,Omega_Max)% s域频率响应的计算Omega=0:1:500*Omega_Max/500;H=freqs(b,a,Omega);mag=abs(H); %幅度响应db=20*log10(mag+eps)/max(mag);%幅度响应的db值pha=angle(H); %相位响应functiondb,mag,pha,w=freqz_m(b,a,w_Max);% z域频率响应的计算w=0:1:499*w_Max/500;H=freqz(b,a,w);mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H)

28、;function b,a=afd_butt(Omegap,Omegar,Ap,Ar); %求H(s)% b 分子多项式系数% a 分母多项式系数% Omegap 通带截止频率(rad/s) Ap通带衰减(dB)% Omegar 阻带截止频率(rad/s) Ar阻带衰减(dB)krp=(10(0.1*Ar)-1)/(10(0.1*Ap)-1)0.5;lemdarp=Omegar/Omegap;N=ceil(log10(krp)/log10(lemdarp);%计算滤波器阶数Omegac=Omegar/(10(0.1*Ar)-1)(1/2/N);%3db截止频率 通带满足指标，阻带指标有剩余b,

29、a=butter(N,Omegac,'s');四. 思考题用双线性变换法设计数字滤波器过程中，变换公式中的取值，对设计结果有无影响？为何？五. 实验报告要求1. 简述实验目的及原理。2. 由输出的曲线及设计过程简述双线性变换的特点。3. 对比滤波前后的心电图信号波形，说明数字滤波器的滤波过程和作用。4. 简要回答思考题。实验四 用窗函数法设计FIR数字滤波器一. 实验目的1. 掌握用窗函数法设计FIR数字滤波器的原理和方法2. 熟悉线性相位FIR数字滤波器特性3. 了解各种窗函数对滤波特性的影响二. 实验原理与方法设所要求的理想数字滤波器的频率响应为，是与其对应的单位脉冲响应，

30、因此：由于是矩形频率特性，故是无限长的非因果序列。而我们所要设计的是FIR 数字滤波器，其单位脉冲响应必然是有限长的，所以要用有限的来逼近无限长的，最有效的方法是截断，用有限长的窗函数来截取，表示为：这种方法即为窗口设计法。三. 实验内容和步骤1. 实验内容（1）用海明窗设计一线性相位低通 FIR数字滤波器，截止频率。窗口长度N=15、33。要求在两种窗口长度情况下，分别求出，输出相应的幅频特性和相频特性曲线，观察3dB和20dB带宽。总结窗口长度N 对滤波特性的影响。%不同N对滤波器的影响figure;N=15; %窗口长度wc=0.25*pi; %截止频率n=0:1:N-1;hd=idea

31、l_lp(wc,N); %求理想滤波器单位脉冲响应subplot(221);stem(n,hd);title('理想单位脉冲响应');w_ham=(hamming(N)'h=hd.*w_ham; %加窗subplot(223);stem(n,h);title('实际单位脉冲响应');H,w=freqz(h); %计算数字滤波器频响函数mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H);subplot(222);plot(w/pi,mag);title('实际的幅度响应');subpl

32、ot(224);plot(w/pi,pha);title('实际的相位响应');figure;plot(w/pi,db);title('实际的幅度db响应');db=20*log10(mag+eps)/max(mag);figure;N=33;wc=0.25*pi;n=0:1:N-1;hd=ideal_lp(wc,N);subplot(221);stem(n,hd);title('理想单位脉冲响应');w_ham=(hamming(N)'h=hd.*w_ham;subplot(223);stem(n,h);title('实际单位脉冲响应');H,w=freqz(h)mag=abs(H);db=20*log10(mag+eps)/max(mag);pha=angle(H);subplot(222);plot(w/pi,mag);title('实际的幅度响应');subplot(224);plot(w/pi,pha);title('实际的相位响应');figure;plot(w/pi,db);title('实际的幅度db响应');ideal_lp函数：function hd=ideal_l