-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathcohere_bootstrap.m
96 lines (83 loc) · 2.79 KB
/
cohere_bootstrap.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
function [sl, slf, df, F] = cohere_bootstrap_signif_level( x, y, p, nb, varargin )
% COHERE_BOOTSTRAP_SIGNIF_LEVEL - simple bootstrap to estimate
% significance level coherence.
%
% Usage: [sig_lev, sig_lev_f, df, freq ] = cohere_bootstrap_signif_level( ...
% x, y, p, nb, VARARGIN )
%
% where x and y are your time series, p is the significance level
% (i.e. 0.95 for a 95% significance level), nb is the number of
% bootstraps to compute. p defaults to 0.95 and nb defaults to 100.
%
% VARARGIN are the extra arguments to be passed to the matlab cohere
% function. The cohere function is in the Signal Processing Toolbox, so you
% will need that for this function to work. Note that the matlab cohere
% function is in the process of being replaced by mscohere. For
% compatibility, I am using cohere here. If this function fails, that
% would be a first thing to check.
%
% If the last argument of varargin is not a string, then the string
% 'mean' is added on the end to detrend your data.
%
% Output arguments:
%
% sig_lev is the significance level. It is computed as the average over all
% frequencies. sig_lev_f gives the actual significance level computed for
% all frequencies. df is the approximate number of degrees of freedom that
% this configuration corresponds to. freq are the frequencies at which the
% coherence was calculated.
%
% This function works by randomizing your original time series to produce
% series that have no correlation and then computing the coherence. This
% is repeated to get a range of values for the coherence at a frequency.
% The p*100% percential is used to determine the significance level.
%
% If you don't have a time series and just want the value try
% x=y=rand(N,1), where N is the number of observations you want.
%
% This method appears to produce values slightly over theoretical value.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% $Id: cohere_bootstrap_signif_level.m,v 1.2 2005/01/08 02:49:29 dmk Exp $
%
% Copyright (C) 2004 David M. Kaplan
% Licence: GPL (Gnu Public License)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin < 3
p = 0.95;
end
if nargin < 4
nb = 100;
end
if length(varargin) == 0 | ~isa(varargin{end},'char')
varargin{end+1} = 'mean';
end
s = length(y);
I = unidrnd( s, s, 1 );
[C,F] = mscohere( x, y(I), varargin{:} );
CC=[ C(:)'; zeros(nb-1,length(C)) ];
for l = 1:nb-1
I = unidrnd( s, s, 1 );
C = mscohere( x, y(I), varargin{:} );
CC(l+1,:) = C(:)';
end
[n,xx]=hist(CC,1000);
n2 = cumsum(n);
slf = zeros(1,size(n2,2));
for k = 1:size(n2,2)
[nn,ii] = unique( n2(:,k) );
slf(k) = interp1( nn, xx(ii), nb*p, 'linear', 0 );
end
sl = mean(slf);
df = 2*s/length(F);
if nargout < 4
clear F
end
if nargout < 3
clear df
end
if nargout < 2
clear slf
end