qr_tls.m

基于QR分解的总体最小二乘方法 基于Matlab的程序

function [ output_args ] = QR_TLS( input_args )
%QR_TLS Summary of this function goes here
%   Detailed explanation goes here
%shuru(n)为实际输入信号，biaozhun(n)为参照信号 
M=1000;%时长 
N=8;%阵元数 
lamda=0.99;%遗忘因子 
w=zeros(N+1,M+1); 
e=zeros(1,M+1); 
w(1,1)=biaozhun(0)/shuru(0); 
 
  for k=1:N 
     for i=1:k 
         sum1=0; 
        for j=1:i 
             sum1=sum1-shuru(j)*w(i-j+1,k+1); 
        end 
         w(i+1,k+1)=(sum1+biaozhun(i))/shuru(0); %式(8.11) 
     end 
  end 
  
 Y=zeros(N+1); 
 for i=1:N+1 
     for j=1:N+2-i 
         Y(i,j)=shuru(N+3-i-j); %式(8.2) 
     end 
 end  
  
  Ul=lamda^(1/2)*Y; %式(8.12) 
  dlq=zeros(N+1,1); 
  for i=1:N+1 
      dlq(i,1)=(lamda^(i/2)*biaozhun(N-i+1)); 
  end 
   
   
  for k=N+1:M 
      rl=1; 
      dl=biaozhun(k); 
      Xl=zeros(1,N+1); 
      for j=1:N+1 
          Xl(1,j)=shuru(k+1-j); 
      end 
       
      for i=0:N 
          c=sqrt(Ul(i+1,N+1-i)^2+Xl(1,N+1-i)^2); 
          cos=Ul(i+1,N+1-i)/c; %式(8.41) 
          sin=Xl(1,N+1-i)/c; %式(8.42)   
           
          s1=sparse([1,i+2],[1,1],[cos,sin],N+2,1); 
          s2=zeros(1,i); 
          s3=eye(i); 
          s4=zeros(N-i+1,i); 
          s5=[s2;s3;s4]; 
          s6=sparse([1,i+2],[1,1],[-sin,cos],N+2,1); 
          s7=zeros(i+2,N-i); 
          s8=eye(N-i); 
          s9=[s7;s8]; 
          Ql=[s1,s5,s6,s9]; %式(8.40) 
      
          A=Ql*[Xl;Ul]; %式(8.39) 
          Xl=A(1,:); 
          Ul=A(2:N+2,:); 
           
           rl=rl*cos; %式(8.54) 
     
           B=Ql*[dl;dlq];%式(8.51) 
           dl=B(1,1); 
           dlq=B(2:N+2,1);  
      end 
  
      dlq=lamda^(1/2)*dlq; 
      Ul=lamda^(1/2)*Ul; 
      r=rl; %式(8.51) 
      e(1,k+1)=dl*r;%式(8.51)%误差 
          
      D=[dl;dlq];%式(8.51) 
      U=Ul; 
      w(1,1)=D(N+2,1)/U(N+1,1); 
 
       
     for i=1:N 
         sum2=0; 
         for j=1:i 
             sum2=sum2-U(N+1-i,i-j+1)*w(i-j+1,k+1); 
         end 
         w(i+1,k+1)=(sum2+D(N+2-i,1))/U(N+1-i,i+1); %式(8.46) 
     end 
      
  end 
 
  figure(8); 
  k=1:M+1; 
  hold on; 
  plot(k,e(1,k),'g'); 
  title('QR-RLS'); 
  grid on; 
  %axis([0,M,-8,8]);