-
Notifications
You must be signed in to change notification settings - Fork 714
/
ft_connectivity_ppc.m
100 lines (95 loc) · 4.25 KB
/
ft_connectivity_ppc.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
function [c, v, n] = ft_connectivity_ppc(inputdata, varargin)
% FT_CONNECTIVITY_PPC computes pairwise phase consistency or weighted pairwise phase
% consistency from a data-matrix containing a cross-spectral density. This implements
% the method described in Vinck M, van Wingerden M, Womelsdorf T, Fries P, Pennartz
% CM. The pairwise phase consistency: a bias-free measure of rhythmic neuronal
% synchronization. Neuroimage. 2010 May 15;51(1):112-22.
%
% Use as
% [c, v, n] = ft_connectivity_ppc(inputdata, ...)
%
% Where the input data input should be organized as:
% Repetitions x Channel x Channel (x Frequency) (x Time)
% or
% Repetitions x Channelcombination (x Frequency) (x Time)
%
% The first dimension should contain repetitions and should not contain an average
% already. Also, it should not consist of leave-one-out averages.
%
% The output c contains the ppc, v is a leave-one-out variance estimate which is only
% computed if dojack = 1,and n is the number of repetitions in the input data.
%
% Additional optional input arguments come as key-value pairs:
% 'dojack' = boolean specifying whether the repetitions represent leave-one-out samples
% 'weighted' = boolean, whether to compute unweighted ppc or weighted ppc, the weighting
% is according to the magnitude of the cross-spectrum
% 'feedback' = 'none', 'text', 'textbar', 'dial', 'etf', 'gui' type of feedback showing progress of computation, see FT_PROGRESS
%
% See also CONNECTIVITY, FT_CONNECTIVITYANALYSIS
% Copyright (C) 2011, Martin Vinck
%
% This file is part of FieldTrip, see http:https://www.fieldtriptoolbox.org
% for the documentation and details.
%
% FieldTrip 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.
%
% FieldTrip 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 FieldTrip. If not, see <http:https://www.gnu.org/licenses/>.
%
% $Id$
feedback = ft_getopt(varargin, 'feedback', 'none');
weighted = ft_getopt(varargin, 'weighted');
dojack = ft_getopt(varargin, 'dojack', false);
siz = size(inputdata);
n = siz(1);
ft_progress('init', feedback, 'computing metric...');
if n>1
if ~weighted
inputdata = inputdata./abs(inputdata); % normalize the crosspectrum
outsum = nansum(inputdata); % compute the sum; this is 1 x size(2:end)
c = (outsum.*conj(outsum) - n)./(n*(n-1)); % do the pairwise thing in a handy way
else
outsum = nansum(inputdata); % normalization of the WPLI
outssq = nansum(inputdata.*conj(inputdata));
outsumw = nansum(abs(inputdata));
c = (outsum.*conj(outsum) - outssq)./(outsumw.*conj(outsumw) - outssq); % do the pairwise thing in a handy way
end
c = reshape(c,siz(2:end)); % remove the first singular dimension
else
c = NaN(siz(2:end)); % for one observation, we should return NaNs
ft_warning('computation ppc requires >1 trial, returning NaNs')
end
[leave1outsum, leave1outssq] = deal(0);
if dojack && n>2 % n needs to be larger than 2 to get a meaningful variance
for k = 1:n
% this code works with both formats of input, also if it is 5-D
s = outsum - inputdata(k,:,:,:,:,:,:); % index up to 7-D, this also works for 5-D then
if ~weighted
num = s.*conj(s) - (n-2);
denom = (n-1)*(n-2);
else
sq = outssq - inputdata(k,:,:,:,:,:,:).*conj(inputdata(k,:,:,:,:,:,:));
sw = outsumw - abs(inputdata(k,:,:,:,:,:,:));
num = s.*conj(s) - sq;
denom = sw.*conj(sw) - sq;
end
leave1outsum = leave1outsum + num./denom;
leave1outssq = leave1outssq + (num./denom).^2;
end
% compute the sem here
v = (n-1).^2*(leave1outssq - (leave1outsum.^2)./n)./(n - 1); % 11.5 efron, sqrt and 1/n done in FT_CONNECTIVITYANALYSIS
v = reshape(v,siz(2:end)); % remove the first singular dimension
elseif dojack && n<=2
v = NaN(siz(2:end));
else
v = [];
end
ft_progress('close');