%% data from 'F:\oxytocin\OXYT011\OXYT011_005_220224_160927',
basepath = pwd;
load('OXYT011_005_220224_160927.eegstates.mat');
info = read_Intan_RHD2000_file(basepath);
analogin = read_Intan_analogin_file('analogin.dat',info.num_board_adc_channels);
lfp_fbtime = find(analogin(1,:)>0.3,1);

rippleCh = 71;
THLFPs = bz_GetLFP(rippleCh,'basepath',basepath,'noPrompts',true);
thetarange = linspace(0.5,250,1000);
lfp = double(THLFPs.data);
lfp(abs(lfp)>8000) = nan;
lfp = fillmissing(lfp,'linear');
[ss,ft,tt] = spectrogram(lfp,2*THLFPs.samplingRate,0.6*THLFPs.samplingRate,thetarange,THLFPs.samplingRate);
[spectrogrammt,tm,fm] = MTSpectrogram(lfp,'window',2);
bands = SpectrogramBands(spectrogrammt,fm);

figure;
colormap('jet')
imagesc(tt,thetarange,abs(ss));clim([0,1e5]);
yticks([0.5,20,50,100,200]); ylabel('Frequency (Hz)'); xlabel('Time (s)');
set(gca,'YDir','normal');
hold on;plot(tm,zscore(smoothdata(bands.spindles,'gaussian',20))*10+20,'color','w','LineWidth',2)

figure
tr = [7600 8100]
subplot(5,1,1); plot(tm,zscore(smoothdata(bands.spindles,'gaussian',20)),'color','r');xlim([tr(1) tr(2)]);
%subplot(5,1,2); plot([1:length(dFoverF)]/100+lfp_fbtime(1)/20000, zscore(smoothdata(dFoverF,'gaussian',20)),'g');xlim([tr(1) tr(2)])
subplot(5,1,3); plot(StateInfo.fspec{1,1}.to, zscore(smoothdata(StateInfo.motion,'gaussian',5)),'b');xlim([tr(1) tr(2)]);
subplot(5,1,4);plot(1:length(SleepState.idx.states*10),SleepState.idx.states);xlim([tr(1) tr(2)]);
subplot(5,1,5);plot(1:length(nREMnoArousals),nREMnoArousals);xlim([tr(1) tr(2)]);

%% OXT spectogram
figure;
tr = [7600 8100];
h = pcolor(tt,log2(thetarange),imgaussfilt(abs(ss),'FilterSize',3));
set(h, 'EdgeColor', 'none');shading interp;
ylabel('Frequency (Hz)'); xlabel('Time (s)');
set(gca,'YDir','normal');
xlim([tr(1) tr(2)]);ylim([log2(4) log2(64)])
clim([-1.5e5,(3e5)]); 
colormap(flipud(hsv));
xt = get(gca, 'YTick');
set (gca, 'YTickLabel', 2.^xt);
box off
set(gcf,'Position',[100 100 900 300])
exportgraphics(gcf,'Fig2A_OXT_Spectrogram.pdf','BackgroundColor','none','ContentType','vector')



%% data from 'F:\LIMP019\LIMP019_004_220303_115500',
%cd F:\LIMP019\LIMP019_004_220303_115500
basepath = pwd;
info = read_Intan_RHD2000_file(basepath);
analogin = read_Intan_analogin_file('analogin.dat',info.num_board_adc_channels);
lfp_fbtime = find(analogin(1,:)>0.3,1);

rippleCh = 21;
THLFPs = bz_GetLFP(rippleCh,'basepath',basepath,'noPrompts',true);
thetarange = linspace(0.5,250,1000);
lfp = double(THLFPs.data);
lfp(abs(lfp)>8000) = nan;
lfp = fillmissing(lfp,'linear');
[ss,ft,tt] = spectrogram(lfp,2*THLFPs.samplingRate,0.6*THLFPs.samplingRate,thetarange,THLFPs.samplingRate);
[spectrogrammt,tm,fm] = MTSpectrogram(lfp,'window',2);
bands = SpectrogramBands(spectrogrammt,fm);

figure;
colormap('jet')
imagesc(tt,thetarange,abs(ss));clim([0,8e4]);
yticks([0.5,20,50,100,200]); ylabel('Frequency (Hz)'); xlabel('Time (s)');
set(gca,'YDir','normal');

figure
tr = [4300 4800]
subplot(5,1,1); plot(tm,zscore(smoothdata(bands.spindles,'gaussian',20)),'color','r');xlim([tr(1) tr(2)]);
subplot(5,1,2); plot([1:length(dFoverF)]/100+lfp_fbtime(1)/20000, zscore(smoothdata(dFoverF,'gaussian',20)),'g');xlim([tr(1) tr(2)])
subplot(5,1,3); plot(accel.timestamps/2, zscore(smoothdata(accel.acceleration.motion,'gaussian',5)),'b');xlim([tr(1) tr(2)])
subplot(5,1,4);plot(1:length(SleepState.idx.states*10),SleepState.idx.states);xlim([tr(1) tr(2)]);
subplot(5,1,5);plot(1:length(nREMnoArousals),nREMnoArousals);xlim([tr(1) tr(2)]);

%%
figure;
tr = [4300 4800];
h = pcolor(tt,log2(thetarange),abs(ss));
set(h, 'EdgeColor', 'none');shading interp;
ylabel('Frequency (Hz)'); xlabel('Time (s)');
set(gca,'YDir','normal');
xlim([tr(1) tr(2)]);ylim([log2(4) log2(64)])
clim([-1.5e5,(3e5)])
colormap(flipud(hsv));
xt = get(gca, 'YTick');
set (gca, 'YTickLabel', 2.^xt);
box off
set(gcf,'Position',[100 100 900 300])
exportgraphics(gcf,'Fig2A_ACH_Spectrogram.pdf','BackgroundColor','none','ContentType','vector')
