信息工程学院电子科学与技术系徐利民课件

信息工程学院电子科学与技术系徐利民用MATLAB进行信号与系统分析1234MATLAB的使用信号与系统的数值化处理工程应用基石：函数与工具箱应用课题的研究、分析与设计MATLABMATLAB是由美国MathWorks公司推出的软件产品。MATLAB是“MatrixLaboratory”的缩写，意即“矩阵实验室”。主界面MATLABMATLAB是一个完整的并可扩展的计算机环境，是一种进行科学和工程计算的交互式程序语言。MathWorks公司为MATLAB配备了涉及到通信系统、自动控制、信息处理、计算机仿真等种类繁多的工具箱(Toolbox)。Matlab集成度高、功能强、使用方便和适用的计算机平台宽，在教学、科研和工程实际中得到广泛应用。Matlab的特点变量以矩阵形式出现，每个元素都可以是复数，适合当前科学计算的需要。如：x=[2pi/2sqrt(3)3+5i]a=[123]+[321]*iMatlab的特点笔算式语言规则和人类书写习惯近似，易写易记，键入命令立即执行，无须编译。程序运行中可设置断点，单步调试，遇有错误立即报告，方便编程。编辑窗口Matlab的特点有丰富的作图指令，能绘制各种坐标系，自动确定坐标，曲线曲面调整方便，可设置各种颜色、线形和视角。t=0:pi/50:4*pi;y0=exp(-t/3);y=exp(-t/3).*sin(3*t);plot(t,y,'-r',t,y0,':b',t,-y0,':b');Matlab的特点Theta=linspace(0,2*pi,361);Beita=pi*sin(Theta);A=abs(sin(3*Beita)./(6*sin(0.5*Beita)));polar(Theta,A);

Matlab的特点功能扩展能力强，有大量工具箱可予以选用，并且还在不断发展中。参考文献张志涌.精通MATLAB6.5版.北京:北京航天航空大学出版社,2003谷源涛.信号与系统-MATLAB综合实验.北京:高等教育出版社,2008梁洪.信号与系统分析与MATLAB实现.北京:电子工业出版社,2001孟桥董志芳王琼信号与系统MATLAB实践北京高等教育出版社2008参考文献(美)James.H.McClellan等.栾晓明译.基于计算机的信号处理实践.北京:电子工业出版社,2006徐利民,舒君等.基于MATLAB的信号与系统实验教程.北京:清华大学出版社,2010用MATLAB进行信号与系统实践1234MATLAB的使用信号与系统的数值化处理工程应用基石：函数与工具箱应用课题的研究、分析与设计二、信号与系统的数值化处理基础1：积分的数值计算例1：计算f=sin(t)在[0，π]的积分。…………t0f(t)矩形法计算积分的思想二、信号与系统的数值化处理矩形法：把积分看作1、定积分从左和从右开始计算误差大小？2、梯形法：把每个矩形变为梯形计算，误差多少？3、抛物线法（Simpson法）的计算原理？误差改善情况？4、先进行插值如何？

二、信号与系统的数值化处理例1：计算f=sin(t)在[0，π]的积分。t=linspace(0,pi,1001);x=sin(t);y1=sum(x)*t(2),%矩形法y2=trapz(x)*t(2),%梯形法f=inline('sin(t)'),y3=quad(f,0,pi),%抛物线(Simpson)法运行结果：y1=2.0000y2=2.0000f=Inlinefunction:f(t)=sin(t)y3=2.0000数值分析指令请参考教材第一部分二、信号与系统的数值化处理可以利用积分处理：卷积→卷积和傅立叶级数傅立叶变换例2：数值计算f1(t)=e-tε(t)，f2(t)=te-tε(t)的卷积。(0<t<6)注：计算离散卷积和用指令conv实现二、信号与系统的数值化处理卷积定义为：第一步：将积分变量数值化：第二步：将参变量t数值化：离散卷积和二、信号与系统的数值化处理例2:数值计算f1(t)=e-tε(t)，f2(t)=te-tε(t)的卷积。(0<t<6)dt=0.01;k1=0:dt:6;k2=k1;f1=exp(-k1);%生成信号f1f2=k2.*exp(-k2);%生成信号f2f=dt*conv(f1,f2);%计算卷积结果fk0=k1(1)+k2(1);%计算序列f非零样值的起点位置k3=length(f1)+length(f2)-2;%计算卷积和f的非零样值的宽度k=k0:dt:k0+k3*dt;%确定卷积和f非零样值的时间向量二、信号与系统的数值化处理subplot(2,2,1);%在子图1绘f1(t)时域波形图plot(k1,f1);title('f1(t)');xlabel('t');subplot(2,2,2);%在子图2绘f2(t)时波形图plot(k2,f2);title('f2(t)');xlabel('t');subplot(2,2,3);%画卷积f(t)的时域波形plot(k,f);h=get(gca,'position');h(3)=2.5*h(3);set(gca,'position',h);%将第三个子图横坐标范围扩至原来的2.5倍title('f(t)=f1(t)*f2(t)');xlabel('t');二、信号与系统的数值化处理运行结果：二、信号与系统的数值化处理例3：计算傅立叶系数。%计算部分程序t=linspace(0,1,101);n=0:9;gt=exp(-i*pi*n'*t);Fn=0.5*trapz(gt,2)*t(2),sig=abs(Fn)>0.00001;Fn=Fn.*sig;1-1T=2f(t)1二、信号与系统的数值化处理%绘图部分程序figure;subplot(2,1,1);stem(n,abs(Fn));xlabel('n'),ylabel('Fn');title('振幅谱');grid;subplot(2,1,2);stem(n,angle(Fn));xlabel('n'),ylabel('Theta');title('相位谱');grid;

二、信号与系统的数值化处理计算Fourier变换有异曲同工之妙：限定有限区间计算：对t采样N个点，采样间隔为：得到：

二、信号与系统的数值化处理设待计算的频率区间为：[ω1,ω2]，计算其间均匀采样的K个值，则有：

注意：1、抽样→频谱混叠2、加窗→频谱泄露二、信号与系统的数值化处理Sa(20t)与频谱二、信号与系统的数值化处理cos(20t)与频谱二、信号与系统的数值化处理实际使用中往往借助函数fft()计算频谱。如果想理解这些工程基本知识，必须切实理解四类对应性：时域连续周期←→频域离散非周期(FS)时域连续非周期←→频域连续非周期(FT)时域离散非周期←→频域连续周期(DTFT)时域离散周期←→频域离散周期(DFS|DFT|FFT)二、信号与系统的数值化处理基础2：微分的数值计算利用以下公式将微分化为差分：后向Euler法二、信号与系统的数值化处理例：数值法求解微分方程，初始条件为：y(0-)=y‘(0-)=1，步长取TS=0.1，求解范围为区间[0,4]。解：把微分方程转化为差分方程为这些系数的转化如何不用人工计算？二、信号与系统的数值化处理根据此方程求出差分方程的初始条件：求解差分方程需要初始条件:y(-1),y(-2),怎么办？上例解得y(-1)=1；y(-2)=0.9。编程如下：

二、信号与系统的数值化处理Ts=0.1;%时间间隔n=1:4/Ts+2;%坐标1,2代表初始条件坐标-2,-1F(n)=1;F(1)=0;F(2)=0;%激励信号初始条件Y(1)=0.9;Y(2)=1;%输入初始条件y(-2)=0.9;y(-1)=1fork=3:length(n)%迭代求解Y(k)=(13*F(k)-10*F(k-1)-100*Y(k-2)+240*Y(k-1))/143;endn1=(n-3)*Ts;subplot(2,1,1);plot(n1,Y);xlabel('t');ylabel('y(t)');title('全响应(数值计算)');gridon;%下图用实际响应进行对比subplot(2,1,2);ezplot('-exp(-3*t)+exp(-t)+1',[0,4]);xlabel('t');ylabel('y(t)');title('全响应');gridon;二、信号与系统的数值化处理运行结果：二、信号与系统的数值化处理思考：1、如果连续系统给出的是初始条件y(0+)；y‘(0+)，那么求离散系统初始条件的方程组应该是什么？2、基于状态方程则只用一阶微分式就足够，请思考原理。3、如果计算微分用起点和终点斜率的平均值计算，则可以是计算精度由二阶提高到三阶(梯形预估-校正法)。二、信号与系统的数值化处理4、如果在每一时间区间内多取几个点的斜率值，将它们进行线性组合，则可能构造出精度更高的算法（龙格-库塔法）这些知识很重要，在MATLAB和SIMULINK中多有用到。请一定理解其原理！！！二、信号与系统的数值化处理例9‑3：求解微分方程初值问题的数值解，求解范围为区间[0,0.5]。，程序如下：fun=inline('-2*y+2*x^2+2*x','x','y');[x,y]=ode23(fun,[0,0.5],1);plot(x,y,'o-')运行结果：二、信号与系统的数值化处理SolverODE类型特点说明ode45非刚性单步算法；4、5阶Runge-Kutta方程；累计截断误差达(Δx)3大部分场合的首选算法ode23非刚性单步算法；2、3阶Runge-Kutta方程；累计截断误差达(Δx)3使用于精度较低的情形ode113非刚性多步法；Adams算法；高低精度均可到10-3~10-6计算时间比

ode45短ode23t适度刚性采用梯形算法适度刚性情形ode15s刚性多步法；Gear's反向数值微分；精度中等若

ode45失效时，可尝试使用ode23s刚性单步法；2阶

Rosebrock算法；低精度当精度较低时，计算时间比

ode15s短ode23tb刚性梯形算法；低精度当精度较低时，计算时间比

ode15s短二、信号与系统的数值化处理SIMULINK中首先要设置的参数：二、信号与系统的数值化处理用MATLAB进行信号与系统实践1234MATLAB的使用信号与系统的数值化处理工程应用基石：函数与工具箱应用课题的研究、分析与设计函数functiony=linspace(d1,d2,n)%LINSPACELinearlyspacedvector.%LINSPACE(X1,X2)generatesarowvectorof100linearly……ifnargin==2n=100;endy=[d1+(0:n-2)*(d2-d1)/(floor(n)-1)d2];函数头输入变量输出变量说明行函数代码，一般要对输出变量赋值H1行函数名，需要与.m文件同名函数1、大量数学函数1）三角函数和双曲函数sin正弦cos余弦tan正切asin反正弦acos反余弦atan反正切sinh双曲正弦cosh双曲余弦tanh双曲正切asinh反双曲正弦acosh反双曲余弦atanh反双曲正切cot余切sec正割csc余割acot反余切asec反正割acsc反余割coth双曲余切sech双曲正割csch双曲余割acoth反双曲余切asech反双曲正割acsch反双曲余割atan2四象限反正切函数2）指数函数exp指数log10常用对数pow22的幂log自然对数log22为底的对数sqrt求平方根3）复数函数abs模；绝对值conj共轭复数real复数的实部imag复数的虚部angle相角（弧度）4）圆整函数和求余函数ceil向+∞圆整函数floor向-∞圆整函数fix向0圆整函数round向最近整数圆整mod模除求余rem求余数sign符号函数函数2、线性代数与矩阵运算det行列式的值null零空间rref转换为行阶梯形diag创建对角阵poly特征多项式tril抽取下三角阵inv矩阵的逆rank秩triu抽取上三角阵jordan<Jordan>分解trace迹eig特征分解svd奇异值分解…3、统计分析指令4、绘图指令5、音频控制信号与系统相关函数1）信号生成sinc抽样函数sawtooth锯齿波square方波stepfun阶跃函数gensig函数产生器sign符号函数2）连续系统响应3）离散系统响应impulse冲激响应impz单位样值响应(同dimpulse)step阶跃响应dstep阶跃函数lsim系统响应dlsim系统响应信号与系统相关函数4）相关（Correlation）corrcoef相关系数cov协方差矩阵subspace子空间之间的角度

5）滤波和卷积（Filteringandconvolution）conv卷积和多项式相乘convnN维卷积deconv解卷和多项式相除conv2二维卷积detrend去除线性分量filter一维数字滤波器filter一维数字滤波器6）傅立叶变换（Fouriertransforms）fft快速离散傅立叶变换fftshift重排fft和fft2的输出fft2二维离散傅立叶变换ifft离散傅立叶反变换fftnN维离散傅立叶变换ifft2二维离散傅立叶反变换ifftnN维离散傅立叶反变换ifftshift反fftshift信号与系统相关函数7）系统分析bode波特图rlocus根轨迹freqsH(s)频率特性freqzH(z)频率特性pzmap零极图(H(s))zplane零极图(H(z))nyquistNyquist图roots求多项式根8）三大变换fourier傅立叶变换ifourier傅立叶逆变换laplaceLaplace变换ilaplaceLaplace逆变换ztransZ变换Iztrans逆Z变换信号与系统相关函数9）系统表示zpk利用零极点增益组生成LTI对象zpkdata从LTI对象获取零极点增益ss利用状态方程组生成LTI对象ssdata从LTI对象获取状态方程tf利用传递函数组生成LTI对象tfdata从LTI对象获取传递函数minreal消去传递函数分子、分母公因子simple表达式化简10）系统形式转换ss2tf状态方程转换为传输函数tf2ss传输函数转换为状态方程tf2zp传输函数转换为零极点zp2tf零极点转换为传输函数ss2zp状态方程转换为零极点zp2ss零极点转换为状态方程工具箱常用工具箱：通信工具箱(comm)信号处理工具箱(signal)滤波器设计工具箱(filterdesign)图像处理工具箱等(images)三、工程化处理：函数与工具箱用MATLAB进行信号与系统实践1234MATLAB的使用信号与系统的数值化处理工程应用基石：函数与工具箱应用课题的研究、分析与设计四、应用课题的研究、分析与设计课题1：洛伦茨方程——混沌现象美国气象学家洛伦兹（E.N.Lorenz）是混沌理论的奠基者之一。20世纪50年代末到60年代初，他的主要工作目标是从理论上进行长期天气预报。他在使用计算机模拟天气时意外发现，对于天气系统，哪怕初始条件的微小改变也会显著影响运算结果。随后，他在同事工作的基础上化简了自己先前的模型，得到了有3个变量的一阶微分方程组，由它描述的运动中存在一个奇异吸引子，即洛伦兹吸引子。洛伦兹的工作结果最初在1963年发表，论文题目为DeterministicNonperiodicFlow，发表在JournaloftheAtmosphericSciences杂志上。如今，这一方程组已成为混沌理论的经典。课题1：洛伦茨方程——混沌现象洛伦茨方程：参数a=10,b=8/3,r=28是混沌理论中一经典方程课题1：洛伦茨方程——混沌现象用下面一段MATLAB程序求解洛伦茨方程：functionxdot=lorenzeq(t,x)rou=10;r=28;b=8/3;xdot=[rou*(x(2)-x(1));x(1)*(r-x(3))-x(2);-b*x(3)+x(1)*x(2)];%这一部分保存成lorenzeq.mt_final=20;x0=[6.0004;-1;10];[t,x]=ode45('lorenzeq',[0,t_final],x0);figure(1);plot(t,x),legend('x1','x2','x3');figure(2);plot3(x(:,1),x(:,2),x(:,3));gridon;课题1：洛伦茨方程——混沌现象解曲线轨道就集中在形式非常复杂的一个吸引子上----Lorenz吸引子。它是三维空间中的一条曲线，这条曲线相互缠绕而互不相交象一只展翅飞翔的蝴蝶课题1：洛伦茨方程——混沌现象系统运行到10秒左右状态变量的绝对误差按照指数规律上升到最大值状态变量的变化曲线误差曲线四、应用课题的研究、分析与设计课题2：根据300年间太阳黑子变化数据，研究太阳黑子运动规律研究300年间太阳黑子变化数据近50年太阳黑子变化数据课题2：太阳黑子运动规律研究计算数据的功率谱Y=fft(wolfer);Y(1)=[];n=length(Y);power=abs(Y(1:n/2)).^2;nyquist=1/2;freq=(1:n/2)/(n/2)*nyquist;plot(freq,power)xlabel('cycles/year')title('Periodogram')结论：太阳黑子最强爆发周期是11.0385年，次强活动还有2个周期，都聚集在10-12年之间。四、应用课题的研究、分析与设计课题3：股票预测根据一断时间的股票数据，用Fourier变换研究其周期性。如果有较强周期，则用相同正弦曲线拟和。研究发现此方法有一定的预测能力。股票交易软件界面课题3：股票预测读入股票数据并进行Fourier变换[A,B]=xlsread('600808.xls');x=A(:,5);R(1)=0;N=length(x);fork=1:Nc(k)=x(k)-x(1)+(x(1)-x(N))*(k-1)/(N-1);endc=c-mean(c);Xc=(fft(c,2048));subplot(211);plot(abs(Xc));title('频谱');subplot(212);stem(abs(Xc(1:50)));title('局部放大频谱');课题3：股票预测低频部分有两个峰值，局部放大后确定峰值位置为k=9和k=25上图是全部频谱，下图是低频部分放大课题3：股票预测由k=9和k=25可以确定周期信号的频率：用两个正弦信号近似股价曲线：其中R1、I1、R2、I2分别是对应k=9和k=25时的傅立叶系数的实部和虚部。

课题3：股票预测k=1:N;R1=real(Xc(9));I1=imag(Xc(9));R2=real(Xc(25));I2=imag(Xc(25));cc=(2/N)*(R1*cos(2*9*pi*k/2048)-I1*sin(2*9*pi*k/2048)+(R2*cos(2*25*pi*k/2048)-I2*sin(2*25*pi*k/2048)));pl