2 function [x,z]=mvfilter(B,A,x,z)
3 % Multi-variate filter function
6 % [Y,Z] = MVFILTER(B,A,X,Z)
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:
13 % a0*Y(n) = b0*X(:,n) + b1*X(:,n-1) + ... + bq*X(:,n-q)
14 % - a1*Y(:,n-1) - ... - ap*Y(:,n-p)
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.
23 % A simulated MV-AR process can be generiated with
24 % Y = mvfilter(eye(M), [eye(M),-AR],randn(M,N));
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);
30 % see also: MVAR, FILTER
32 % $Id: mvfilter.m 6981 2010-03-02 23:38:34Z schloegl $
33 % Copyright (C) 1996-2003 by Alois Schloegl <a.schloegl@ieee.org>
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.
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.
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/>.
54 fprintf(2,'ERROR MVFILTER: number of rows of A and B do not fit\n');
63 fprintf(2,'Warning MVFILTER: dimensions fit only to transposed data. X has been transposed.\n');
65 %[x,z] = mvfilter(B,A,x,z); x = x.'; return;
67 fprintf(2,'ERROR MVFILTER: dimensions do not fit\n');
79 if any(size(z)~=[M,oo])
80 fprintf('Error MVFILTER: size of z does not fit\n');
86 %%%%% normalization to A{1}=I;
90 A(:,k*M+(1:M)) = A(:,k*M+(1:M)) / A(:,1:M);
96 B(:,k*M+(1:M)) = B(:,k*M+(1:M)) / A(:,1:M);
101 acc = B(:,1:M) * x(:,k) + z(:,1); % / A{1};
102 z = [z(:,2:oo), zeros(M,1)];
104 z(:,l) = z(:,l) + B(:,l*M+(1:M)) * x(:,k);
107 z(:,l) = z(:,l) - A(:,l*M+(1:M)) * acc;