《传递函数频域辨识及m序列生成指南》由会员分享，可在线阅读，更多相关《传递函数频域辨识及m序列生成指南（12页珍藏版）》请在金锄头文库上搜索。
1、传递函数频域辨识及M序列生成指南12自动化 摘要：系统辨识、状态估计和控制理论是现代控制论中相互渗透的三个领域。控制理论的应用离不开系统辨识技术，实际中，许多控制系统的模型在工作中是变化的，为了实现自适应控制，需要系统辨识技术不断更新模型参数。通过学习使用MATLAB软件，初步体验系统辨识方法。关键字：系统辨识，控制理论，MATLAB。Abstract: System identification, State estimation and The Principle of Automatic Control are three different disciplines of the mod
2、ern control theory, which are interpenetratedwithoneanother. In practice, the model of system is changing all the time. To control adaptively, the system model should be update its parameters, by the method of System identification. By learning the using of MATLAB, we are supposed to practice the me
3、thod of system identification.Key Words: System identification, System identification, MATLAB目录一、引言21.1介绍21.2实验目的3二、实验内容和方法32.1实验内容32.1.1实验一:32.1.2实验二:32.2实验步骤32.2.1实验一：3非周期激励信号52.2.2实验二：M序列611三、实验分析分析11一、 引言1.1介绍在自然科学和社会科学的许多领域中，人们越来越重视对系统进行定量的系统分析、系统综合、仿真、控制和预测。将被控对象模型化，是开展工作的前提和基础。基于传递函数的经典控制理论分析
4、前需要确定模型的各个参数，忽略工作中的变化，适用于精度要求不高的情况下。而实际中许多控制系统内部参数是未知的，而且工作中有可能会变化，这就需要系统辨识来确定模型参数。通常为机理分析法和测试法相结合来辨识，对于系统机理已知的部分采用机理分析法，机理未知的部分采用测试法。1.2实验目的本次实验通过使用 MATLAB 仿真系统辨识的过程，采用频域特性拟合的Levy方法，按要求完成传递函数的辨识。描述实验验证的数据准备、基本过程和实验结果。结合课堂学习的M序列的特性和生成方法，使用matlab语言产生一个M序列并分析器自相关特性。二、 实验内容和方法2.1实验内容2.1.1实验一:1 自己设定一个稳定
5、系统，采用周期测试信号，测定系统的频率响应。2 对题1中的系统，采用非周期测试信号确定系统的频率响应，并与题1的结果对比。3 基于题1或题2 产生的频率响应数据，采用课堂讲授的频域特性拟合方法，辨识传递函数的参数。将辨识结果与Matlab工具库中的等价的功能函数invfreqs产生的结果做对比。2.1.2实验二: 根据最大长度现行反馈寄存器M序列生成机制，编写M序列生成的生成程序。1 自己设定移位寄存器的级数和初值，产生响应的M序列。2 绘制题1产生的M序列的自相关函数和功率谱密度图形。2.2实验步骤2.2.1实验一：周期性信号：周期性信号最典型的就是正弦周期函数，在不同频率下会产生不同的周期
6、，进而有不同的频率响应。将不同频率的正弦波组合起来作为输入信号。U=Asin(w*t+)Y=Asin(w*t+)所以： Am=A/A Phi= -设计系统函数：num=4 3;den=4 5 6;sys=tf(num,den);Transfer function: 4 s + 3-4 s2 + 5 s + 6根据这个原理，设计离散的正弦周期信号。1、一个周期内取1000个采样点，采样是个周期。单位采样间隔二0.001*T2、由于 T=2*pi./w，所以每一次频率改变都黑改变周期，进而有不同的激励和响应图像。3、在十个周期内，对第五和第六个采样周期进行分析，得到周期内的最大值和最大值的坐标。4
7、、对比激励和响应的最大值的峰值比和相位差，进而得到不同频率的一系列频谱响应。Matlab代码如下：w=logspace(-1,2,1000); %w为0.1到100期间的1000的频率点T=2*pi./w;dt = 0.001*T;tmax =round(10*T);for i=1:1:1000t = (0:dt(i):tmax(i); u1 = sin(w(i)*t);u=lsim(sys,u1,t); %产生在激励下的相应Am(i) n=max(u(5*1000:6*1000);%n为在50006000个点中峰值点的坐标。Phi(i)=(n-T(i)/4/dt(i)*0.001*2*pi;
8、 %得到相频特性。endP=Am.*cos(-Phi);Q=Am.*sin(-Phi);H=P+1j*Q;n1 d1=freq2levy(H,w,2,1);n2 d2=freq2levy(h,w,2,1);bode(tf(n2,d2),r);hold on;bode(tf(n1,d1),g);结果如下：绿色为辨识出来的系统的波特图，红色为用系统函数用matlab自带freqs函数画出来的图像。可以看出这里二者的图像已经较为接近，如果继续增大w的密集程度，图像将更为接近。但计算结果较慢。非周期激励信号周期激励信号有很多种选择，如单位阶跃信号和脉冲信号，以及高斯信号和指数函数信号。在选择中我纠结了
9、很长时间，开始选择用阶跃信号来做，结果不是很理想。后来想起了指数函数信号，决定拿来试一试。非周期激励信号顾名思义就是不能用周期信号的周期性方法来解决这个问题。像正弦信号这种周期和频率密切相关的信号可以直接由激励响应得到关于不同频率的输出，而非周期信号如指数函数则不能这么做。要进行傅里叶变换将其变换到频域上去。Matlab在矩阵运算方面有着得天独厚的优势，所以这里的傅里叶变换我采用了矩阵的计算方法。，代码非常简单容易理解。y=y2*exp(-1j*q*w*dt);u=u2*exp(-1j*q*w*dt);num=3 2;den=2 3 4;sys=tf(num,den); dt=0.001; t
10、=0:0.001:30; u2=exp(-t); y2 = lsim(sys,u2,t); w=0:0.01:2*2*pi; n=length(w);q=1:30001;y=y2*exp(-1j*q*w*dt);u=u2*exp(-1j*q*w*dt);l=y./u;n d=freq2levy(l,w,2,1); bode(tf(n,d),r)hold onbode(tf(n2,d2),g)可以看出，辨识效果非常好，非常接近于系统自带函数freqs()的频响图。2.2.2实验二：M序列m序列码发生器是一种反馈移位型结构的电路，它由n位移位寄存器加异或反馈网络组成，其序列长度21nM=-，只有一
11、个多余状态即全0状态，所以称为最大线性序列码发生器。由于其结构已定型，且反馈函数和连接形式都有一定的规律，因此利用查表的方式就设计出m序列码。列出部分m序列码的反馈函数F和移存器位数n的对应关系。如果给定一个序列信号长度M，则根据12-=nM求出n，由n查表2-1便可以得到相应的反馈函数F。根据m序列的特征方程：移位寄存器结构为自相关特性： 1、 定义： h: 0,1 1,-1, 即 h(0)=1, h(1)=-1 2、 M序列an与经 迟延移位序列an- 分别表示： an=a1,a2,ap, an-=a1 + ,a2 +,ap,a1,a 3、 序列h(an)的自相关定义为 R()=h(a1)
12、h(a1+)+h(a2)h(a2+) + h(ap-)h(ap)+ h(ap)h(a)/Np =h(a1a1+)+h(a2a2+)+h(ap-ap)+h(apa)/Np =(0的数目-1的数目)/Np 注意到 h(anan+)=h(an)h(an+)4、根据移位相加特性，anan+仍是M序列中的元素。M序列中一个周期中 0 的数目与 1 的数目之差。又由M序列的均衡性可知，0 比 1 的个数少一个 R(0)= 1； R( )=-1/Np 生成M序列的代码如下：function new_sequence=M_sequence(init_register,type_of_connection)N=
13、length(init_register);%判断初始值的长度length1=2N-1; %可产生最大序列的长度register=init_register;temp_register=ones(1,N);new_sequence=ones(1,length1);new_sequence(1)=init_register(N); %输出第一个M序列值for i=2:1:length1 temp_register(1)=mod(sum(register.*type_of_connection),2); for t=2:1:N temp_register(t)=register(t-1); %移位
14、 end register=temp_register; new_sequence(i)=register(N);endend如：输入：M=M_sequence(0 0 0 0 1,0 1 0 0 1);结果如下：下面为自相关函数的代码：function f=Autocorrelation(M)Np=length(M);sum=0;M1=ones(1,Np);for k=1:1:150 sum=0; k1=mod(k,Np); for i=1:1:Np-k1 M1(i)=M(i+k1); end for j=1:1:k1 M1(Np-k1+j)=M(j); end for t=1:1:Np sum=sum+juan(M(t),M1(t); end R(k)=sum/
