octave_packages/tsa-4.2.4/mvfilter.m

   1
   2 function [x,z]=mvfilter(B,A,x,z)
   3 % Multi-variate filter function
   4 %
   5 % Y = MVFILTER(B,A,X)
   6 % [Y,Z] = MVFILTER(B,A,X,Z)
   7 %
   8 %  Y = MVFILTER(B,A,X) filters the data in matrix X with the
   9 %    filter described by cell arrays A and B to create the filtered
  10 %    data Y.  The filter is a 'Direct Form II Transposed'
  11 %    implementation of the standard difference equation:
  12 %
  13 %    a0*Y(n) = b0*X(:,n) + b1*X(:,n-1) + ... + bq*X(:,n-q)
  14 %                        - a1*Y(:,n-1) - ... - ap*Y(:,n-p)
  15 %
  16 %  A=[a0,a1,a2,...,ap] and B=[b0,b1,b2,...,bq] must be matrices of
  17 %  size  Mx((p+1)*M) and Mx((q+1)*M), respectively.
  18 %  a0,a1,...,ap, b0,b1,...,bq are matrices of size MxM
  19 %  a0 is usually the identity matrix I or must be invertible
  20 %  X should be of size MxN, if X has size NxM a warning
  21 %  is raised, and the output Y is also transposed.
  22 %
  23 % A simulated MV-AR process can be generiated with
  24 %       Y = mvfilter(eye(M), [eye(M),-AR],randn(M,N));
  25 %
  26 % A multivariate inverse filter can be realized with
  27 %       [AR,RC,PE] = mvar(Y,P);
  28 %       E = mvfilter([eye(M),-AR],eye(M),Y);
  29 %
  30 % see also: MVAR, FILTER
  31
  32 %       $Id: mvfilter.m 6981 2010-03-02 23:38:34Z schloegl $
  33 %       Copyright (C) 1996-2003 by Alois Schloegl <a.schloegl@ieee.org>
  34 %
  35 %    This program is free software: you can redistribute it and/or modify
  36 %    it under the terms of the GNU General Public License as published by
  37 %    the Free Software Foundation, either version 3 of the License, or
  38 %    (at your option) any later version.
  39 %
  40 %    This program is distributed in the hope that it will be useful,
  41 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  42 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  43 %    GNU General Public License for more details.
  44 %
  45 %    You should have received a copy of the GNU General Public License
  46 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  47
  48
  49 [ra, ca] = size(A);
  50 [rb, cb] = size(B);
  51 [M,  N ] = size(x);
  52
  53 if (ra~=rb),
  54         fprintf(2,'ERROR MVFILTER: number of rows of A and B do not fit\n');
  55         return;
  56 end;
  57 if nargin<4,
  58         z = []; %zeros(M,oo);
  59 end;
  60
  61 if (M~=ra),
  62         if (N==ra),
  63                 fprintf(2,'Warning MVFILTER: dimensions fit only to transposed data. X has been transposed.\n');
  64                 x = x.';
  65                 %[x,z] = mvfilter(B,A,x,z); x = x.'; return;
  66         else
  67                 fprintf(2,'ERROR MVFILTER: dimensions do not fit\n');
  68                 return;
  69         end;
  70 end;
  71
  72 p  = ca/M-1;
  73 q  = cb/M-1;
  74 oo = max(p,q);
  75
  76 if isempty(z)
  77         z = zeros(M,oo);
  78 else
  79         if  any(size(z)~=[M,oo])
  80                 fprintf('Error MVFILTER: size of z does not fit\n');
  81                 [size(z),oo,M]
  82                 return;
  83         end;
  84 end;
  85
  86 %%%%% normalization to A{1}=I;
  87 if p<=q,
  88         for k=1:p,
  89                 %A{k}=A{k}/A{1};
  90                 A(:,k*M+(1:M)) = A(:,k*M+(1:M)) / A(:,1:M);
  91         end;
  92         A(:,1:M) = eye(M);
  93 else
  94         for k=0:q,
  95                 %B{k}=B{k}/A{1};
  96                 B(:,k*M+(1:M)) = B(:,k*M+(1:M)) / A(:,1:M);
  97         end;
  98 end;
  99
 100 for k = 1:N,
 101         acc = B(:,1:M) * x(:,k) + z(:,1);  % / A{1};
 102         z   = [z(:,2:oo), zeros(M,1)];
 103         for l = 1:q,
 104                 z(:,l) = z(:,l) + B(:,l*M+(1:M)) * x(:,k);
 105         end;
 106         for l = 1:p,
 107                 z(:,l) = z(:,l) - A(:,l*M+(1:M)) * acc;
 108         end;
 109         x(:,k) = acc;
 110 end;
 111