This shows you the differences between two versions of the page.
analysis:amplipex:phaseplot [2014/01/21 17:34] mbonsma created |
analysis:amplipex:phaseplot [2018/07/07 10:19] |
||
---|---|---|---|
Line 1: | Line 1: | ||
- | <code matlab> | ||
- | %% channels specified by Eric | ||
- | |||
- | |||
- | cd('C:\Users\mvdmlab\Documents\Data\promoted\R036\R036-2013-05-23'); | ||
- | fname = 'R036-2013-05-23-pre.dat'; %this is some of ERic's data | ||
- | |||
- | %% | ||
- | channels = [33,57,1,25,36,58,2,26,38,62,6,30,37,61,5,29]; % specify channels in the APX data file to load | ||
- | %% | ||
- | data = AMPX_loadData(fname,channels_to_load,2); % note decimation factor of 1 | ||
- | |||
- | %% frequency ranges of interest | ||
- | |||
- | % low gamma: 45-65 Hz | ||
- | lgamma = [45 65]; | ||
- | stopbandgamm = 15; | ||
- | Fs = data.hdr.Fs; | ||
- | |||
- | %% filter for low gamma | ||
- | Wp2 = lgamma * 2 / Fs; | ||
- | Ws2 = [lgamma(1)-stopbandgamm lgamma(2)+stopbandgamm] * 2 / Fs; | ||
- | [N2,Wn2] = cheb1ord( Wp2, Ws2, 3, 20); % use a chebyshev because we want a steep cutoff and can sacrifice the passband quality | ||
- | [b2,a2] = cheby1(N2,2,Wn2); | ||
- | |||
- | [N2b,Wn2b] = buttord( Wp2, Ws2, 3, 20); % try a butterworth filter too. | ||
- | [b2b,a2b] = butter(N2b,Wn2b); | ||
- | |||
- | fvtool(b2,a2,b2b,a2b) | ||
- | |||
- | %use the butterworth here for a cleaner passband | ||
- | %% filter data | ||
- | ygamm = cell(1,length(channels)); | ||
- | |||
- | %% | ||
- | for ii = 1:length(channels); | ||
- | ygamm{ii} = filtfilt(b2b,a2b,data.channels{ii}); | ||
- | end | ||
- | |||
- | %% plot to check (optional) | ||
- | figure; | ||
- | plot(data.tvec,data.channels{1},'b'); | ||
- | hold on; | ||
- | plot(data.tvec,ygamm{1},'r'); | ||
- | |||
- | plot(data.tvec,data.channels{2},'m'); | ||
- | plot(data.tvec,data.channels{3},'g'); | ||
- | plot(data.tvec,data.channels{4},'k'); | ||
- | set(gca,'XLim',[124.95 125]); | ||
- | |||
- | |||
- | %% get phase | ||
- | phi = cell(1,length(channels)); | ||
- | for ii = 1:length(channels); | ||
- | phi{ii} = angle(hilbert(ygamm{ii})); | ||
- | end | ||
- | |||
- | %% make nice fancy plot - this one plots filtered signal and phase | ||
- | xwin = [124.8 125]; | ||
- | |||
- | figure; | ||
- | [AX,~,~] = plotyy(data.tvec,(ygamm{1}-mean(ygamm{1})+1*2),data.tvec,(phi{1}-mean(phi{1})+1*20)); | ||
- | freq = AX(1); | ||
- | ph = AX(2); | ||
- | hold all; | ||
- | set(AX,'XLim',xwin); | ||
- | hold(AX(1),'on'); hold(AX(2),'on'); | ||
- | |||
- | for ii = 2:length(channels); | ||
- | plot(freq,data.tvec,ygamm{ii}-mean(ygamm{ii})+ii*2); | ||
- | plot(ph,data.tvec,phi{ii}-mean(phi{ii})+ii*20); | ||
- | hold(freq,'on'); hold(ph,'on'); | ||
- | count = ii; | ||
- | end | ||
- | set(AX(1),'YLim',[0,count*2+2]); | ||
- | set(AX(2),'YLim',[0,count*20+20]); | ||
- | |||
- | %% make data for colormap testing | ||
- | i1 = find(data.tvec >= xwin(1)); | ||
- | |||
- | %% | ||
- | y = cell(1,length(channels)); | ||
- | cdata = cell(1,length(channels)); | ||
- | |||
- | x = data.tvec(i1(1):i1(1000)); | ||
- | |||
- | for ii = 1:length(channels) | ||
- | y{ii} = ygamm{ii}(i1(1):i1(1000)); | ||
- | cdata{ii} = phi{ii}(i1(1):i1(1000)); | ||
- | end | ||
- | |||
- | %% plot filtered signal with phase as a colour | ||
- | |||
- | scalefactor = 2; | ||
- | |||
- | figure; | ||
- | colormap(hsv); | ||
- | |||
- | for ii = 1:length(channels); | ||
- | | ||
- | scatter(x,y{ii}-mean(y{ii})+ii*scalefactor,10,cdata{ii},'fill'); | ||
- | | ||
- | hold on; | ||
- | count = ii; | ||
- | end | ||
- | ytick = [count:scalefactor:count*scalefactor]; | ||
- | set(gca,'YLim',[0,count*scalefactor+scalefactor],'YTick',ytick,'YTickLabel',channels); | ||
- | |||
- | |||
- | </code> |