1 function [S,h,PDC,COH,DTF,DC,pCOH,dDTF,ffDTF, pCOH2, PDCF, coh,GGC,Af,GPDC,GGC2]=mvfreqz(B,A,C,N,Fs)
2 % MVFREQZ multivariate frequency response
3 % [S,h,PDC,COH,DTF,DC,pCOH,dDTF,ffDTF,pCOH2,PDCF,coh,GGC,Af,GPDC] = mvfreqz(B,A,C,f,Fs)
4 % [...] = mvfreqz(B,A,C,N,Fs)
8 % A, B multivariate polynomials defining the transfer function
10 % a0*Y(n) = b0*X(n) + b1*X(n-1) + ... + bq*X(n-q)
11 % - a1*Y(n-1) - ... - ap*Y(:,n-p)
13 % A=[a0,a1,a2,...,ap] and B=[b0,b1,b2,...,bq] must be matrices of
14 % size Mx((p+1)*M) and Mx((q+1)*M), respectively.
16 % C is the covariance of the input noise X (i.e. D'*D if D is the mixing matrix)
17 % N if scalar, N is the number of frequencies
18 % if N is a vector, N are the designated frequencies.
19 % Fs sampling rate [default 2*pi]
21 % A,B,C and D can by obtained from a multivariate time series
22 % through the following commands:
23 % [AR,RC,PE] = mvar(Y,P);
24 % M = size(AR,1); % number of channels
27 % C = PE(:,M*P+1:M*(P+1));
29 % Fs sampling rate in [Hz]
30 % (N number of frequencies for computing the spectrum, this will become OBSOLETE),
31 % f vector of frequencies (in [Hz])
37 % h transfer functions, abs(h.^2) is the non-normalized DTF [11]
38 % PDC partial directed coherence [2]
39 % DC directed coupling
40 % COH coherency (complex coherence) [5]
41 % DTF directed transfer function
42 % pCOH partial coherence
43 % dDTF direct Directed Transfer function
44 % ffDTF full frequency Directed Transfer Function
45 % pCOH2 partial coherence - alternative method
46 % GGC a modified version of Geweke's Granger Causality [Geweke 1982]
47 % !!! it uses a Multivariate AR model, and computes the bivariate GGC as in [Bressler et al 2007].
48 % This is not the same as using bivariate AR models and GGC as in [Bressler et al 2007]
49 % Af Frequency transform of A(z), abs(Af.^2) is the non-normalized PDC [11]
50 % PDCF Partial Directed Coherence Factor [2]
51 % GPDC Generalized Partial Directed Coherence [9,10]
53 % see also: FREQZ, MVFILTER, MVAR
56 % [1] H. Liang et al. Neurocomputing, 32-33, pp.891-896, 2000.
57 % [2] L.A. Baccala and K. Samashima, Biol. Cybern. 84,463-474, 2001.
58 % [3] A. Korzeniewska, et al. Journal of Neuroscience Methods, 125, 195-207, 2003.
59 % [4] Piotr J. Franaszczuk, Ph.D. and Gregory K. Bergey, M.D.
60 % Fast Algorithm for Computation of Partial Coherences From Vector Autoregressive Model Coefficients
61 % World Congress 2000, Chicago.
62 % [5] Nolte G, Bai O, Wheaton L, Mari Z, Vorbach S, Hallett M.
63 % Identifying true brain interaction from EEG data using the imaginary part of coherency.
64 % Clin Neurophysiol. 2004 Oct;115(10):2292-307.
65 % [6] Schlogl A., Supp G.
66 % Analyzing event-related EEG data with multivariate autoregressive parameters.
67 % (Eds.) C. Neuper and W. Klimesch,
68 % Progress in Brain Research: Event-related Dynamics of Brain Oscillations.
69 % Analysis of dynamics of brain oscillations: methodological advances. Elsevier.
70 % Progress in Brain Research 159, 2006, p. 135 - 147
71 % [7] Bressler S.L., Richter C.G., Chen Y., Ding M. (2007)
72 % Cortical fuctional network organization from autoregressive modelling of loal field potential oscillations.
73 % Statistics in Medicine, doi: 10.1002/sim.2935
75 % J.Am.Stat.Assoc., 77, 304-313.
76 % [9] L.A. Baccala, D.Y. Takahashi, K. Sameshima. (2006)
77 % Generalized Partial Directed Coherence.
78 % Submitted to XVI Congresso Brasileiro de Automatica, Salvador, Bahia.
79 % [10] L.A. Baccala, D.Y. Takahashi, K. Sameshima.
80 % Computer Intensive Testing for the Influence Between Time Series,
81 % Eds. B. Schelter, M. Winterhalder, J. Timmer:
82 % Handbook of Time Series Analysis - Recent Theoretical Developments and Applications
85 % On the evaluation of informatino flow in multivariate systems by the directed transfer function
86 % Biol. Cybern. 94: 469-482, 2006.
88 % $Id: mvfreqz.m 8141 2011-03-02 08:01:58Z schloegl $
89 % Copyright (C) 1996-2008 by Alois Schloegl <a.schloegl@ieee.org>
90 % This is part of the TSA-toolbox. See also
91 % http://hci.tugraz.at/schloegl/matlab/tsa/
92 % http://octave.sourceforge.net/
93 % http://biosig.sourceforge.net/
95 % This program is free software: you can redistribute it and/or modify
96 % it under the terms of the GNU General Public License as published by
97 % the Free Software Foundation, either version 3 of the License, or
98 % (at your option) any later version.
100 % This program is distributed in the hope that it will be useful,
101 % but WITHOUT ANY WARRANTY; without even the implied warranty of
102 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
103 % GNU General Public License for more details.
105 % You should have received a copy of the GNU General Public License
106 % along with this program. If not, see <http://www.gnu.org/licenses/>.
122 f = (0:N-1)*(Fs/(2*N));
125 fprintf(1,'Warning MVFREQZ: The forth input argument N is a scalar, this is ambigous.\n');
126 fprintf(1,' In the past, N was used to indicate the number of spectral lines. This might change.\n');
127 fprintf(1,' In future versions, it will indicate the spectral line.\n');
128 f = (0:N-1)*(Fs/(2*N));
133 s = exp(i*2*pi*f/Fs);
143 %COH2=zeros(K1,K1,N);
145 %PDC3=zeros(K1,K1,N);
146 PDCF = zeros(K1,K1,N);
147 pCOH = zeros(K1,K1,N);
159 ddc2 = diag(diag(C).^(-1/2));
163 atmp = atmp + A(:,k*K1+(1-K1:0))*exp(z*(k-1)*f(n));
166 % compensation of instantaneous correlation
171 btmp = btmp + B(:,k*K1+(1-K1:0))*exp(z*(k-1)*f(n));
173 h(:,:,n) = atmp\btmp;
174 Af(:,:,n) = atmp/btmp;
175 S(:,:,n) = h(:,:,n)*C*h(:,:,n)'/Fs;
176 S1(:,:,n) = h(:,:,n)*h(:,:,n)';
177 ctmp = ddc2*atmp; %% used for GPDC
179 tmp = squeeze(atmp(:,k1));
180 tmp1(k1) = sqrt(tmp'*tmp);
181 tmp2(k1) = sqrt(tmp'*invC*tmp);
183 %tmp = squeeze(atmp(k1,:)');
184 %tmp3(k1) = sqrt(tmp'*tmp);
186 tmp = squeeze(ctmp(:,k1));
187 tmp3(k1) = sqrt(tmp'*tmp);
190 PDCF(:,:,n) = abs(atmp)./tmp2(ones(1,K1),:);
191 PDC(:,:,n) = abs(atmp)./tmp1(ones(1,K1),:);
192 GPDC(:,:,n) = abs(ctmp)./tmp3(ones(1,K1),:);
193 %PDC3(:,:,n) = abs(atmp)./tmp3(:,ones(1,K1));
196 G(:,:,n) = g'*invC*g;
197 detG(n) = det(G(:,:,n));
200 if nargout<4, return; end;
202 %%%%% directed transfer function
204 DEN=sum(abs(h(k1,:,:)).^2,2);
206 %COH2(k1,k2,:) = abs(S(k1,k2,:).^2)./(abs(S(k1,k1,:).*S(k2,k2,:)));
207 COH(k1,k2,:) = (S(k1,k2,:))./sqrt(abs(S(k1,k1,:).*S(k2,k2,:)));
208 coh(k1,k2,:) = (S1(k1,k2,:))./sqrt(abs(S1(k1,k1,:).*S1(k2,k2,:)));
209 %DTF(k1,k2,:) = sqrt(abs(h(k1,k2,:).^2))./DEN;
210 DTF(k1,k2,:) = abs(h(k1,k2,:))./sqrt(DEN);
211 ffDTF(k1,k2,:) = abs(h(k1,k2,:))./sqrt(sum(DEN,3));
212 pCOH2(k1,k2,:) = abs(G(k1,k2,:).^2)./(G(k1,k1,:).*G(k2,k2,:));
214 %M(k2,k1,:) = ((-1)^(k1+k2))*squeeze(G(k1,k2,:))./detG; % oder ist M = G?
220 if nargout<6, return; end;
224 DC = DC + A(:,k*K1+(1:K1)).^2;
228 if nargout<13, return; end;
232 % Bivariate Granger Causality (similar to Bressler et al. 2007. )
233 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,:));
234 %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,:));
240 if nargout<7, return; end;
244 M(k2,k1,:) = ((-1)^(k1+k2))*squeeze(G(k1,k2,:))./detG; % oder ist M = G?
251 pCOH(k1,k2,:) = abs(M(k1,k2,:).^2)./(M(k1,k1,:).*M(k2,k2,:));