Sandbox to load LFP data, filter data, compute and plot phases - MB Jan 21 2014

%% 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);