Scribd Logo
Upload

Welcome to Scribd!

  • 9 views

    Code Mat

    Uploaded by

    nsh
    sdf

    Copyright:

    © All Rights Reserved
    Download as TXT, PDF, TXT or read online from Scribd

    Code Mat

    Uploaded by

    nsh
    9 views2 pages
    sdf

    Original Title

    code mat

    Copyright

    © © All Rights Reserved

    Available Formats

    TXT, PDF, TXT or read online from Scribd

    Did you find this document useful?

    Is this content inappropriate?

    sdf

    Copyright:

    © All Rights Reserved
    Download as TXT, PDF, TXT or read online from Scribd
    Download as txt, pdf, or txt
    9 views2 pages

    Code Mat

    Uploaded by

    nsh
    sdf

    Copyright:

    © All Rights Reserved
    Download as TXT, PDF, TXT or read online from Scribd
    Download as txt, pdf, or txt
    You are on page 1of 2

    clear all

    close all
    clc
    V=[1 1 1 1 1];
    EP=[];

    for i=1:31
    PRBS(i,1)=V(5);
    EP=[EP;V(5);V(5);V(5);V(5)];
    f=V(3)V(5)(-1);
    V(5)=V(4);
    V(4)=V(3);
    V(3)=V(2);
    V(2)=V(1);
    V(1)=f;
    end
    figure(1)
    plot(PRBS,'o')
    figure(2)
    plot(EP,'*')

    %y(i)=1.5*y(i-1)-0.75*y(i-2)+1*EP(i-1)+0.25*EP(i-2)
    y(1,1)=0;
    y(2,1)=1*EP(1);
    for i=3:124
    y(i,1)=1.5*y(i-1)-0.75*y(i-2)+1*EP(i-1)+0.25*EP(i-2);
    end;
    figure(3)
    plot(y,'r')
    title('ideal output')

    %calculte a1,a2,b1 and b2 using LSE method, define Y and X==>


    %theta=inv(X'*X)*X'*Y
    Y=y(3:124);
    X=[y(2:123) y(1:122) EP(2:123) EP(1:122)];
    theta=inv(X'*X)*X'*Y;

    %output with noise


    yn=y+randn(124,1);
    figure(4)
    plot(1:124,y,'r',1:124,yn,'k')
    legend('y','yn')

    %LSE method with noise


    Yn=yn(3:124);
    Xn=[yn(2:123) yn(1:122) EP(2:123) EP(1:122)];
    thetan =inv(Xn'*Xn)*Xn'*Yn;

    %test of the identified model

    yst(1,1)=0;
    yst(2,1)=1*1;
    for i=3:100
    yst(i,1)=1.5*yst(i-1)-0.75*yst(i-2)+1*1+0.25*1;
    end;
    ystLSE(1,1)=0;
    ystLSE(2,1)=thetan(3)*1;
    for i=3:100
    ystLSE(i,1)=thetan(1)*yst(i-1)+thetan(2)*yst(i-2)+thetan(3)*1+thetan(4)*1
    end;
    legend('yst','ystLSE');
    uGLS=EP; yGLS=yn; precision=1e-4; e11=1000; e22=1;
    while abs(e11-e22)>precision
    %calculate thetaGLS
    YGLS=yGLS(3:124);
    XGLS=[yGLS(2:123) yGLS(1:122) uGLS(2:123) uGLS(1:122)];
    thetaGLS=inv(XGLS'*XGLS)*XGLS'*YGLS;
    %estimate the residu
    e=YGLS-XGLS*thetaGLS; %e:122x1
    %e(i)+f1*e(i-1)=xi(i)==>e(i)=-f1e(i-1)+xi(i)
    E1=e(2:122);
    Xe=-e(1:121);
    f=inv(Xe'*Xe)*Xe'*E1;
    %update U and Y
    ys(i,1)=yGLS(i);
    us(i,1)=uGLS(i);
    for i=2:124
    ys(i,1)=yGLS(i)+f*yGLS(i-1);
    us(i,1)=uGLS(i)+f*uGLS(i-1);
    end;
    yGLS=ys; uGLS=us; e11=e22; e22=e'*e;
    end;

    ystGLS(1,1)=0;
    ystGLS(2,1)=thetaGLS(3)*1;
    for i=3:100
    ystGLS(i,1)=thetaGLS(1)*yst(i-1)+thetaGLS(2)*yst(i-2)+thetaGLS(3)*1+
    (thetaGLS(4))*1;
    end;

    %gradient methode
    yGRAD=yn; uGRAD=EP; parGRAD=thetan; epsGRAD=1e-5; alpha=1e-5; eGRAD=10;
    while eGRAD>epsGRAD
    grad=[0;0;0;0];
    ymGRAD(1,1)=0;
    ymGRAD(2,1)=parGRAD(3)*uGRAD(1);
    for i=3:124
    ymGRAD(i,1)=parGRAD(1)*ymGRAD(i-1,1)+parGRAD(2)*ymGRAD(i-
    2,1)+parGRAD(3)*uGRAD(i-1,1)+parGRAD(4)*uGRAD(i-2,1);
    sensib=2*[-ymGRAD(i-1,1); -ymGRAD(i-2,1); -uGRAD(i-1,1); -uGRAD(i-2,1)];
    end;
    parGRADkp1=parGRAD-alpha*grad;
    eGRAD=abs(max(parGRADkp1-parGRAD));
    parGRAD=parGRADkp1;
    end;
    parGRAD

    y_gr(1,1)=0;
    y_gr(2,1)=1;
    for i=3:100
    y_gr(i,1)=parGRAD(1)*yst(i-1)+parGRAD(2)*yst(i-2)+parGRAD(3)*1+parGRAD(4)*1;

    end;
    figure(5)
    plot(1:100,yst,'r',1:100,ystLSE,'k',1:100,ystGLS,'b',1:100,y_gr,'g')
    legend('process','LSE','GLS','gradient')

    You might also like