clc;clear ;
global p1 p2; 
p1=0.3,c=150,e=0.01
% 系統隨機信號產生
u=20*rand(1,100)-10;
t=20*rand(1,100)-10;
for k=1:100
    u(k+100)=t(k);
end
%作輸入信號圖形和未知系統輸出圖形
h1=figure;
for n=1:200
    k=1:200;
    subplot(2,1,1);
    plot(k,u(k),'g');
    ylabel('u(k)');title('input rand data');
    hold on
    subplot(2,1,2);
    for t=2:199
    y(1)=0;y(2)=0;
    y(t+1)=y(t).*y(t-1).*(y(t)+1.5)./(y(t).^2+y(t-1).^2+1)+u(t);
end
    plot(k,y(k),'r');
    ylabel('y(k+1)'),title('output of the noknown sysytem')
end
%SVM的訓練
v=zeros(3,100);v(2,1)=0; v(1,1)=0;%v(3,2)=0;v(1,2)=u(1);v(3,1)=0,v(2,2)=0;
%v(1,99)=u(99),v(1,100)=u(100);
for i=3:100;
    v(1,i-1)=u(i-1);
    v(2,i-1)=y(i-1);
    v(3,i-1)=y(i-2);
     s(i)=y(i)';
     s(1)=y(1);s(2)=y(2);
end
ker='rbf';X=v';Y=s'; %Y=v(3,:)';
clear n;clear t;
% c=150*ones(size(X,1),1);
[nsv,beta,bias]=svr(X,Y,ker,c,'einsensitive',e);
% SVM的測試
h2=figure;
for i=2:101
    r(1,i-1)=u(99+i);
    r(2,i-1)=y(i+99);
    r(3,i-1)=y(98+i);%r(1,100)=u(200),r(2,100)=0,r(3,100)=y(200);
end
tstX=r';trnX=v';
tstY=svroutput(trnX,tstX,ker,beta,bias);
for n=1:100
    t=1:100;
    subplot(2,1,1);
    %subplot(2,1,1);
    plot(t,y(t+100)','b',t,tstY(t),'r',t,tstY(t),'k+')
    ylabel('trnY and tstY')
    legend('the actual line','the prediction line','the prediction values')  
    title ('the actual and prediction values')
end
subplot(2,1,2);
%svrplot(X,Y,ker,beta,bias,0.01);
%輸出誤差函數圖形
for n=1:100
    t=1:100 ;
    error(t)=y(t+100)'-tstY(t);
end
plot(t,error,'m');
ylabel('error');
title ('error=actul-prediction')
hold on
tf=find(abs(error)>0.2);
plot(tf,error(tf),'k*')
legend('abs(error)<0.2 ','the fault data',2);
num=length(tf);
fprintf('the fault data number: %d (%3.1f%%)\n',num,100*num/max(t));
%非線性系統辨識(訓練數據的再輸入）
h3=figure;
trnY=svroutput(trnX,trnX,ker,beta,bias);
for n=1:100
    t=1:100;
    subplot(2,1,1);
    plot(t,y(t)','b',t,trnY(t),'r',t,trnY(t),'k+')
    ylabel('y and trnY')
    legend('the actual line','the prediction line','the prediction values')  
    title ('the actual and prediction values')
end
%輸出誤差函數圖形（真實訓練數據和SVM辨識數據之間的誤差）
subplot(2,1,2)
for n=1:95
    w=1:95 ;
    error2(w)=y(w)'-trnY(w);
end
plot(w,error2,'m');
ylabel('error2');
title ('error2=actul-estemt')