# HG changeset patch
# User schloegl
# Date 1056482557 0
# Node ID 9c2a30bbaa0bdd25bc2b6a7f76792aeb1c1bc548
# Parent  8fe7ccc464348c225306857a2882cef64287e997
adim and mvfreqz added

diff -r 8fe7ccc46434 -r 9c2a30bbaa0b extra/tsa/adim.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/extra/tsa/adim.m	Tue Jun 24 19:22:37 2003 +0000
@@ -0,0 +1,131 @@
+function [IR, CC, D] = adim(U, UC, IR, CC, arg5); 
+% ADIM adaptive information matrix. Estimates the inverse
+%   correlation matrix in an adaptive way. 
+%
+% [IR, CC] = adim(U, UC [, IR0 [, CC0]]); 
+%   U 	input data  
+%   UC 	update coefficient 0 < UC << 1
+%   IR0	initial information matrix
+%   CC0 initial correlation matrix
+%   IR	information matrix (inverse correlation matrix)
+%   CC  correlation matrix
+% 	
+%  The algorithm uses the Matrix Inversion Lemma, also known as 
+%     "Woodbury's identity", to obtain a recursive algorithm.  
+%     IR*CC - UC*I should be approx. zero. 
+%
+% Reference(s):
+% [1] S. Haykin. Adaptive Filter Theory, Prentice Hall, Upper Saddle River, NJ, USA 
+%     pp. 565-567, Equ. (13.16), 1996.
+
+
+%       $Revision$ 
+%       $Id$
+%       Copyright (c) 1998-2003 by  Alois Schloegl <a.schloegl@ieee.org>
+%
+
+%
+% This library is free software; you can redistribute it and/or
+% modify it under the terms of the GNU Library General Public
+% License as published by the Free Software Foundation; either
+% version 2 of the License, or (at your option) any later version.
+% 
+% This library is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+% Library General Public License for more details.
+%
+% You should have received a copy of the GNU Library General Public
+% License along with this library; if not, write to the
+% Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+% Boston, MA  02111-1307, USA.
+
+
+[ur,uc] = size(U);
+p = uc;
+
+Mode_E = 1;
+if nargin < 4,
+        m  = zeros(size(U)+[0,Mode_E]);
+end;
+
+if nargin<2,
+        fprintf(2,'Error ADIM: missing update coefficient\n');
+        return;
+else	
+        if ~((UC > 0) & (UC <1)),
+                fprintf(2,'Error ADIM: update coefficient not within range [0,1]\n');
+                return;
+        end;
+        if UC < 1/uc,
+                fprintf(2,'Warning ADIM: update coefficient should be smaller than 1/number_of_dimensions\n');
+        end;
+end;
+
+if nargin<3,
+        IR = [];
+end;
+if nargin<4,
+        CC = [];
+end;
+if nargin<5,
+        arg5 = 1;
+end;
+if isempty(IR),
+        IR = eye(p+Mode_E);
+end;
+if isempty(CC),
+        CC = eye(p+Mode_E);
+end;
+
+D = zeros(ur,(p+Mode_E)^2);
+%D = zeros(ur,1);
+W = eye(p+Mode_E)*UC/(p+Mode_E);
+W2= eye(p+Mode_E)*UC*UC/(p+Mode_E);
+
+for k = 1:ur,
+        if ~Mode_E,
+                % w/o mean term 
+                d  = U(k,:);	
+        else
+                % w/ mean term 
+                d  = [1,U(k,:)];
+        end;
+        
+        if ~any(isnan(d)),
+                CC = (1-UC)*CC + UC*(d'*d);
+                
+                %K = (1+UC)*IR*d'/(1+(1+UC)*d*IR*d');
+                v  = IR*d';
+                %K = v/(1-UC+d*v);
+                
+                if arg5==0;  		% original version accoring to [1], can become unstable 
+                        IR = (1+UC)*IR - (1+UC)/(1-UC+d*v)*v*v';
+                        
+                elseif arg5==1;  	% DEFAULT: this provides IR*CC == I, this seems to be stable 
+                        IR = IR - (1+UC)/(1-UC+d*v)*v*v' + sum(diag(IR))*W;
+                        
+                        		% more variants just for testing, do not use them. 
+                elseif arg5==2;
+                        IR = IR - (1+UC)/(1-UC+d*v)*v*v' + sum(diag(IR))*W2;
+                        
+                elseif arg5==3;
+                        IR = IR - (1+UC)/(1-UC+d*v)*v*v' + eps*eye(p+Mode_E);
+                        
+                elseif arg5==4;
+                        IR = (1+UC)*IR - (1+UC)/(1-UC+d*v)*v*v' + eps*eye(p+Mode_E);
+                        
+                elseif arg5==5;
+                        IR = IR - (1+UC)/(1-UC+d*v)*v*v' + eps*eye(p+Mode_E);
+                        
+                end;
+        end;
+        
+        %D(k) = det(IR);
+        D(k,:) = IR(:)';
+end;
+
+% IR = IR / UC;
+% IR * CC == UC * I;
+
+
diff -r 8fe7ccc46434 -r 9c2a30bbaa0b extra/tsa/contents.m
--- a/extra/tsa/contents.m	Tue Jun 24 17:19:23 2003 +0000
+++ b/extra/tsa/contents.m	Tue Jun 24 19:22:37 2003 +0000
@@ -1,8 +1,10 @@
-% Time Series Analysis (Ver 3.10)
-% Schloegl A. (1996-2003) Time Series Analysis - A Toolbox for the use with Matlab. 
+% Time Series Analysis (Ver 3.20)
+% Schloegl A. (1996-2003) Time Series Analysis - A toolbox for the use with Matlab. 
 % WWW: http://www.dpmi.tu-graz.ac.at/~schloegl/matlab/tsa/
 %
-%	Version 3.10    Date: 25. Feb 2003
+%	$Revision$ 
+%	$Id$
+%	Version 3.21    Date: 23 Juni 2003
 %	Copyright (C) 1996-2003 by Alois Schloegl <a.schloegl@ieee.org>
 %
 %  Time Series Analysis - a toolbox for the use with Matlab
@@ -26,13 +28,16 @@
 %   ar_spa	(*) spectral analysis based on the autoregressive model
 %   detrend 	(*) removes trend, can handle missing values, non-equidistant sampled data       
 %   flix	floating index, interpolates data for non-interger indices
-%   quantiles   calculates quantiles 
+%
 %
-% Multivariate analysis (planned in future)
+% Multivariate analysis 
+%   adim	adaptive information matrix (inverse correlation matrix) 
 %   mvar	multivariate (vector) autoregressive estimation 
 %   mvfilter	multivariate filter
+%   mvfreqz	multivariate spectra 	
 %   arfit2	provides compatibility to ARFIT [Schneider and Neumaier, 2001]
-       	
+%
+%   	
 %  Conversions between Autocorrelation (AC), Autoregressive parameters (AR), 
 %             	prediction polynom (POLY) and Reflection coefficient (RC)  
 %   ac2poly 	(*) transforms autocorrelation into prediction polynom
diff -r 8fe7ccc46434 -r 9c2a30bbaa0b extra/tsa/mvfreqz.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/extra/tsa/mvfreqz.m	Tue Jun 24 19:22:37 2003 +0000
@@ -0,0 +1,202 @@
+function [S,h,PDC,COH,DTF,DC,pCOH,dDTF,ffDTF, pCOH2]=mvfreqz(B,A,C,N,Fs)
+% MVFREQZ multivariate frequency response
+% [S,h,PDC,COH,DTF,DC,pCOH,dDTF,ffDTF,pCOH2] = mvfreqz(B,A,C,N,Fs)
+%
+% INPUT: 
+% ======= 
+% A, B	multivariate polynomials defining the transfer function
+%
+%    a0*Y(n) = b0*X(n) + b1*X(n-1) + ... + bq*X(n-q)
+%                          - a1*Y(n-1) - ... - ap*Y(:,n-p)
+%
+%  A=[a0,a1,a2,...,ap] and B=[b0,b1,b2,...,bq] must be matrices of
+%  size  Mx((p+1)*M) and Mx((q+1)*M), respectively. 
+%
+%  X must be of size N*M
+%  a0,a1,...,ap, b0,b1,...,bq are matrices of size MxM
+%  a0 is usually the identity matrix I or must be invertible 
+%	
+%  C is the covariance of the input noise X
+%  N if scalar, N is the number of frequencies 
+%    if N is a vector, N are the designated frequencies. 
+%  Fs sampling rate [default 2*pi]
+% 
+%
+% OUTPUT: 
+% ======= 
+% S   	power spectrum
+% PDC 	partial directed coherence
+% DC  	directed coupling	
+% COH 	coherence
+% DTF 	directed transfer function
+% pCOH 	partial coherence
+% dDTF 	direct Directed Transfer function
+% ffDTF full frequency Directed Transfer Function 
+% pCOH2 partial coherence -alternative method 
+%
+%
+% see also: FREQZ, MVFILTER, MVAR
+%
+% 
+% REFERENCE(S):
+% H. Liang et al. Neurocomputing, 32-33, pp.891-896, 2000. 
+% L.A. Baccala and K. Samashima, Biol. Cybern. 84,463-474, 2001. 
+% A. Korzeniewska, et al. Journal of Neuroscience Methods, 125, 195-207, 2003. 
+% Piotr J. Franaszczuk, Ph.D. and Gregory K. Bergey, M.D.
+% 	Fast Algorithm for Computation of Partial Coherences From Vector Autoregressive Model Coefficients
+%	World Congress 2000, Chicago. 
+
+%	$Revision$
+%	$Id$
+%	Copyright (C) 1996-2003 by Alois Schloegl <a.schloegl@ieee.org>	
+
+% This library is free software; you can redistribute it and/or
+% modify it under the terms of the GNU Library General Public
+% License as published by the Free Software Foundation; either
+% Version 2 of the License, or (at your option) any later version.
+%
+% This library is distributed in the hope that it will be useful,
+% but WITHOUT ANY WARRANTY; without even the implied warranty of
+% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+% Library General Public License for more details.
+%
+% You should have received a copy of the GNU Library General Public
+% License along with this library; if not, write to the
+% Free Software Foundation, Inc., 59 Temple Place - Suite 330,
+% Boston, MA  02111-1307, USA.
+
+[K1,K2] = size(A);
+p = K2/K1-1;
+%a=ones(1,p+1);
+[K1,K2] = size(B);
+q = K2/K1-1;
+%b=ones(1,q+1);
+if nargin<3
+        C = ones(K1,K1);
+end;
+if nargin<4,
+        N = 512;
+end;
+if nargin<5,
+	Fs= 1;        
+end;
+if all(size(N)==1),	
+        f = (0:N-1)/N;
+else
+        f = N;
+        N = length(N);
+end;
+s = exp(i*2*pi*f/Fs);
+z = i*2*pi/Fs;
+
+h=zeros(K1,K1,N);
+g=zeros(K1,K1,N);
+SP=zeros(K1,K1,N);
+DTF=zeros(K1,K1,N);
+COH=zeros(K1,K1,N);
+%COH2=zeros(K1,K1,N);
+PDC=zeros(K1,K1,N);
+PDCF = zeros(K1,K1,N);
+pCOH = zeros(K1,K1,N);
+invC=inv(C);
+tmp1=zeros(1,K1);
+tmp2=zeros(1,K1);
+
+M = zeros(K1,K1,N);
+detG = zeros(N,1);
+
+for n=1:N,
+        atmp = zeros(K1);
+        for k = 1:p+1,
+                atmp = atmp + A(:,k*K1+(1-K1:0))*exp(z*(k-1)*f(n));
+        end;        
+        btmp = zeros(K1);
+        for k = 1:q+1,
+                btmp = btmp + B(:,k*K1+(1-K1:0))*exp(z*(k-1)*f(n));
+        end;        
+        h(:,:,n) = atmp\btmp;        
+        S(:,:,n) = h(:,:,n)*C*h(:,:,n)';        
+        
+        for k1 = 1:K1,
+                tmp = squeeze(atmp(:,k1));
+                tmp1(k1) = sqrt(tmp'*tmp);
+                tmp2(k1) = sqrt(tmp'*invC*tmp);
+        end;
+        
+        %PDCF(:,:,n,kk) = abs(atmp)./tmp2(ones(1,K1),:);
+        %PDC(:,:,n,kk)  = abs(atmp)./tmp1(ones(1,K1),:);
+        PDCF(:,:,n) = abs(atmp)./tmp2(ones(1,K1),:);
+        PDC(:,:,n)  = abs(atmp)./tmp1(ones(1,K1),:);
+        
+        g = atmp/btmp;        
+        G(:,:,n) = g'*invC*g;
+        detG(n) = det(G(:,:,n));        
+        
+end;
+
+size(detG)
+size(squeeze(M(2,1,:)))
+size(G(2,1,:))
+
+if nargout<4, return; end;
+        
+%%%%% directed transfer function
+for k1=1:K1;
+        DEN=sum(abs(h(k1,:,:)).^2,2);	        
+for k2=1:K2;
+        %COH2(k1,k2,:) = abs(S(k1,k2,:).^2)./(abs(S(k1,k1,:).*S(k2,k2,:)));
+        COH(k1,k2,:) = abs(S(k1,k2,:))./sqrt(abs(S(k1,k1,:).*S(k2,k2,:)));
+	%DTF(k1,k2,:) = sqrt(abs(h(k1,k2,:).^2))./DEN;	        
+	DTF(k1,k2,:) = abs(h(k1,k2,:))./sqrt(DEN);
+        ffDTF(k1,k2,:) = abs(h(k1,k2,:))./sqrt(sum(DEN,3));
+        pCOH2(k1,k2,:) = abs(G(k1,k2,:).^2)./(G(k1,k1,:).*G(k2,k2,:));
+        
+        M(k2,k1,:) = ((-1)^(k1+k2))*squeeze(G(k1,k2,:))./detG; % oder ist M = G?
+end;
+end;
+
+for k1=1:K1;
+for k2=1:K2;
+        pCOH(k1,k2,:) = abs(M(k1,k2,:).^2)./(M(k1,k1,:).*M(k2,k2,:));
+end;
+end;
+
+dDTF = pCOH2.*ffDTF; 
+
+
+
+if nargout<6, return; end;
+        
+DC = zeros(K1);
+for k = 1:p,
+        DC = DC + A(:,k*K1+(1:K1)).^2;
+end;        
+
+
+
+return;
+
+if nargout<7, return; end;
+
+for n=1:N,
+        %COH2(k1,k2,:) = abs(S(k1,k2,:).^2)./(abs(S(k1,k1,:).*S(k2,k2,:)));
+        M(k1,k2,n) = det(squeeze(S([1:k1-1,k1+1:K1],[1:k2-1,k2+1:K2],n)));
+        
+end;
+
+for k1=1:K1;
+for k2=1:K2;
+for n=1:N,
+        %COH2(k1,k2,:) = abs(S(k1,k2,:).^2)./(abs(S(k1,k1,:).*S(k2,k2,:)));
+         M(k1,k2,n) = det(squeeze(S([1:k1-1,k1+1:K1],[1:k2-1,k2+1:K2],n)));
+        
+end;
+end;
+end;
+
+for k1=1:K1;
+for k2=1:K2;
+        pCOH(k1,k2,:) = abs(M(k1,k2,:).^2)./(M(k1,k1,:).*M(k2,k2,:));
+end;
+end;
+