octave_packages/tsa-4.2.4/adim.m

   1 function [IR, CC, D] = adim(U, UC, IR, CC, arg5);
   2 % ADIM adaptive information matrix. Estimates the inverse
   3 %   correlation matrix in an adaptive way.
   4 %
   5 % [IR, CC] = adim(U, UC [, IR0 [, CC0]]);
   6 %   U   input data
   7 %   UC  update coefficient 0 < UC << 1
   8 %   IR0 initial information matrix
   9 %   CC0 initial correlation matrix
  10 %   IR  information matrix (inverse correlation matrix)
  11 %   CC  correlation matrix
  12 %
  13 %  The algorithm uses the Matrix Inversion Lemma, also known as
  14 %     "Woodbury's identity", to obtain a recursive algorithm.
  15 %     IR*CC - UC*I should be approx. zero.
  16 %
  17 % Reference(s):
  18 % [1] S. Haykin. Adaptive Filter Theory, Prentice Hall, Upper Saddle River, NJ, USA
  19 %     pp. 565-567, Equ. (13.16), 1996.
  20
  21
  22 %       $Id: adim.m 5090 2008-06-05 08:12:04Z schloegl $
  23 %       Copyright (C) 1998-2003 by Alois Schloegl <a.schloegl@ieee.org>
  24 %
  25 %    This program is free software: you can redistribute it and/or modify
  26 %    it under the terms of the GNU General Public License as published by
  27 %    the Free Software Foundation, either version 3 of the License, or
  28 %    (at your option) any later version.
  29 %
  30 %    This program is distributed in the hope that it will be useful,
  31 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  32 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  33 %    GNU General Public License for more details.
  34 %
  35 %    You should have received a copy of the GNU General Public License
  36 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  37
  38
  39
  40 [ur,p] = size(U);
  41
  42 Mode_E = 1;
  43 %if nargin < 4,
  44 %        m  = zeros(size(U)+[0,Mode_E]);
  45 %end;
  46
  47 if nargin<2,
  48         fprintf(2,'Error ADIM: missing update coefficient\n');
  49         return;
  50 else
  51         if ~((UC > 0) & (UC <1)),
  52                 fprintf(2,'Error ADIM: update coefficient not within range [0,1]\n');
  53                 return;
  54         end;
  55         if UC > 1/p,
  56                 fprintf(2,'Warning ADIM: update coefficient should be smaller than 1/number_of_dimensions\n');
  57         end;
  58 end;
  59
  60 if nargin<3,
  61         IR = [];
  62 end;
  63 if nargin<4,
  64         CC = [];
  65 end;
  66 if nargin<5,
  67         arg5 = 6;
  68 end;
  69 if isempty(IR),
  70         IR = eye(p+Mode_E);
  71 end;
  72 if isempty(CC),
  73         CC = eye(p+Mode_E);
  74 end;
  75
  76 D = zeros(ur,(p+Mode_E)^2);
  77 %D = zeros(ur,1);
  78 W  = eye(p+Mode_E)*UC/(p+Mode_E);
  79 W2 = eye(p+Mode_E)*UC*UC/(p+Mode_E);
  80
  81 for k = 1:ur,
  82         if ~Mode_E,
  83                 % w/o mean term
  84                 d  = U(k,:);
  85         else
  86                 % w/ mean term
  87                 d  = [1,U(k,:)];
  88         end;
  89
  90         if ~any(isnan(d)),
  91                 CC = (1-UC)*CC + UC*(d'*d);
  92
  93                 %K = (1+UC)*IR*d'/(1+(1+UC)*d*IR*d');
  94                 v  = IR*d';
  95                 %K = v/(1-UC+d*v);
  96
  97                 if arg5==0;             % original version according to [1], can become unstable
  98                         IR = (1+UC)*IR - (1+UC)/(1-UC+d*v)*v*v';
  99
 100                 elseif arg5==1;         % this provides IR*CC == I, this seems to be stable
 101                         IR = IR - (1+UC)/(1-UC+d*v)*v*v' + sum(diag(IR))*W;
 102
 103                 elseif arg5==6;         % DEFAULT:
 104                         IR = ((1+UC)/2)*(IR+IR') - ((1+UC)/(1-UC+d*v))*v*v';
 105
 106                         % more variants just for testing, do not use them.
 107                 elseif arg5==2;
 108                         IR = IR - (1+UC)/(1-UC+d*v)*v*v' + sum(diag(IR))*W2;
 109
 110                 elseif arg5==3;
 111                         IR = IR - (1+UC)/(1-UC+d*v)*v*v' + eps*eye(p+Mode_E);
 112
 113                 elseif arg5==4;
 114                         IR = (1+UC)*IR - (1+UC)/(1-UC+d*v)*v*v' + eps*eye(p+Mode_E);
 115
 116                 elseif arg5==5;
 117                         IR = IR - (1+UC)/(1-UC+d*v)*v*v' + eps*eye(p+Mode_E);
 118
 119                 end;
 120         end;
 121
 122         %D(k) = det(IR);
 123         D(k,:) = IR(:)';
 124 end;
 125
 126 IR = IR / UC;
 127 if Mode_E,
 128         IR(1,1) = IR(1,1) - 1;
 129 end;
 130
 131