load 'pilotdata.mat'
U=U(300:1190);
Y=Y(300:1190);
z=dtrend([Y,U]);
ze=z(300:650,:);
zv=z;
datose=iddata(ze(:,1),ze(:,2),5);
datosv=iddata(zv(:,1),zv(:,2),5);
%
% Linear model with N4SID
%
modelo1=n4sid(datosv,2,'focus','prediction','n4weight','CVA');
[A,B,C,D]=ssdata(modelo1);
yhat1=dlsim(A,B,C,D,zv(:,2));
mserror1=norm(zv(:,1)-yhat1)^2/length(yhat1);
plot([zv(:,1),yhat1])
%
% Linear model with orthonormal bases
%
psi=abs(eig(A));
figure
[theta,mserror2,yhat2]=orthobases(zv,zv,psi);
figure
[A2,B2,C2,D2] = ob2ss(psi,theta);
yhat2=dlsim(A2,B2,C2,D2,zv(:,2));
mserror2=norm(zv(:,1)-yhat2)^2/length(yhat2);
plot([zv(:,1),yhat2])
%
% Wiener model
% Nonlinear block: 3rd order polynomal  (r=3)
% Linear block: second order system   (p=2)
%
g1=zv(:,1);
g2=zv(:,1).^2;
g3=zv(:,1).^3;
y1=orthofilter(psi,zv(:,2));
phi=[-g2,-g3,y1]';
theta2=pinv(phi')*zv(:,1);
ahat=[1;theta2(1:2)];
bhat=theta2(3:4);
g1v=zv(:,1);
g2v=zv(:,1).^2;
g3v=zv(:,1).^3;
vv=ahat(1)*g1v+ahat(2)*g2v+ahat(3)*g3v;
[A3,B3,C3,D3] = ob2ss(psi,bhat);
vvu=dlsim(A3,B3,C3,D3,zv(:,2));
yy=[-0.4:0.01:0.4]';
vvy=ahat(1)*yy+ahat(2)*yy.^2+ahat(3)*yy.^3;
yhat3=interp1(vvy,yy,vvu,'spline');
mserror3=norm(zv(:,1)-yhat3)^2/length(yhat3);
figure
plot([zv(:,1),yhat3])
