在線等！！！MATLAB編程利用廣義最小二乘法進行系統辨識

clear all;

close all;

clc;

% M array construction

Np=15;r=4;

X1=1;X2=1;X3=1;X4=1;

m_length = r*Np;

a=1;

for i=1:1:m_length

Y4=X4;Y3=X3;Y2=X2;Y1=X1;

X4=Y3;X3=Y2;X2=Y1;

X1=xor(Y3,Y4);

if Y4==0

M(i)=-a;

else

M(i)=a;

end

end

figure;

i=1:1:m_length;

plot(i,M);

% 白噪聲

noise = zeros(1,m_length);

for i=1:1:m_length

temp = noise + 0.5*rands(1,m_length);

noise = temp;

end

noise = noise/12;

%noise = temp;

% parameter of system

n=2;d=1;a1=-1;a2=0.5;b1=1;b2=0.5;

S_U0=0.2;S_Y0=0.2;

% generate u,y

u0=ones(1,m_length)*S_U0;

U = M + u0 + noise;

figure;

i=1:1:m_length;

plot(i,U);

%formulation

y(1)=0;y(2)=0;y(3)=0;

Y(1)=S_Y0+y(1)+noise(1);Y(2)=S_Y0+y(2)+noise(2);Y(3)=S_Y0+y(3) +noise(3);

for k=4:m_length

y(k) = b1*U(k-1-d)+b2*U(k-2-d)-a1*y(k-1)-a2*y(k-2);

Y(k)=S_Y0+y(k)+noise(k);

end

figure;

i=1:1:m_length;

plot(i,Y);

%辨識

% premitive value

c=2;

P = (c^3)*eye(3*n+d);

sita(:,3) = [0,0,0,0,0,0,0]';

alf = 0.995;

% M

%compute U0,Y0

sum_U = 0;sum_Y = 0;

for k=1:1:m_length

sum_U = sum_U + U(k);

sum_Y = sum_Y + Y(k);

end

t_U0 = sum_U/m_length;t_Y0 = sum_Y/m_length;

for k=1:1:m_length

t_u(k) = U(k) - t_U0;

t_y(k) = Y(k) - t_Y0;

end

figure;

i=1:1:m_length;

plot(i,t_u);

figure;

i=1:1:m_length;

plot(i,t_y);

v(1)=0;v(2)=0;v(3)=0;

for k=4:1:m_length

fai = [-t_y(k-1),-t_y(k-2),t_u(k-1-d),t_u(k-2-d),v(k-1),v(k-2),v(k-3)]';

v(k) = t_y(k) - fai'*sita(:,k-1);

W = P*fai;

dan = 1/(alf + fai'*W);

sita(:,k) = sita(:,k-1) + dan*W*v(k);

P = ( P - dan * W*W')/alf;

end

figure;

i = 3:1:m_length;

plot(i,sita(1,3:m_length),'r',i,sita(2,3:m_length),'g',i,sita(3,3:m_length),'b',i,sita(4,3:m_length),'y',i,sita(5,3:m_length),'m',i,sita(6,3:m_length),'c',i,sita(7,3:m_length),'k');