ls_car.m

此程序为基于CAR的递推最小二乘算法

%---------------------------------------------------------------------------------
% LS_CAR.m ; Book, Chapter 3 for  A(z)y(t)=B(z)u(t) +v(t) , for given A(z) and B(z)
%Example 3.2.1, August 15 Sunday, 2004
% sigma=0.1, 1.0, LS_CAR_sigma0110.*   k0=17
% st=40, FF=1, sigma=0.10, 1.00
%fprintf('\n----------------------------------------------------------------------');

%----------------------------------------------------------------------

for st=40 %1, 7, 18, 20,22,40 %Set initial state for random variable
    save st st
clear;clf;format short g
load st
Method = 'CAR-RLS';

FF=1;    %  Forgetting factor
sigma=0.1;  %  sigma=0.1, 1.0  % variance of noise


figure(1);  % figure 创建图形窗口
length1=3100; %  length1什么含义
na = 2; nb = 2;  dd=1; %  dd: delay =1
a=[-1.6, 0.8]; b=[0.412, -0.309]; d=[1];  % d=[1]代表什么
%a=[1, 0.412, 0.309];  b=[0, 0.6804,0.6303];

par0=[a,b]';  roots([1,a]);  %  roots 给定一个分量为多项式 P（s）系数的行向量，返回P（s）=0的解
k0=16; %k0=16 
n=length(par0); PP = eye(n)*1e6;   RR = 1;  h = ones(n,1)*1e-6;  par1 = ones(n,1)*1e-6;  % 1e-6 是10的六次方吗，PP是设定P的初值
% Compute the noise-to-signal ratio--------------------------
a1=[1,a]; b1=[0,b];   % a1，b1 表示什么
sy=f_integral(a1,b1);      sv=f_integral(a1,d);   % f_integral 是什么指令
delta_ns = sqrt(sv/sy)*100*sigma;  % delta_ns 是信噪比吗
    ns = [sv,sy,delta_ns]   % ns 表示什么
    
    %fprintf('$\\sigma_v^2=%6.2f^2$, $\\delta_{\\ns}=%6.2f%s \n',sigma,delta_ns,'\%$');  % 文件形式怎么保存？

% ----Generate the input-output data-----------------------------------------
% rand('state',1);    randn('state',1);
% u=(rand(length1,1) - 0.5)*sqrt(12); v=randn(length1,1);   
rand('state',st);  %  randn('state',1); rand（)产生随机数列的格式
eta=f_rn(length1);  %   f_rn 什么函数
u=eta(:,1); v=eta(:,2);   %  eta(:,1)表示什么

Gz=tf(b1,a1,1); Gn=tf(d,a1,1); 
y=lsim(Gz,u) + lsim(Gn,v * sigma); % lsim(Gz,u)返回时间相应

% you can also use the following

%for k = n:length1   
%    y(k)=-par0(1:na)'*y(k-1:-1:k-na) + par0(na+1:n)'*u(k-1:-1:k-nb)+v(k)*sigma;
%end
%----RLS----------------------------------------------
jj = 0; j1 = 0; j2 = 0;
for k = 20:length1 %　length1=3100  
    jj=jj+1;
    h = [-y(k-1:-1:k-na); u(k-1:-1:k-nb)];   %delay = 1
    yy=y(k);
    [par1,delta,PP,RR] = f_ls(par0,par1,h,yy,FF,'LS',PP,RR);  % f_ls 是什么意思
    % [par1,delta,PP,RR] = f_ls(par0,par1,h,yy,0.7,'SG',PP,RR); 
    ls(jj,:)=[jj, par1',delta];
    if (jj==100)|(jj==200)|(jj==300)|mod(jj,500)==0
        j1 = j1+1;
        ls_100(j1,:)=[jj, par1', delta*100];
    end
    ls_100(j1+1,:)=[ 0, par0', 0];
end
fprintf('\nst = %d,  Method = %s, ($\\sigma_v^2=%5.2f^2$, $\\delta_{\\ns}=%6.2f\\%s$)', st, Method, sigma,delta_ns,'%');  
fprintf('\n %s  \n','$k$  &  $a_1$  &  $a_2$  &  $b_1$  &  $b_2$  &  $\delta\ (\%)$ \\');  
fprintf('%5d &%10.5f &%10.5f &%10.5f &%10.5f &%10.5f  \\\\\n',ls_100');


% st=1:40, fprintf('%d   %9.5f   %9.2f   \\\\\n',st,ls_100(9,n+2),sigma);


%figure(1)
jk=(17:10:2990)';
plot(ls(jk,1), ls(jk,n+2));

end

fprintf('=====================================================================\n')

if sigma==0.1
    z1=[ls(:,1), ls(:,n+2)];
    save z1 z1
else  %sigma==1.0
    load z1
    z0=[z1, ls(:,n+2)];
    figure(2)
    plot(z0(jk,1),z0(jk,2),'k',z0(jk,1),z0(jk,3),'b')
    axis([0, 3000, 0, 0.33])
    text(600,0.058,'{\it\sigma_v^2} = 1.00^2')
    text(600,0.13,'{\it\sigma_v^2} = 0.10^2')
    line([247,620],[0.024,0.119])
%    legend('{\it\sigma_v} = 0.10','{\it\sigma_v} = 1.00')       
    end
xlabel('\it        t'); ylabel('{\it    \delta}');