系统辨识相关程序
% 最小二乘一次计算方法
%令z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+v(k), 则参数Theta=[-1.5;0.7;1;0.5]
clear
clc
%产生M序列作为输入
x=[0 1 0 1 1 0 1 1 1];
n=403; %n为脉冲数目
M=[]; %M存放M序列
for i=1:n
temp=xor(x(4),x(9));
M(i)=x(9);
for j=9:-1:2
x(j)=x(j-1);
end
x(1)=temp;
end
v=randn(1,400); %生成高斯白噪声
%生成观测序列z
z=zeros(402,1); %初始化
z(1)=-1;
z(2)=0;
for i=3:402
z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i-2);
end
H=zeros(400,4);
for i=3:402
H(i,1)=-z(i-1);
H(i,2)=-z(i-2);
H(i,3)=M(i-1);
H(i,4)=M(i-2);
end
theta=inv(H'*H)*H'*z
%限定记忆最小二乘
%z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+v(k)
clc
%产生M序列作为输入
x=[0 1 0 1 1 0 1 1 1];
n=403; %n为脉冲数
M=[]; %存放M序列
for i=1:n
temp=xor(x(4),x(9));
M(i)=x(9);
for j=9:-1:2
x(j)=x(j-1);
end
x(1)=temp;
end
%产生均值为0，方差为1的高斯白噪声
v=randn(1,402);
%产生观测序列
z=zeros(402,1);
z(1)=-1;
z(2)=0;
for i=3:402
z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i);
end
%递推求解
P_a=10^6*eye(4);
theta_a=[0;0;0;0];
L=20; %记忆长度
for i=3:L-1
h=[-z(i-1);-z(i-2);M(i-1);M(i-2)];
K=P_a*h*inv(h'*P_a*h+1);
theta_a=theta_a+K*(z(i)-h'*theta_a);
P_a=(eye(4)-K*h')*P_a;
end
for k=0:380
hl=[-z(k+L-1);-z(k+L-2);M(k+L-1);M(k+L-2)];
K_b=P_a*hl*inv(hl'*P_a*hl+1);
theta_b=theta_a+K_b*(z(k+L)-hl'*theta_a);
P_b=(eye(4)-K_b*hl')*P_a;
hk=[-z(k+L);-z(k+L-1);M(k+L);M(k+L-1)];
K_a=P_b*hk*inv(hk'*P_b*hk+1);
theta_a=theta_b-K_a*(z(k+L+1)-hk'*theta_b);
P_a=(eye(4)+K_a*hk')*P_b;
end;
disp('参数 a1,a2,b1,b2的估计值为：')
theta_a
%遗忘因子最小二乘一次计算方法
%z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+v(k)
clc
%产生M序列作为输入
x=[0 1 0 1 1 0 1 1 1];
n=403; %n为脉冲数
M=[]; %存放M序列
for i=1:n
temp=xor(x(4),x(9));
M(i)=x(9);
for j=9:-1:2
x(j)=x(j-1);
end
x(1)=temp;
end
v=randn(1,400);
z=zeros(402,1);
z(1)=-1;
z(2)=0;
u=0.98; %遗忘因子
L=400; %数据长度
for i=3:402
z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i-2);
Z(i)=z(i)*u^(L-i+2);
end
H=zeros(400,4);
for i=1:400
H(i,1)=-z(i+1)*u^(L-i);
H(i,2)=-z(i)*u^(L-i);
H(i,3)=M(i+1)*u^(L-i);
H(i,4)=M(i)*u^(L-i);
end
Estimate=inv(H'*H)*H'*(Z(3:402))'
%增广最小二乘的递推算法
%z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)-v(k-1)+0.2*v(k-2)
clc
%产生M序列作为输入
x=[0 1 0 1 1 0 1 1 1];
n=403; %n为脉冲数
M=[]; %存放M序列
for i=1:n
temp=xor(x(4),x(9));
M(i)=x(9);
for j=9:-1:2
x(j)=x(j-1);
end
x(1)=temp;
end
v=randn(1,402);
z=zeros(402,1);
z(1)=-1;
z(2)=0;
for i=3:402
z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.
5*M(i-2)-v(i-1)+0.2*v(i-2);
end
P=10^6*eye(6);
theta=zeros(6,1);
for i=3:402
h=[-z(i-1);-z(i-2);M(i-1);M(i-2);v(i-1);v(i-2)];
K=P*h*inv(h'*P*h+1);
theta=theta+K*(z(i)-h'*theta);
P=(eye(6)-K*h')*P;
end
disp(theta);
%递推最小二乘算法，并显示参数估值的过渡过程
%令z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+v(k), 则参数Theta=[-1.5;0.7;1;0.5]
clear
clc
%产生M序列作为输入
x=[0 1 0 1 1 0 1 1 1];
n=403; %n为脉冲数目
M=[]; %M存放M序列
for i=1:n
temp=xor(x(4),x(9));
M(i)=x(9);
for j=9:-1:2
x(j)=x(j-1);
end
x(1)=temp;
end
v=randn(1,400); %生成高斯白噪声
%生成观测序列z
z=zeros(402,1); %初始化
z(1)=-1;
z(2)=0;
for i=3:402
z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i-2);
end
%递推求解
P=10^6*eye(4); %初始化P
Theta=zeros(4,401); %存放中间过程估值
Theta(:,1)=[0;0;0;0];
for i=3:402
H=[-z(i-1);-z(i-2);M(i-1);M(i-2)];
K=P*H*inv(H'*P*H+1);
Theta(:,i-1)=Theta(:,i-2)+K*(z(i)-H'*Theta(:,i-2));
P=(eye(4)-K*H')*P;
end
i=1:401;
%plot(i,Theta(1,:),i,Theta(2,:),i,Theta(3,:),i,Theta(4,:))
plot(i,Theta(1,:),'r',i,Theta(2,:),'b',i,Theta(3,:),'g',i,Theta(4,:),'m')
legend('theta_1','theta_2','theta_3','theta_4')
%grid on;
title('待估参数过渡过
程')
% 加权最小二乘
%令z(k)=1.5*z(k-1)-0.7*z(k-2)+u(k-1)+0.5*u(k-2)+v(k), 则参数Theta=[-1.5;0.7;1;0.5]
clear
clc
%产生M序列作为输入
x=[0 1 0 1 1 0 1 1 1];
n=403; %n为脉冲数目
M=[]; %M存放M序列
for i=1:n
temp=xor(x(4),x(9));
M(i)=x(9);
for j=9:-1:2
x(j)=x(j-1);
end
x(1)=temp;
end
v=randn(1,400); %生成高斯白噪声
%生成观测序列z
z=zeros(402,1); %初始化
z(1)=-1;
z(2)=0;
for i=3:402
z(i)=1.5*z(i-1)-0.7*z(i-2)+M(i-1)+0.5*M(i-2)+v(i-2);
end
H=zeros(402,4);
W=randn(1,402);
V=diag(W);
for i=3:402
H(i,1)=-z(i-1);
H(i,2)=-z(i-2);
H(i,3)=M(i-1);
H(i,4)=M(i-2);
end
theta=inv(H'*V*H)*H'*V*z
%以下为四种参数辨识方法,都是基于最小二乘原理
t=[24.0 26.5 31.1 35.0 40.3 45.0 49.7 54.9]; %温度t
R=[2.897 2.919 2.969 3.003 3.059 3.107 3.155 3.207]; %电阻R
figure(1)
plot(t,R,'*');
hold on
%方法一：简单方法实现线性拟合
a=polyfit(t,R,1);
b=polyval(a,t);
plot(t,b,'r');
theta=a' %结果为：R=2.6534+0.0101t
%方法二
figure(2)
t=[24.0 26.5 31.1 35.0 40.3 45.0 49.7 54.9]; %温度t
R=[2.897 2.919 2.969 3.003 3.059 3.107 3.155 3.207]; %电阻R
plot(t,R,'*')
hold on
m=7; %m=采样数据总数-1
n=1; %n为需要拟合的阶数
A=zeros(n+1);
for j=1:n+1
for i=1:n+1
for k=1:m+1
A(j,i)=A(j,i)+t(k)^(j+i-2);
end
end
end
B=[0 0];
for j=1:n+1
for i=1:m+1
B(j)=B(j)+R(i)*t(i)^(j-1);
end
end
B=B';
theta=inv(A)*B
t=[24.0 26.5 31.1 35.0 40.3 45.0 49.7 54.9];
R
2=theta(1)+theta(2)*t;
plot(t,R2,'g') %结果为：R=2.6534+0.0101t
%方法三:一般最小二乘算法
figure(3)
t=[24.0 26.5 31.1 35.0 40.3 45.0 49.7 54.9]; %温度t
R=[2.897 2.919 2.969 3.003 3.059 3.107 3.155 3.207]; %电阻R
plot(t,R,'*')
hold on
Z=R';
H=zeros(8,2);
for i=1:8
H(i,1)=1;
H(i,2)=t(i);
end
theta=inv(H'*H)*H'*Z
R3=theta(1)+theta(2)*t;
plot(t,R3,'b') %结果为：R=2.6534+0.0101t
%方法四:递推最小二乘
figure(4)
t=[24.0 26.5 31.1 35.0 40.3 45.0 49.7 54.9]; %温度t
R=[2.897 2.919 2.969 3.003 3.059 3.107 3.155 3.207]; %电阻R
plot(t,R,'*')
hold on
Z=R';
P=10^6*eye(2); %初始化P
theta=zeros(2,1); %初始化Theta
for i=1:8
H=[1;t(i)];
K=P*H*inv(H'*P*H+1);
theta=theta+K*(Z(i)-H'*theta);
P=(eye(2)-K*H')*P;
end
R4=theta(1)+theta(2)*t;
plot(t,R4,'black')
disp('theta=')
disp(theta) %结果为：R=2.6534+0.0101t
%以下为M序列生成程序、自相关函数计算程序、功率谱绘制程序
x=[0 1 0 1 1 0 1 1 1];
n=403; %n为脉冲数
M=[]; %存放M序列
for i=1:n
temp=xor(x(4),x(9));
M(i)=x(9);
for j=9:-1:2
x(j)=x(j-1);
end
x(1)=temp;
end
plot(M);
%计算M序列的自相关函数
figure(2)
R=zeros(1,n);
for t=1:n
for k=1:n-t+1
R(t)=R(t)+M(k)*M(k+t-1);
end
R(t)=R(t)/(n-t+1);
end
plot(R)
xlabel('delay')
ylabel('Amplitude')
figure(3)
c1=xcorr(M,n-1,'coeff') %计算M序列的自相关函数，并作归一化处理
plot(c1)
figure(4)
periodogram(M,[]) %绘制M序列的功率谱
%以下为曲线拟合程序
x=0:0.1:1;
y=[ -0.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.30 11.2 ];
plot(x,y)
grid on
figure(2)
p=polyfit(x,y,3)
z=polyval(p,x);
plot(x,z)
grid on
figure(3)
x=[1:1:30];
y=[9,10,11,12,13,14,13,12,11,9,10,11,12,13,14,12,11,10,9,8,7,8,9,11,9,7,6,5,3,1];
a1=polyfit(x,y,3) %三次多项式拟合%
a2= polyfit(x,y,9) %九次多项式拟合%
a3= polyfit(x,y,15) %十五次多项式拟合%
b1= polyval(a1,x)
b2= polyval(a2,x)
b3= polyval(a3,x)
r1= sum((y-b1).^2) %三次多项式误差平方和%
r2= sum((y-b2).^2) %九次次多项式误差平方和%
r3= sum((y-b3).^2) %十五次多项式误差平方和%
plot(x,y,'*') %用*画出x,y图像%
hold on
plot(x,b1, 'r') %用红色线画出x,b1图像%
hold on
plot(x,b2, 'g') %用绿色线画出x,b2图像%
hold on
plot(x,b3, 'b:o') %用蓝色o线画出x,b3图像%
figure(4)
x=-1.0:0.5:2.0;y=[-4.447,-0.452,0.551,0.048,-0.447,0.549,4.552];
plot(x,y,'*')
xlabel 'x轴'
ylabel 'y轴'
title '散点图'
hold on
m=6;n=3;
A=zeros(n+1);
for j=1:n+1
for i=1:n+1
for k=1:m+1
A(j,i)=A(j,i)+x(k)^(j+i-2)
end
end
end;
B=[0 0 0 0];
for j=1:n+1
for i=1:m+1
B(j)=B(j)+y(i)*x(i)^(j-1)
end
end
B=B';
a=in
v(A)*B
x=[-1.0:0.0001:2.0];
z=a(1)+a(2)*x+a(3)*x.^2+a(4)*x.^3;
plot(x,z)
legend('离散点','y=a(1)+a(2)*x+a(3)*x.^2+a(4)*x.^3')
title('拟合图')
hold on
x=-1.0:0.5:2.0;
y=[-4.447,-0.452,0.551,0.048,-0.447,0.549,4.552];
p=polyfit(x,y,10)
z=polyval(p,x);
plot(x,z,'r')
%所谓高斯白噪声中的高斯是指概率分布是正态函数，而白噪声是指它的二阶矩不相关，一阶矩为常数，是指先后信号在时间上的相关性。这是考查一个信号的两个不同方面的问题。
%高斯白噪声：如果一个噪声，它的幅度分布服从高斯分布，而它的功率谱密度又是均匀分布的，则称它为高斯白噪声。
%产生均值是1，方差为0.5的高斯白噪声序列
N=300;
mean_v=1;
var_v=0.5;
x=mean_v+sqrt(var_v)*randn(1,N);
disp(x)
xm=mean(x) %计算均值
xv=var(x) %计算方差
Cx=cov(x) %计算自协方差
subplot(221)
plot(x)
xlabel('time')
ylabel('Amplitude')
%计算x的自相关函数
subplot(222)
R=zeros(1,N);
for t=1:N
for k=1:N-t+1
R(t)=R(t)+x(k)*x(k+t-1);
end
R(t)=R(t)/(N-t+1);
end
plot(R)
xlabel('delay')
ylabel('Amplitude')
%产生高斯白噪声之方法2
y=normrnd(1,sqrt(0.5),1,N) %生成均值是1，方差是0.5的高斯随机序列
ym=mean(y) %计算均值
yv=var(y) %计算方差
Cy=cov(y) %计算自协方差
subplot(223)
plot(y)
xlabel('time')
ylabel('Amplitude')
%计算y的自相关函数
subplot(224)
M=zeros(1,N);
for t=1:N
for k=1:N-t+1
M(t)=M(t)+y(k)*y(k+t-1);
end
M(t)=M(t)/(N-t+1);
end
plot(M)
xlabel('delay')
ylabel('Amplitude')
figure(2)
subplot(121)
c1=xcorr(x,N-1,'coeff') %计算x的自相关函数，并作归一化处理
plot(c1)
grid on
subplot(122)
c2=xcorr(y,N-1,'coeff') %计算y的自相关函数，并作归一化处理
plot(c2)
grid on
%绘制x和y的功率谱
figure(3)
subplot(211)
periodogram(x,[]) %绘制x的功率谱
subplot(212)
periodogram(y,[]) %绘制y的功率谱
%绘制x和y的功率谱之方法2
figure(4)
subplot(211)
pwelch(x,[]) %绘制x的功率谱
subplot(212)
pwelch(y,[]) %绘制x的功率谱