baizaosheng.m

可以产生白噪声

clear;
randn('state',0)          % set the state of randn
T = 20;                      %计算时间
N =1000;                    %总步数
dt=T/N;                     %计算步长
t=[dt:dt:T];                %计算的时间序列
t2=[dt:dt:T];
n=2;
M=1;
% T is time, N is the number of step,
% M is times with difference Brownian Motion

%%%%%%%%%%%%%%%%%%%% matrix value of system %%%%%%
% A=[0.9 0;0 1.6];
% B=[-0.5 1;-0.2 1];
% W=[0.1 0.5;0.2 0.4];
% C=[0.2 -0.1;0.1 0.2];
% D=C;
% G=[0.1;0.2];
% Ea=[0.2 0.3];
% Eb=[0.3 0.2];
% Ew=[0.2 0.1];
% ================================================
% A=[1.5 0;0 2.1];
% B=[-0.5 1;-0.2 1];
% W=[0.1 0.5;0.2 0.4];
% C=[0.2 -0.1;0.1 0.2];
% D=C;
% G=[0.1;0.2];
% Ea=[0.2 0.3];
% Eb=[0.3 0.2];
% Ew=[0.2 0.1];
%===============================================
A=[1 0;0 1];
B=[2.0 -0.1;-5 2.8];
W=[-1.6 -0.1;-0.3 -2.5];
C=[0 0;0 0];
D=C;
G=[0;0];
Ea=[0 0];
Eb=[0 0];
Ew=[0 0];
%===============================================
r=2;
dB = 0*sqrt(dt)*randn(M,N);   % This is dB,
X1t=zeros(r,M,N);            % format Xt
X1zero(1,:)=0.5*ones(1,M);            %set the inital state
X1zero(2,:)=-0.6*ones(1,M);
x1=zeros(M,N);
x2=zeros(M,N);
%%%%%%%%%%%%%%%%%%%% calculate state of system  %%%%%%%%%%

for i=1:M
    X1temp=X1zero(:,i);
    for j=1:N
        Ft=sin(j);
        Tautemp=1;
        Tau=ceil(Tautemp/dt);
        kk=j-Tau;
        if(kk<=0)
            X1tempDelay=[0.5;-0.6];
        else
            X1TempDelay=X1t(:,i,kk);
        end
        %%%%%%%%%%%%%%%%%%%%%%%%% sys state %%%%%%%%%%% 
        X1temp=X1temp+(-(A+G*Ft*Ea)*X1temp+(B+G*Ft*Eb)*Mytanh(X1temp)+(W+G*Ft*Ew)*Mytanh(X1tempDelay)).*dt+(C*X1temp+D*X1tempDelay).*dB(i,j);  
        X1t(:,i,j)=X1temp;   % get Xt
        x1(i,j)=X1t(1,i,j); % get x1
        x2(i,j)=X1t(2,i,j); % get x2      
    end 
end

figure(3)
lineX1=plot([0,t],[X1zero(1,1),x1(1,:)],'r--'); hold on; % plot x1 against t
lineX2=plot([0,t],[X1zero(2,1),x2(1,:)]); hold on; % plot x2 against t
legend('x_1(t)', 'x_2(t)');
xlabel('Time (Sec)');
ylabel('x(t)');
% figure(2);
% plot(x1,x2);