octave_packages/tsa-4.2.4/mvfreqz.m

   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)
   5 %
   6 % INPUT:
   7 % =======
   8 % A, B  multivariate polynomials defining the transfer function
   9 %
  10 %    a0*Y(n) = b0*X(n) + b1*X(n-1) + ... + bq*X(n-q)
  11 %                          - a1*Y(n-1) - ... - ap*Y(:,n-p)
  12 %
  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.
  15 %
  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]
  20 %
  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
  25 %       A = [eye(M),-AR];
  26 %       B = eye(M);
  27 %       C = PE(:,M*P+1:M*(P+1));
  28 %
  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])
  32 %
  33 %
  34 % OUTPUT:
  35 % =======
  36 % S     power spectrum
  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]
  52 %
  53 % see also: FREQZ, MVFILTER, MVAR
  54 %
  55 % REFERENCE(S):
  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
  74 % [8] Geweke J., 1982
  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
  83 %       Wiley, p.413, 2006.
  84 % [11] M. Eichler
  85 %       On the evaluation of informatino flow in multivariate systems by the directed transfer function
  86 %       Biol. Cybern. 94: 469-482, 2006.
  87
  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/
  94 %
  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.
  99 %
 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.
 104 %
 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/>.
 107
 108 [K1,K2] = size(A);
 109 p = K2/K1-1;
 110 %a=ones(1,p+1);
 111 [K1,K2] = size(B);
 112 q = K2/K1-1;
 113 %b=ones(1,q+1);
 114 if nargin<3
 115         C = eye(K1,K1);
 116 end;
 117 if nargin<5,
 118         Fs= 1;
 119 end;
 120 if nargin<4,
 121         N = 512;
 122         f = (0:N-1)*(Fs/(2*N));
 123 end;
 124 if all(size(N)==1),
 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));
 129 else
 130         f = N;
 131 end;
 132 N = length(f);
 133 s = exp(i*2*pi*f/Fs);
 134 z = i*2*pi/Fs;
 135
 136 h=zeros(K1,K1,N);
 137 Af=zeros(K1,K1,N);
 138 g=zeros(K1,K1,N);
 139 S=zeros(K1,K1,N);
 140 S1=zeros(K1,K1,N);
 141 DTF=zeros(K1,K1,N);
 142 COH=zeros(K1,K1,N);
 143 %COH2=zeros(K1,K1,N);
 144 PDC=zeros(K1,K1,N);
 145 %PDC3=zeros(K1,K1,N);
 146 PDCF = zeros(K1,K1,N);
 147 pCOH = zeros(K1,K1,N);
 148 GGC=zeros(K1,K1,N);
 149 GGC2=zeros(K1,K1,N);
 150 invC=inv(C);
 151 tmp1=zeros(1,K1);
 152 tmp2=zeros(1,K1);
 153
 154 M = zeros(K1,K1,N);
 155 detG = zeros(N,1);
 156
 157 %D = sqrtm(C);
 158 %iD= inv(D);
 159 ddc2 = diag(diag(C).^(-1/2));
 160 for n=1:N,
 161         atmp = zeros(K1);
 162         for k = 1:p+1,
 163                 atmp = atmp + A(:,k*K1+(1-K1:0))*exp(z*(k-1)*f(n));
 164         end;
 165
 166         % compensation of instantaneous correlation
 167         % atmp = iD*atmp*D;
 168
 169         btmp = zeros(K1);
 170         for k = 1:q+1,
 171                 btmp = btmp + B(:,k*K1+(1-K1:0))*exp(z*(k-1)*f(n));
 172         end;
 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
 178         for k1 = 1:K1,
 179                 tmp = squeeze(atmp(:,k1));
 180                 tmp1(k1) = sqrt(tmp'*tmp);
 181                 tmp2(k1) = sqrt(tmp'*invC*tmp);
 182
 183                 %tmp = squeeze(atmp(k1,:)');
 184                 %tmp3(k1) = sqrt(tmp'*tmp);
 185
 186                 tmp = squeeze(ctmp(:,k1));
 187                 tmp3(k1) = sqrt(tmp'*tmp);
 188         end;
 189
 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));
 194
 195         g = atmp/btmp;
 196         G(:,:,n) = g'*invC*g;
 197         detG(n) = det(G(:,:,n));
 198 end;
 199
 200 if nargout<4, return; end;
 201
 202 %%%%% directed transfer function
 203 for k1=1:K1;
 204         DEN=sum(abs(h(k1,:,:)).^2,2);
 205         for k2=1:K2;
 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,:));
 213
 214                 %M(k2,k1,:) = ((-1)^(k1+k2))*squeeze(G(k1,k2,:))./detG; % oder ist M = G?
 215         end;
 216 end;
 217
 218 dDTF = pCOH2.*ffDTF;
 219
 220 if nargout<6, return; end;
 221
 222 DC = zeros(K1);
 223 for k = 1:p,
 224         DC = DC + A(:,k*K1+(1:K1)).^2;
 225 end;
 226
 227
 228 if nargout<13, return; end;
 229
 230 for k1=1:K1;
 231         for k2=1:K2;
 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,:));
 235         end;
 236 end;
 237
 238 return;
 239
 240 if nargout<7, return; end;
 241
 242 for k1=1:K1;
 243         for k2=1:K2;
 244                 M(k2,k1,:) = ((-1)^(k1+k2))*squeeze(G(k1,k2,:))./detG; % oder ist M = G?
 245         end;
 246 end;
 247
 248
 249 for k1=1:K1;
 250         for k2=1:K2;
 251                 pCOH(k1,k2,:) = abs(M(k1,k2,:).^2)./(M(k1,k1,:).*M(k2,k2,:));
 252         end;
 253 end;
 254