[f1,f2] = butter(4,[55/1250/2 65/1250/2],'stop');
[f3,f4] = butter(4,[105/1250/2 115/1250/2],'stop');
rawEeg2{1} = filtfilt(f1,f2,double(rawEeg{1}));
rawEeg2{1} = filtfilt(f3,f4,(rawEeg2{1}));
FO.eeg{i} = rawEeg2{i};

FO.maxFreq = 180; %default starting frequency extent
FO.hanningW = 25; %default hanning smoothing window
colors = {};
FO.resolution = 0.25;
f = StateInfo.fspec;
f1 = round(sum(f{1}.fo >= 2 & f{1}.fo <= 4)/(2./FO.resolution));
spec = {};
for I = 1:FO.nCh
    fo = [];
    s = [];
    for i = 1:f1:(size(f{I}.spec, 2) - f1)
        s = [s, mean(f{I}.spec(:, i:(i + f1 - 1), :), 2)];
        fo = [fo; mean(f{I}.fo(i:(i + f1 - 1)))];
    end
    spec{I} = s;
end

for i = 1:length(f)
    n = prctile(reshape(f{i}.spec, 1, []), [1, 99]);
    f{i}.spec(f{i}.spec < n(1)) = n(1);
    f{i}.spec(f{i}.spec > n(2)) = n(2);
end

FO.unsmoothedSpec = spec;
FO.originalFO = f{1}.fo;

for i = 1:nCh
    FO.spec{i} = log10(convWithIn_misi(spec{i}, hanning(FO.hanningW)));
end

FO.fo = fo;
FO.to = f{1}.to;
FO.lims = [min(FO.to), max(FO.to)];

if nCh > 1
    min1 = min(min(FO.spec{1}(:, FO.fo <= FO.maxFreq)));
    max1 = max(max(FO.spec{1}(:, FO.fo <= FO.maxFreq)));
    for i = 2:nCh
        min2 = min(min(FO.spec{i}(:, FO.fo <= FO.maxFreq)));
        FO.spec{i} = FO.spec{i} - min2;
        max2 = max(max(FO.spec{i}(:, FO.fo <= FO.maxFreq)));
        FO.spec{i} = FO.spec{i}./max2;
        FO.spec{i} = FO.spec{i}.*(max1 - min1);
        FO.spec{i} = FO.spec{i} + min1;
    end
end
%%
Clim = [0.4 2.3];
figure
FO.iSpec{i} = imagesc(FO.to, FO.fo(FO.fo < 50), FO.spec{i}(:, FO.fo < 50)',Clim);
set(gca,'YDir','normal')
colormap('jet');
hold on
for ii = 1:size(SleepState.ints.REMstate,1)
    plot([SleepState.ints.REMstate(ii,1) SleepState.ints.REMstate(ii,2)],[45 45],'r-');
end
clear ii
for ii = 1:size(SleepState.ints.NREMstate,1)
    plot([SleepState.ints.NREMstate(ii,1) SleepState.ints.NREMstate(ii,2)],[46 46],'b-');
end
clear ii
for ii = 1:size(SleepState.ints.WAKEstate,1)
    plot([SleepState.ints.WAKEstate(ii,1) SleepState.ints.WAKEstate(ii,2)],[47 47],'k-');
end
clear ii
ylimit = get(gca,'ylim');
plot([1.216e4 1.516e4],[ylimit(2) ylimit(2)],'color','m','linewidth',3);
hold off
ylabel('Frequency (Hz)')
yyaxis right
plot(tempKell,'k');
box off
set(gca,'xtick',0:1800:2e4,'xticklabel',0:0.5:12);
xlabel('Time (hours)')
ylabel('Temperature (C)')

set(gcf,'renderer','Painters');
saveas(gcf,[basepath filesep basename '-spectrum-full-LOW.pdf'])
saveas(gcf,[basepath filesep basename '-spectrum-full-LOW.png'])
close all

%%
figure
FO.iSpec{i} = imagesc(FO.to, FO.fo(FO.fo > 70), FO.spec{i}(:, FO.fo > 70)',Clim);
set(gca,'YDir','normal')
colormap('jet');
ylim([70 150])
ylabel('Frequency (Hz)')
yyaxis right
plot(tempKell,'k');
box off
set(gca,'xtick',0:1800:2e4,'xticklabel',0:0.5:12);
xlabel('Time (hours)')
ylabel('Temperature (C)')

set(gcf,'renderer','Painters');
saveas(gcf,[basepath filesep basename '-spectrum-full-HIGH.pdf'])
saveas(gcf,[basepath filesep basename '-spectrum-full-HIGH.png'])
close all


%%

i = 1;
figure
FO.iSpec{i} = imagesc(FO.to, FO.fo(FO.fo < 50), FO.spec{i}(:, FO.fo < 50)');
set(gca,'YDir','normal')
ylabel('Frequency (Hz)')
colormap('jet');
hold on
for ii = 1:size(SleepState.ints.REMstate,1)
    plot([SleepState.ints.REMstate(ii,1) SleepState.ints.REMstate(ii,2)],[45 45],'r-');
end
clear ii
for ii = 1:size(SleepState.ints.NREMstate,1)
    plot([SleepState.ints.NREMstate(ii,1) SleepState.ints.NREMstate(ii,2)],[46 46],'b-');
end
clear ii
for ii = 1:size(SleepState.ints.WAKEstate,1)
    plot([SleepState.ints.WAKEstate(ii,1) SleepState.ints.WAKEstate(ii,2)],[47 47],'k-');
end
clear ii
hold off
yyaxis right
plot(tempKell,'k');
set(gca,'xtick',0:1800:2e4,'xticklabel',0:0.5:12);
xlabel('Time (hours)')
ylabel('Temperature (C)')
xlim([1.216e4 1.516e4])
box off
set(gcf,'renderer','Painters');
saveas(gcf,[basepath filesep basename '-spectrum-zoom-LOW.pdf'])
saveas(gcf,[basepath filesep basename '-spectrum-zoom-LOW.png'])
close all

%%
i = 1;
figure
FO.iSpec{i} = imagesc(FO.to, FO.fo(FO.fo > 70), FO.spec{i}(:, FO.fo > 70)');
set(gca,'YDir','normal')
colormap('jet');
ylabel('Frequency (Hz)')
ylim([70 150])
hold on
for ii = 1:size(SleepState.ints.REMstate,1)
    plot([SleepState.ints.REMstate(ii,1) SleepState.ints.REMstate(ii,2)],[45 45],'r-');
end
clear ii
for ii = 1:size(SleepState.ints.NREMstate,1)
    plot([SleepState.ints.NREMstate(ii,1) SleepState.ints.NREMstate(ii,2)],[46 46],'b-');
end
clear ii
for ii = 1:size(SleepState.ints.WAKEstate,1)
    plot([SleepState.ints.WAKEstate(ii,1) SleepState.ints.WAKEstate(ii,2)],[47 47],'k-');
end
clear ii
hold off
yyaxis right
plot(tempKell,'k');
set(gca,'xtick',0:1800:2e4,'xticklabel',0:0.5:12);
xlabel('Time (hours)')
ylabel('Temperature (C)')
xlim([1.216e4 1.516e4])
box off
set(gcf,'renderer','Painters');
saveas(gcf,[basepath filesep basename '-spectrum-zoom-HIGH.pdf'])
saveas(gcf,[basepath filesep basename '-spectrum-zoom-HIGH.png'])
close all





