-
Notifications
You must be signed in to change notification settings - Fork 2
/
nr_make_comodulogram_plv.m
executable file
·72 lines (63 loc) · 3.17 KB
/
nr_make_comodulogram_plv.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
function [Comodulogram] = make_comodulogram_plv(data,Fs,PhaseLow,PhaseHigh,AmpLow,AmpHigh,bad_times,Ampwidth)
%%used PLV method to calculate PAC. No statistics (but you can modify to do, use Tort script as a reference),
%will exclude noisy segments of data if you specify "badtimes"
%inputs:
%data= data (from 1 channel)you want to analyze as a row vector
%Fs = sampling rate
%PhaseLow = Low frequency value for phase in Hz (ex. 4).
%PhaseHigh = high value of phase frequency vector in Hz (ex. 45).
%bad_times = vector of data to exclude (b/c of artifacts etc), in units of sample points, leave empty if none. Note removing the bad
%segments of data is done after filtering to avoid edge artifacts.
%Ampwidth = the width of each band in the amplitude component of the PAC (the phase component is 1 by default)
%5 for amp and 1 for phase is what Canolty 2006 used (ex. 5)
%outputs
%Comodulogram = Comodulogram of PLV values for each phase and
%ampltide frequency examines (matrix phase x amp)
%initialize variables
PhaseFreqVector=[PhaseLow:PhaseHigh];
AmpFreqVector=[AmpLow:Ampwidth:AmpHigh];
Comodulogram = NaN(length(PhaseFreqVector),length(AmpFreqVector));
%loop through phase frequencies
phaseCount=1;
for p=PhaseLow:PhaseHigh
ampCount=1;
%filter data at phase frequencies with a bandwidth of 1.
filtered1=eegfilt(data,Fs,p,[]);
filtered1=eegfilt(filtered1,Fs,[],p+1);
%take angle of filtered data
slow = angle(hilbert(filtered1));
%loop through amplitude frequencies
for a=AmpLow:Ampwidth:AmpHigh
%take amplitude of filtered data for amplitude frequencies
fast = abs(hilbert(eegfilt(data,Fs,a,a+Ampwidth-1)));
%filter amplitude envelope at low frequency with bandwidth of 1
fast = eegfilt(fast,Fs,p,[]);
fast = eegfilt(fast,Fs,[],p+1);
%take angle of amplitude envelope filtered signal.
fast = angle(hilbert(fast));
x = slow;
y = fast;
%remove time points with noise
if (~isempty(bad_times))
x(bad_times)=[];
y(bad_times)=[];
end
%calculate the plv value as the difference in phase angle between low freq filtered signal and amplitude envelope
%filtered signal. Note that the difference is multipled by i and taken to the exponent b/c values are circular.
%then sum plv values over time points
plvs=exp(i*(y(1)-x(1)));
for ee=2:length(x)
plvs=plvs+exp(i*(y(ee)-x(ee)));
end
%take magnitide of plv (vector length), and divide by number of
%time points (i.e. average plv)
PLV=abs(plvs)/length(x);
%put in correct place in matrix
Comodulogram(phaseCount,ampCount)=PLV;
clear PLV
ampCount=ampCount+1;
end
% fprintf('\n*****\n\n\n\nfinished about %2.0f percent of electrode %d\n\n\n\n*****\n', phaseCount/gencnt*100, e)
phaseCount=phaseCount+1;
end
end