2．对题1中的系统，采用非周期测试信号确定系统的频率响应，并与题1的结果对比。
3．基于题1或题2 产生的频率响应数据，采用课堂讲授的频域特性拟合方法，辨识传递函数的参数。将辨识结果与Matlab工具库中的等价的功能函数invfreqs产生的结果做对比。
2.1.2
根据最大长度现行反馈寄存器M序列生成机制，编写M序列生成的生成程序。
T=2*pi./w;
dt = 0.001*T;
tmax =round(10*T);
fori=1:1:1000
t = (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为在5000~6000个点中峰值点的坐标。
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')
holdon
bode(tf(n2,d2),'g')
可以看出，辨识效果非常好，非常接近于系统自带函数freqs()的频响图。
register=init_register;
temp_register=ones(1,N);
new_sequence=ones(1,length1);
new_sequence(1)=init_register(N);%输出第一个M序列值
fori=2:1:length1
temp_register(1)=mod(sum(register.*type_of_connection),2);
y=y2'*exp(-1j*q'*w*dt);
u=u2*exp(-1j*q'*w*dt);
num=[32];
den=[234];
sys=tf(num,den);
dt=0.001;
t=0:0.001:30;
u2=exp(-t);
y2 = lsim(sys,u2,t);
w=[0:0.01:2]*2*pi;
Key Words:System identification, System identification,MATLAB
一、
1.1
在自然科学和社会科学的许多领域中，人们越来越重视对系统进行定量的系统分析、系统综合、仿真、控制和预测。将被控对象模型化，是开展工作的前提和基础。
基于传递函数的经典控制理论分析前需要确定模型的各个参数，忽略工作中的变化，适用于精度要求不高的情况下。而实际中许多控制系统内部参数是未知的，而且工作中有可能会变化，这就需要系统辨识来确定模型参数。通常为机理分析法和测试法相结合来辨识，对于系统机理已知的部分采用机理分析法，机理未知的部分采用测试法。
根据m序列的特征方程：移位寄存器结构为
自相关特性：
1、定义：: {0,1}{1,-1},即(0)=1,(1)=-1
2、M序列{an}与经τ迟延移位序列{an-τ}分别表示：
{an}={a1,a2,…,ap}, {an-τ}={a1 +τ,a2 +τ,…,ap,a1,…,aτ}
3、序列{(an)}的自相关定义为
1．自己设定移位寄存器的级数和初值，产生响应的M序列。
2．绘制题1产生的M序列的自相关函数和功率谱密度图形。
2.2
2.2.1
周期性信号：
周期性信号最典型的就是正弦周期函数，在不同频率下会产生不同的周期，进而有不同的频率响应。将不同频率的正弦波组合起来作为输入信号。
U=
所以：Am=A’/A
Phi=
设计系统函数：
下面为自相关函数的代码：
functionf=Autocorrelation(M)
Np=length(M);
sum=0;
M1=ones(1,Np);
fork=1:1:150
sum=0;
k1=mod(k,Np);
fori=1:1:Np-k1
M1(i)=M(i+k1);
end
forj=1:1:k1
M1(Np-k1+j)=M(j);
结果如下：
绿色为辨识出来的系统的波特图，红色为用系统函数用matlab自带freqs函数画出来的图像。
可以看出这里二者的图像已经较为接近，如果继续增大w的密集程度，图像将更为接近。但计算结果较慢。
非周期
周期激励信号有很多种选择，如单位阶跃信号和脉冲信号，以及高斯信号和指数函数信号。在选择中我纠结了很长时间，开始选择用阶跃信号来做，结果不是很理想。后来想起了指数函数信号，决定拿来试一试。
Abstract:System identification, State estimation and The Principle of Automatic Control are three different disciplines of the modern control theory, which are interpenetrated with one another. 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 method of system identification.
R(τ)={(a1)(a1+τ)+(a2)(a2+τ)+…+(ap-τ)(ap)+…+(ap)(aτ)}/Np
={(a1a1+τ)+(a2a2+τ)+…+(ap-τap)+…+(apaτ)}/Np
=(0的数目-1的数目)/Np
注意到(anan+τ)=(an)(an+τ)
4、根据移位相加特性，{an}{an+τ}仍是M序列中的元素。M序列中一个周期中0的数目与1的数目之差。又由M序列的均衡性可知，0比1的个数少一个
Phi(i)=(n-T(i)/4/dt(i))*0.001*2*pi;%得到相频特性。
end
P=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');
fort=2:1:N
temp_register(t)=register(t-1);%移位
end
register=temp_register;
new_sequence(i)=register(N);
end
end
如：
输入：M=M_sequence([0 0 0 0 1],[0 1 0 0 1]);
结果如下：
2、由于T=2*pi./w，所以每一次频率改变都黑改变周期，进而有不同的激励和响应图像。
3、在十个周期内，对第五和第六个采样周期进行分析，得到周期内的最大值和最大值的坐标。
4、对比激励和响应的最大值的峰值比和相位差，进而得到不同频率的一系列频谱响应。
Matlab代码如下：
w=logspace(-1,2,1000);%w为0.1到100期间的1000的频率点
end
fort=1:1:Np
sum=sum+juan(M(t),M1(t));
end
R(k)=sum/Np;
end
h=1:1:150;
f=R;figure;
p= fft(R,Np);
Pxx=abs(Hpsd);
plot_Pxx=10*log10(Pxx);
step(0:1:Np-1,Pxx);
title('Gonglv')
end
functionf=juan(a,b)
ifa==b
f=1;
else
f=-1;
end
end
结果如下：
自相关函数：
功率谱：
三、
1、实验中发现频域法的频率取点非常重要，如果取的点不适合结果就差别很大。另外，我编写的辨识程序有一定的局限性，比如非周期的激励信号，在实验报告中的系统的辨识系统辨识效果非常好。如果我把系统函数中的各个参数都放大数倍，或者增加阶数，辨识效果可能就没有这么好了，在这个方面可能还需要改进，要达到系统自带函数的辨识度，还有一点差距。
1
本次实验通过使用MATLAB仿真系统辨识的过程，采用频域特性拟合的Levy方法，按要求完成传递函数的辨识。描述实验验证的数据准备、基本过程和实验结果。
结合课堂学习的M序列的特性和生成方法，使用matlab语言产生一个M序列并分析器自相关特性。
二、
2.1
2.1.1
1．自己设定一个稳定系统，采用周期测试信号，测定系统的频率响应。
num=[4 3];
den=[4 5 6];
sys=tf(num,den);
Transfer function:
4 s + 3
---------------