% RBF01.m
clear all;
close all;
ts=0.001;
yout_1=0;
b=1.5*ones(1,4);
c=0.5*ones(2,4);
w=rands(1,4);
b_1=b;b_2=b_1;
c_1=c;c_2=c_1;
w_1=w;w_2=w_1;
eta=0.5;
alfa=0.05;
x=[0,0]';

for k=1:1:2000
    time(k)=k*ts;
    % input signal
    u(k)=0.5*sin(2*pi*k*ts);
    % nonlinear plant output
    yout(k)=u(k)^3+yout_1/(1+yout_1^2);
    yout_1=yout(k);      
    % RBF h(j)
    for j=1:1:4
        h(j)=exp(-norm(x-c(:,j))^2/(2*b(j)*b(j)));
    end
    % RBF output
    ynout(k)=w*h';
    % update RBF dw,db,dc
    em(k)=yout(k)-ynout(k);
    for j=1:1:4
        dw(j)=eta*em(k)*h(j);
        db(j)=eta*em(k)*w(j)*h(j)*(b(j)^-3)*norm(x-c(:,j))^2;
        for i=1:1:2
            dc(i,j)=eta*em(k)*w(j)*h(j)*(x(i)-c(i,j))*(b(j)^-2);
        end
    end
    w=w_1+dw+alfa*(w_1-w_2);
    b=b_1+db+alfa*(b_1-b_2);
    c=c_1+dc+alfa*(c_1-c_2);
    w_2=w_1;w_1=w;
    b_2=b_1;b_1=b;
    c_2=c_1;c_1=c;
    % update RBF input vector
    x=[u(k),yout(k)]'; 
end
figure(1);
plot(time,yout,'b',time,ynout,'r');
xlabel('time(s)');ylabel('yout and ynout');
figure(2);
plot(time,em);
xlabel('time(s)');ylabel('error');