基于MATLAB的图形仿真与矫正.doc

武汉理工大学学生实验报告书实验课程名称：测试信号分析与处理开 课 学 院 ：机 电 工 程 学 院指导教师姓名：陈 祯 学 生 姓 名 ：乾 引 学生专业班级：测 控1303 学号班级序号：0121304940718 测试信号分析与处理实验实验一 离散系统的时频域分析一、 实验目的 通过该实验熟悉matlab软件的基本操作指令，掌握 matlab软件的使用方法，掌握离散系统的时频域分析的基本原理、方法以及matlab函数的调用。二、 实验设备1、 微型计算机1台；2、 matlab软件1套三、 实验原理 matlab 软件是由mathworks公司于1984年推出的一套科学计算软件，分为总包和若干个工具箱，其中包含用于信号分析与处理的sptool工具箱和用于滤波器设计的fdatool工具箱。它具有强大的矩阵计算和数据可视化能力，是广泛应用于 信号分析与处理中的功能强大且使用简单方便的成熟软件。matlab软件中已有大量的关于数字信号处理的运算函数可供调用，本实验主要是针对离散系统的时频域分析的的matlab函数的熟悉和应用。（备注：可上网查询sptool、fdatool工具箱的使用方法，如下网址是使用sptool的ppt介绍。/link?url=g9lkzhnedxcojj6in61zjfm6fvprpmijuxa0wjiybqrxrdfymy227utt_4qxqti6qalbqmphuwqdkcvhb7xfx9lunqqdnxokbw19pljtpde） 相关matlab函数1. 数字信号的表示 (1)已知数字信号采样值的数字信号的表示 x = 0.2 0.4 0.3 0.5 0.9 1.2 0.4; (2)已知数字信号函数表达式的数字信号的表示 t = tmin:1/fs:tmax; x = sin(2*t)-1; 其中tmin-起始时间,tmax-结束时间, fs-采样频率(3)常用数字函数单位脉冲函数delta = imp_fun(n,n0)%n为采样点数，n0为延迟采样周期数阶跃函数u = step_fun(n,n0),exp函数x = 2*exp(-0.5*n);% sinandcos函数x = 3*sin(pi/8*n) x = -cos(5*pi/7*n-pi/3) 2.求数字卷积y = conv(x,h)%x为数字信号，h为滤波器的脉冲响应 或 y = conv(h,x)3.滤波：filter()函数y=filter(b,a,x) 由向量b和a组成的系统对输入x进行滤波4. 求解单位脉冲响应：impz()函数（1）h=impz(b,a) 计算单位脉冲响应的序列值，取样点个数由matlab自动选取（2）h=impz(b,a,n)计算指定范围内（0: n-1）的单位脉冲响应的序列值（3）impz(b,a)绘制单位脉冲响应的时域波形5. 求解频率响应：freqz()函数（1）h,w=freqz(b,a,n)可得到n点频率响应，这n个点均匀地分布在上半单位圆（即 ），并将这n点频率记录在w中，相应的频率响应记录在h中。n最好能取2的幂次方，如果缺省，则n=512。（2）h,w=freqz(b,a,n,whole)均匀选取n个点计算频率响应。（3）h,w=freqz(b,a,n,fs)fs为采样频率（以hz为单位），在0fs/2频率范围内选取n个频率点，计算相应的频率响应。（4）h,w=freqz(b,a,n,whole,fs) 在0fs之间均匀选取n个点计算频率响应。 （5）freqz(b,a)可以直接得到系统的幅频和相频特性曲线。其中幅频特性以分贝的形式给出，频率特性曲线的横轴采用的是归一化频率，即fs/2=1。6. 系统函数的零极点图的绘制：zplane(b,a)7. (1)求取非周期数字信号的频谱 dtft函数 mag,phase,w = dtft(x);%x为数字时域信号x，mag为x的幅度频谱，phase为x的相位频谱，w为数字频率 函数dtft4和函数dtft的功能和格式相同,只是其计算幅度频谱和相位频谱时是以/4为间隔. 数字信号频谱图绘制 plotdtft(mag,phase,w,option)option1、2 plotdtft(mag,phase,w,option,fs)option3、4 (2)求取周期数值信号频谱dfs函数mag,phase,k=dfs(x,period)/x为数字周期信号，period为x的周期绘制频谱图plotdfs(mag,phase,k,option)-option1、2或plotdfs(mag,phase,k,option,fs)-option3、48.求取数字信号频谱 dft函数mag,phase,k = dft(x)x-数字信号mag,phase,k = dft(h)计算数字滤波器频率响应h-数字滤波器脉冲响应9.快速傅里叶变换 fft函数xx = fft(x);mag = abs(xx); hase = angle(xx); k = 0:(length(xx)-1); plotfft(mag,phase,k,1); 四、实验步骤1、熟悉matlab软件基本操作指令。读懂下列matlab程序指令，键入程序并运行，观察运行结果。（1）conv.m% 计算两个序列的线性卷积；clear;n=5;m=6;l=n+m-1;x=1,2,3,4,5;h=6,2,3,6,4,2;y=conv(x,h);nx=0:n-1;nh=0:m-1;ny=0:l-1;subplot(131);stem(nx,x,.k);xlabel(n);ylabel(x(n);grid on;subplot(132);stem(nh,h,.k);xlabel(n);ylabel(h(n);subplot(133); stem(ny,y,.k);xlabel(n);ylabel(y(n);grid on;（2）filter.m;%求利用数字滤波器对信号进行滤波处理后的输出；clear;x=ones(100);t=1:100;b=.001836,.007344,.011016,.007374,.001836; a=1,-3.0544,3.8291,-2.2925,.55075;y=filter(b,a,x);stem(t,y,.);ylabel(y);xlabel(n);（3）impz .m% 计算滤波器的脉冲响应clear;b=.001836,.007344,.011016,.007374,.001836;a=1,-3.0544,3.8291,-2.2925,.55075;h,t=impz(b,a,40); stem(t,h,.);grid on;ylabel(h(n)xlabel(n)（4）filter.m% 计算滤波器的阶跃响应x=ones(100);t=1:100;b=.001836,.007344,.011016,.007374,.001836;a=1,-3.0544,3.8291,-2.2925,.55075;y=filter(b,a,x);plot(t,x,g.,t,y,k-);grid on;ylabel(x(n) and y(n); xlabel(n);（5）dtft.m%求取非周期数字信号频谱并绘图；n = 0:50;x = sin(n*pi/7).*(stepfun(n,0)-stepfun(n,40);mag,phase,w = dtft(x);plotdtft(mag,phase,w,1);（6）dfs.m%求取周期数字信号频谱并绘图n = 0:30;x = sin(n*pi/8);mag,phase,k = dfs(x,16);plotdfs(mag,phase,k,1); figureplotdfs(mag,phase,k,2); figureplotdfs(mag,phase,k,3,8000); figure;plotdfs(mag,phase,k,4,8000); （7）zplane(b,a); % 画出所给系统的极零图b=1 -1.7 1.53 -0.68;a=1;zplane(b,a); （8）dft.m%求取数字信号频谱x = 1 -2 3 0 4 -1 5 0;mag,phase,k = dft(x);plotdft(mag,phase,k,1); （9）fft.m% 应用fft 求频谱；x = rand(1,32) - 0.5; xx = fft(x); mag = abs(xx); phase = angle(xx); k = 0:(length(xx)-1); plotfft(mag,phase,k,1); 2、编程求出下列问题的解 1）、滤波器的差分方程为：yn=xn-0.8xn-1-0.5yn-1求出此滤波器脉冲响应和阶跃响应的前十个采样值。程序：脉冲响应clear;b=1,-0.8;a=1,0.5;h,t=impz(b,a,10);stem(t,h,.);grid on;ylabel(x(n) and y(n);xlabel(n); 阶跃响应x=ones(10);t=1:10;b=1,-0.8;a=1,0.5;y=filter(b,a,x);plot(t,x,g-,t,y,k-);grid on; 2）、系统的脉冲响应为hn=e-n(un-un-3)，用卷积求系统的阶跃响应。 程序：clear;n=10;m=3;l=n+m-1;n=0:2;x=ones(1,10);h=exp(-n);y=conv(x,h);nx=0:n-1;nh=0:m-1;ny=0:l-1;subplot(131);stem(nx,x,.k);xlabel(n);ylabel(x(n);grid on;subplot(132);stem(nh,h,.k);xlabel(n);ylabel(h(n);grid on;subplot(133);stem(ny,y,.k);xlabel(n);ylabel(y(n);grid on; 五、实验讨论和分析程序：1、 差分方程、卷积、z变换和傅里叶变换之间如何进行转换？答：差分方程;a0yn+a1yn-1+a2yn-2+anyn-n=b0xn+b1xn-1+bmxn-m卷积是由输入xn所引起的全部输出yn是所有这些加权脉冲相应之和。即yn=xn*hn只要知道脉冲响应和输入就可以得到输出z变换是把时域信号向频域进行转换x（z）=xnz-n y（z）=ynz-n 傅里叶变换x（）=xne-jn(把z变换中的z-n换成e-jn)2、 边界效应是如何产生的？它对信号的滤波效果有何影响？答：多数情况下，采样开始之前的输入情况是未知的，当脉冲响应与未知的的输入采样点重叠时，由于实际的输出值可能受采样开始之前输入信号的影响，所以无法准确的计算输出。计算的开始和末尾都存在这种现象。仅当输入序列与脉冲响应完全重叠时，计算才有意义，这种现象就是边界效应。当一个系统开始运行或条件改变时，输出需要一些时间过渡到新的稳态。边界效应会产生输出的暂态部分和稳态部分，会影响滤波效果，并且会导致失真现象出现。 实验二 iir和fir数字滤波器设计一、 实验目的通过该设计实验掌数字滤波器设计的一般步骤，掌握利用matlab 软件设计数字滤波器的方法，熟悉sptool工具箱的使用方法。二、 实验设备1、 微型计算机1台；2、 matlab软件1套三、 实验原理一）、滤波器的形状及重要参数理想滤波器的形状是矩形，图 1 给出非理想滤波器。 图 1 滤波器形状 通带：增益高的频率范围，信号可以通过，称为滤波器的通带。阻带：增益低的频率范围，滤波器对信号有衰减或阻塞作用，称滤波器的阻带。滤波器截止频率：增益为最大值的0.707倍时所对应的频率为滤波器截止频率增益通常用分贝（db）表示。增益（db）= 20log(增益）增益为 0.707 时对应 -3db，因此截止频率常被称为 -3db。滤波器的带宽：对于低通滤波器宽带是从0 - 3db 对于高通滤波器宽带是从 - 3db采样频率的一半对于带通滤波器带宽是截止频率之间的频率距离二）加窗低通 fir 滤波器的设计1）在过渡带宽度的中间，选择通带边缘频率(hz): f1=所要求的通带边缘频率+(过渡带宽度)/22）计算 1=2f1/fs，并将此值代入理想低通滤波器的脉冲响应 h1n 中： h1n = sin(n1)/n3）从表中选择满足阻带衰减及其他滤波器要求的窗函数，用表中 n 的公式计算所需要的非零项数目。选择奇数项，这样脉冲响应可以完全对称，避免了滤波器产生相位失真，对于|n|(n-1)/2，计算窗函数wn。4）对于|n|(n-1)/2，从式 hn=h1nwn计算（有限）脉冲响应，对于其他 n 值hn=0，此脉冲响应是非因果的。5）将脉冲响应右移 (n-1)/2，确保第一个非零值在n=0处，使此低通滤波器为因果的。三）、设计低通巴特沃斯滤波器：1) 确定待求通带边缘频率 fp1 hz 、待求阻带边缘频率fs1 hz 和待求阻带衰减 - 20logsdb(或待求阻带增益 20logsdb)。通带边缘频率对应 3db增益。2) 用式 =2f/fs 把由 hz 表示的待求边缘频率转成由弧度表示的数字频率，得到 p1 和s1 。3) 计算预扭曲模拟频率以避免双线性变化带来的失真。由 =2fs tan(/2) 求得 p1和 s1，单位是弧度/秒。4)由已给定的阻带衰减 - 20logs(或增益- 20logs)确定阻带边缘增益 s 。5) 计算所需滤波器的阶数n 取整数。6）把 p1代入 n 阶模拟巴特沃斯滤波器传输函数h(s)中，并对 h(s) 进行双线性变换得到 n 阶数字传输函数 h(z)。滤波器实现所需的差分方程可直接从传输函数 h(s) 求出。四）、低通切比雪夫型滤波器的设计：1）确定待求的通带与阻带边缘频率 fp1 和fs1 、待求的通带边缘增益 20log(1- p) 和待求的阻带衰减-20logs（或待求的阻带增益 20logs ）。2）用公式 =2f/fs 将待求的边缘频率转换为数字频率（用弧度表示），得到 p1 和 s1 。3）对数字频率采用预扭曲以避免双线性变换引起的误差。由 =2fs tan(/2) 得到p1和 s1，单位是弧度/秒。4）由指定的通带边缘增益 20log(1- p) ，确定通带边缘增益 1- p 。计算参数。5）由指定的衰减-20logs（或增益 20logs），确定阻带边缘增益 s 。6）计算所需的阶数n。7）将 p1 和 p 代入 n 阶模拟切比雪夫型滤波器的传输函数 h(s)，并对其进行双线性变换，得到 n 阶数字滤波器传输函数 h(z)。实现滤波器所需的差分方程可由传输函数 h(z) 直接得到。五）主要的matlab函数l fir滤波器设计相关函数：matlab提供了几个函数来实现这些窗函数： w=rect(n); n点的矩形窗 w=triang(n); n点的三角窗 w=hanning(n); n点的汉宁窗 w=hamming(n); n点的哈明窗 w=blackman(n); n点的布莱克曼窗 w=kaiser(n,beta); n点的凯塞窗(beta=0.7865) 画窗函数stem(n,w)窗函数的幅度响应图和相位响应图的绘制:mag,phase,ww = dtft(w);plotdtft(mag,phase,ww,2);理想低通滤波器的脉冲响应：matlab中表示为h1 = omega1*sinc(n*omega1);其中omega1= , w1 是低通滤波器设计实际采用的数字通带边缘频率，单位为弧度，即设计要求的通边缘频率与阻带边缘频率的中点。函数makelp用于得到加窗fir滤波器的脉冲响应:h = makelp(n,omega1,window);h = makelp(n,omega1,kaiser,beta);n为滤波器的项数同理b = makelp(n,omega1,window);b = makelp(n,omega1,kaiser,beta);相应低通脉冲响应的频率响应：mag,phase,w = dtft(h);plotdtft(mag,phase,w,2)选项2用于绘制增益为db ，相位为度。举例：h = makelp(101,0.25,hamming);mag,phase,w = dtft(h);plotdtft(mag,phase,w,2);幅度响应图 相位响应图带通和高通滤波器设计函数makebp. h = makebp(n,omega1,omega0,window);h = makebp(n,omega1,omega0,kaiser,beta);mag,phase,w = dtft(h);plotdtft(mag,phase,w,option)举例：h = makebp(81,0.1,0.6,kaiser,8);mag,phase,w = dtft(h);plotdtft(mag,phase,w,1)举例：h = makebp(51,0.1,1,hanning);mag,phase,w = dtft(h);plotdtft(mag,phase,w,2)带阻滤波器设计函数bandfilt h = bandfilt(n,omegal,omegah,type,window);其中type=0为带阻滤波器，type = 1为带通滤波器。举例：h = bandfilt(101,0.5,0.85,0,hamming);mag,phase,w = dtft(h);plotdtft(mag,phase,w,2)举例：h=bandfilt(101,0.3,0.6,1,hamming);mag,phase,w = dtft(h);plotdtft(mag,phase,w,2)l iir滤波器设计相关函数双线形变换函数b,a = bilinear(num, den, fs)举例：h(s) = 2000/(s2+ 100s + 500) ，采样频率为500 hz.b,a = bilinear(2000,1 100 500,500)结果:b = 0.0018 0.0036 0.0018a = 1.0000 -1.8164 0.8183定义的数字滤波器为:h(z) = (0.0018 +0.0036z-1+ 0.0018z-2)/(1 - 1.8164z-1+ 0.8183z-2)获得巴特沃斯滤波器的系数bka和 ak 的函数:b,a = butter(n,omegap1);其中n为滤波器的阶数，omegap1 = , 为滤波器截止频率。举例b,a = butter(3,0.3);得到b =0.0495 0.1486 0.1486 0.0495a =1.0000 -1.1619 0.6959 -0.1378滤波器的幅度响应图和相位响应图的绘制:mag,phase,w = dtft(b,a);plotdtft(mag,phase,w,1);巴特沃思滤波器阶数的求取函数:n = buttord(omegap1,omegas1,3,stopatt)omegap1 = , 是数字通带边缘频率，omegas1 = , 是数字阻带边缘频率. 参数3指通带边缘频率在 3 db ，stopatt是阻带衰减，单位为 db. 举例n = buttord(0.5,0.6,3,30)得到n = 11带通、高通、带阻巴特沃思滤波器生成函数:b,a = butter(n,omega1 omega2)b,a = butter(n,omegah,high)b,a = butter(n,omega1 omega2,stop)其中omega1,omega2和omegah,w1,w2,和wh是数字通带边缘频率，单位为弧度.切比雪夫1型滤波器的系数 bkand ak获取函数:b,a = cheby1(n,passatt,omegap1);passatt为通带边缘衰减举例 b,a = cheby1(3,0.5,0.3);得到b =0.0444 0.1332 0.1332 0.0444a =1.0000 -1.3867 1.0542 -0.3122滤波器的幅度响应图和相位响应图的绘制:mag,phase,w = dtft(b,a);plotdtft(mag,phase,w,1);切比雪夫1型滤波器的阶数获取函数:n = cheb1ord(omegap1,omegas1,passatt,stopatt)举例n = cheb1ord(0.5,0.6,3,30)得到n = 5切比雪夫1型带通、高通、带阻滤波器生成函数:b,a = cheby1(n,passatt,omega1 omega2)b,a = cheby1(n,passatt,omegah,high)b,a = cheby1(n,passatt,omega1 omega2,stop)数字滤波器对信号进行滤波的matlab函数调用如下.滤波器滤波的调用方式 1）：计算输入信号的滤波输出，向量b和a分别是所采用的滤波器的分子系数向量和分母系数向量。 2）：参数zi指定滤波器的初始条件值，3）或：参数dim指定滤波的维数。四、实验步骤1、第九章和第十章课分别任选一题进行滤波器设计，并画出所设计滤波器的幅度响应图、相位响应图和脉冲响应图。例9.10 程序：h=makebp(151,0.06818,0.363,hamming);stem(h);mag,phase,w=dtft(h);plotdtft(mag,phase,w,2);stem(0:150,h,.k);ylabel(h(n);xlabel(n); 10.10 b.n=6 程序: 巴特沃斯高通 b,a=butter(6,.591,high);mag,phase,w=dtft(b,a);plotdtft(mag,phase,w,2);ylabel(h(n);xlabel(n); 2、仿真一带噪声信号，针对其进行滤波器设计。 带噪声信号：f1=4;f2=25;dt=0.02;t=0:dt:3;x=sin(4*pi*f1*t)+cos(4*pi*f2*t); mag1,phase1,k1=dtft(x);plotdtft(mag1,phase1,k1,3,50);%不带噪声的信号z=awgn(x,10,measured);%加噪声mag1,phase1,k1=dtft(z);figure;plotdtft(mag1,phase1,k1,3,50);%带噪声的信号加噪声之前的信号 加噪声之后的信号 低通滤波器： 程序：h = makelp(25,0.64,blackman);mag,ang,w=dtft(h);plotdtft(mag,ang,w,3,50); 3、用设计得到的数字滤波器对该信号进行滤波处理并画出信号滤波器前后的时域和频域图，进行比较。程序：f1=3;f2=25;dt=0.02;t=0:dt:3;x=sin(2*pi*f1*t)+cos(2*pi*f2*t); %原始信号figure;stem(t,x);%时域分析z=awgn(x,10,measured);%加噪声figure;stem(t,z);%时域分析h = makelp(25,0.64,blackman);%加入低通滤波器y1=filter(h,1,x);figure;subplot(2,1,1);plot(t,x);title(input signal);subplot(2,1,2);plot(t,y1);title(output signal);%加噪声之前的输入输出y2=filter(h,1,z);figure;subplot(2,1,1);plot(t,z);title(input signal);subplot(2,1,2);plot(t,y2);title(output signal);%加噪声之后的输入输出figure;stem(t,y1);figure;stem(t