wiki

===Phase Plot Script === Sandbox to load LFP data, filter data, compute and plot phases - MB Jan 21 2014 <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>