octave-tsa 4.2.7-1build1 / usr / share / octave / packages / tsa-4.2.7 / mvfreqz.m

function [S,h,PDC,COH,DTF,DC,pCOH,dDTF,ffDTF, pCOH2, PDCF, coh,GGC,Af,GPDC,GGC2, DCOH]=mvfreqz(B,A,C,N,Fs)
% MVFREQZ multivariate frequency response
% [S,h,PDC,COH,DTF,DC,pCOH,dDTF,ffDTF,pCOH2,PDCF,coh,GGC,Af,GPDC,GGC2,DCOH] = mvfreqz(B,A,C,f,Fs)
% [...]  = 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. 
%
%  C is the covariance of the input noise X (i.e. D'*D if D is the mixing matrix)
%  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]
% 
%  A,B,C and D can by obtained from a multivariate time series 
%       through the following commands: 
%  [AR,RC,PE] = mvar(Y,P);
%       M = size(AR,1); % number of channels       
%       A = [eye(M),-AR];
%       B = eye(M); 
%       C = PE(:,M*P+1:M*(P+1)); 
%
% Fs 	sampling rate in [Hz]
% (N 	number of frequencies for computing the spectrum, this will become OBSOLETE),  
% f	vector of frequencies (in [Hz])  
%
%
% OUTPUT: 
% ======= 
% S   	power spectrum
% h	transfer functions, abs(h.^2) is the non-normalized DTF [11]
% PDC 	partial directed coherence [2]
% DC  	directed coupling [13] 	
% COH 	coherency (complex coherence) [5]
% DTF 	directed transfer function [3,13]
% pCOH 	partial coherence 
% dDTF 	direct Directed Transfer function
% ffDTF full frequency Directed Transfer Function 
% pCOH2 partial coherence - alternative method 
% GGC	a modified version of Geweke's Granger Causality [Geweke 1982]
%	   !!! it uses a Multivariate AR model, and computes the bivariate GGC as in [Bressler et al 2007]. 
%	   This is not the same as using bivariate AR models and GGC as in [Bressler et al 2007]
% Af	Frequency transform of A(z), abs(Af.^2) is the non-normalized PDC [11]
% PDCF 	Partial Directed Coherence Factor [2]
% GPDC 	Generalized Partial Directed Coherence [9,10]
% DCOH  directed coherence or Generalized DTF (GDTF) [12] (equ. 11a)
%
% see also: FREQZ, MVFILTER, MVAR
% 
% REFERENCE(S):
% [1] H. Liang et al. Neurocomputing, 32-33, pp.891-896, 2000. 
% [2] L.A. Baccala and K. Samashima, Biol. Cybern. 84,463-474, 2001. 
% [3] A. Korzeniewska, et al. Journal of Neuroscience Methods, 125, 195-207, 2003. 
% [4] 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. 
% [5] Nolte G, Bai O, Wheaton L, Mari Z, Vorbach S, Hallett M.
%	Identifying true brain interaction from EEG data using the imaginary part of coherency.
%	Clin Neurophysiol. 2004 Oct;115(10):2292-307. 
% [6] Schlogl A., Supp G.
%       Analyzing event-related EEG data with multivariate autoregressive parameters.
%       (Eds.) C. Neuper and W. Klimesch, 
%       Progress in Brain Research: Event-related Dynamics of Brain Oscillations. 
%       Analysis of dynamics of brain oscillations: methodological advances. Elsevier. 
%       Progress in Brain Research 159, 2006, p. 135 - 147
% [7] Bressler S.L., Richter C.G., Chen Y., Ding M. (2007)
%	Cortical fuctional network organization from autoregressive modelling of loal field potential oscillations.
%	Statistics in Medicine, doi: 10.1002/sim.2935 
% [8] Geweke J., 1982	
%	J.Am.Stat.Assoc., 77, 304-313.
% [9] L.A. Baccala, D.Y. Takahashi, K. Sameshima. (2006) 
% 	Generalized Partial Directed Coherence. 
%	Submitted to XVI Congresso Brasileiro de Automatica, Salvador, Bahia.  
% [10] L.A. Baccala, D.Y. Takahashi, K. Sameshima. 
% 	Computer Intensive Testing for the Influence Between Time Series, 
%	Eds. B. Schelter, M. Winterhalder, J. Timmer: 
%	Handbook of Time Series Analysis - Recent Theoretical Developments and Applications
%	Wiley, p.413, 2006.
% [11] M. Eichler
%	On the evaluation of informatino flow in multivariate systems by the directed transfer function
%	Biol. Cybern. 94: 469-482, 2006. 	
% [12] L. Faes, S. Erla, and G. Nollo, (2012)
%	Measuring Connectivity in Linear Multivariate Processes: Definitions, Interpretation, and Practical Analysis
%	Computational and Mathematical Methods in Medicine Volume 2012 (2012), Article ID 140513, 18 pages
%	doi:10.1155/2012/140513
% [13] Maciej Kaminski, Mingzhou Ding, Wilson A. Truccolo, Steven L. Bressler
%	Evaluating causal relations in neural systems: Granger causality, 
%       directed transfer function and statistical assessment of significance.
% 	Biol. Cybern. 85, 145-157 (2001)
%

%	$Id: mvfreqz.m 12366 2013-11-25 22:25:52Z schloegl $
%	Copyright (C) 1996-2008 by Alois Schloegl <alois.schloegl@gmail.com>	
%	Copyright (C) 2013 Martin Billinger	
%       This is part of the TSA-toolbox. See also 
%       http://pub.ist.ac.at/~schloegl/matlab/tsa/
%       http://octave.sourceforge.net/
%       http://biosig.sourceforge.net/
%
%    This program is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    This program 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 General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with this program.  If not, see <http://www.gnu.org/licenses/>.

[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 = eye(K1,K1);
end;
if nargin<5,
        Fs= 1;        
end;
if nargin<4,
        N = 512;
        f = (0:N-1)*(Fs/(2*N));
end;
if all(size(N)==1),	
	fprintf(1,'Warning MVFREQZ: The forth input argument N is a scalar, this is ambigous.\n');
	fprintf(1,'   In the past, N was used to indicate the number of spectral lines. This might change.\n');
	fprintf(1,'   In future versions, it will indicate the spectral line.\n');
        f = (0:N-1)*(Fs/(2*N));
else
        f = N;
end;
N = length(f);
s = exp(i*2*pi*f/Fs);
z = i*2*pi/Fs;

h=zeros(K1,K1,N);
Af=zeros(K1,K1,N);
g=zeros(K1,K1,N);
S=zeros(K1,K1,N);
S1=zeros(K1,K1,N);
DTF=zeros(K1,K1,N);
COH=zeros(K1,K1,N);
%COH2=zeros(K1,K1,N);
PDC=zeros(K1,K1,N);
%PDC3=zeros(K1,K1,N);
PDCF = zeros(K1,K1,N);
pCOH = zeros(K1,K1,N);
GGC=zeros(K1,K1,N);
GGC2=zeros(K1,K1,N);
DCOH=zeros(K1,K1,N);
invC=inv(C);
tmp1=zeros(1,K1);
tmp2=zeros(1,K1);

M = zeros(K1,K1,N);
detG = zeros(N,1);

%D = sqrtm(C);
%iD= inv(D);
ddc2 = diag(diag(C).^(-1/2));
ddc2i = diag(diag(C).^(1/2)); 
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;   
        
        % compensation of instantaneous correlation 
        % atmp = iD*atmp*D;
        
        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;        
        Af(:,:,n)  = atmp/btmp;        
        S(:,:,n)  = h(:,:,n)*C*h(:,:,n)'/Fs;        
        S1(:,:,n) = h(:,:,n)*h(:,:,n)';        
        ctmp = ddc2*atmp;	%% used for GPDC 
        dtmp = h(:,:,n) * ddc2i; %% used for directed coherence (DCOH)
        for k1 = 1:K1,
                tmp = squeeze(atmp(:,k1));
                tmp1(k1) = sqrt(tmp'*tmp);
                tmp2(k1) = sqrt(tmp'*invC*tmp);

                %tmp = squeeze(atmp(k1,:)');
                %tmp3(k1) = sqrt(tmp'*tmp);

                tmp = squeeze(ctmp(:,k1));
                tmp3(k1) = sqrt(tmp'*tmp);
                
                tmp = dtmp(k1,:);
                tmp4(k1) = sqrt(tmp*tmp');
        end;
        
        PDCF(:,:,n)  = abs(atmp)./tmp2(ones(1,K1),:);
        PDC(:,:,n)   = abs(atmp)./tmp1(ones(1,K1),:);
        GPDC(:,:,n)  = abs(ctmp)./tmp3(ones(1,K1),:);
        %PDC3(:,:,n) = abs(atmp)./tmp3(:,ones(1,K1));
        
        DCOH(:,:,n) = abs(dtmp) ./ tmp4(ones(1,K1),:)';
        
        g = atmp/btmp;        
        G(:,:,n) = g'*invC*g;
        detG(n) = det(G(:,:,n));        
end;

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,:)   = (S(k1,k2,:))./sqrt(abs(S(k1,k1,:).*S(k2,k2,:)));
                coh(k1,k2,:)   = (S1(k1,k2,:))./sqrt(abs(S1(k1,k1,:).*S1(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;

dDTF = pCOH2.*ffDTF; 

if nargout<6, return; end;

DC = zeros(K1);
for k = 1:p,
        DC = DC + A(:,k*K1+(1:K1)).^2;
end;        


if nargout<13, return; end;

for k1=1:K1;
        for k2=1:K2;
		% Bivariate Granger Causality (similar to Bressler et al. 2007. )
                GGC(k1,k2,:) = ((C(k1,k1)*C(k2,k2)-C(k1,k2)^2)/C(k2,k2))*real(h(k1,k2,:).*conj(h(k1,k2,:)))./abs(S(k2,k2,:));
                %GGC2(k1,k2,:) = -log(1-((C(k1,k1)*C(k2,k2)-C(k1,k2)^2)/C(k2,k2))*real(h(k1,k2,:).*conj(h(k1,k2,:)))./S(k2,k2,:));
        end;
end;

return;

if nargout<7, return; end;

for k1=1:K1;
        for k2=1: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;