%%% PD STN Neurofeedback %%
%%% Shenghong He %%%
%%% University of Oxford %%%
%%% 04/11/2020 %%%
%%% shenghong.he@ndcn.ox.ac.uk
function NFB_analysis_all()
clc 
clear
global fileName;
fileName = {'P1L','P1R','P2L','P2R','P3L','P3R','P4L','P4R','P5L','P5R','P6L1','P6L2','P6R1','P6R2',...
    'P7L1','P7L2','P7R1','P7R2','P8L','P9L','P10L1','P10L2','P10R1','P10R2','P11L1','P11L2','P11R1',...
    'P11R2','P12R'};

myColor.RdBu = cbrewer('div', 'RdBu', 64);
myColor.RdBu = flipud(myColor.RdBu)
myColor.colors4 = [146 197 226;244 166 128;26 26 26;186 186 186]/255;
load('MLMTableTrial.mat')


%% Figure 2A
fun_RestMovePSD();

%% Figure 2B
fun_MovRelPow();

%% Figure 3A
fun_plot_Dot_Bar(MLMTableTrial,myColor,{'Ball'});

%% Figure 3B
fun_compareEMG(MLMTableTrial)

%% Figure 3C
fun_plotPosTorN()

%% Figure 3D
fun_corrMovBeta_TrainPerf(MLMTableTrial)

%% Figure 4A
fun_plot_Dot_Bar(MLMTableTrial,myColor,{'wBeta4sLFP_Nor'});
%% Figure 4B-D
fun_plot_Dot_Bar(MLMTableTrial,myColor,{'Burst4sDur1LFP','Burst4sDur2LFP','Burst4sNumLFP'});

%% Figure 4E
fun_plot_Dot_Bar(MLMTableTrial,myColor,{'wBeta4sEEG_Nor'});
%% Figure 4F-H
fun_plot_Dot_Bar(MLMTableTrial,myColor,{'Burst4sDur1EEG','Burst4sDur2EEG','Burst4sNumEEG'});

%% Figure 4I-J
fun_plot_Dot_Bar(MLMTableTrial,myColor,{'PIS','Coh'});

%% Figure 4-figure supplement 1
fun_plot_BurstProfile(MLMTableTrial);

%% Figure 5
fun_plot_Dot_Bar(MLMTableTrial,myColor,{'wGamma4sLFP_Nor','wPowBetaLFP_minus8_Nor','wPowBetaLFP_plus8_Nor'});

%% Figure 6A
fun_plot_Dot_Bar(MLMTableTrial,myColor,{'RT'});

%% Figure 6B
fun_plot_raw_RT(MLMTableTrial);

%% Figure 6C-D
fun_plot_Dot_Bar(MLMTableTrial(MLMTableTrial.idTremor==1,:),myColor,{'wPowTremorACC_Nor','wTremor4sLFP_Nor'});

%% Figure 6-figure supplement 1
fun_plot_MatchedBeta(MLMTableTrial,myColor,{'wBeta4sLFP_Nor','RT'});
fun_plot_MatchedBeta(MLMTableTrial,myColor,{'wBeta4sLFP_Nor','wGamma4sLFP_Nor'});
fun_plot_MatchedBeta(MLMTableTrial(MLMTableTrial.idTremor==1,:),myColor,{'wBeta4sLFP_Nor','wPowTremorACC_Nor'});

%% Figure 6-figure supplement 2
fun_Tremor_LFPtheta(MLMTableTrial(MLMTableTrial.idTremor==1,:))

%% Figure 7A,B,C,E
MLMTableSession = fun_mlmSession(MLMTableTrial);
fun_day1day2Learning(MLMTableSession);

%% 7D
fun_day1day2Baseline();

%% Table II
fun_glme(MLMTableTrial)

function fun_glme(MLMTableTrial)
myDistribution = 'InverseGaussian';
linkFun = 'identity';
residualsType = 'fitted';
corrType = 'Pearson';
idxRTNaN = find(isnan(MLMTableTrial.RT) | MLMTableTrial.RT<=0.2 | MLMTableTrial.RT>=2);
MLMTableTrial.RT(idxRTNaN) = nan;

%1. RT ~ 1 + k * TorN + 1 | SubID
lme1 = fitglme(MLMTableTrial,'RT~ idTorN + (1|idSub)',...
    'Distribution',myDistribution,'Link',linkFun)
R2 = lme1.Rsquared

%2. RT ~ 1 + k * betaLFP + 1 | SubID
lme2 = fitglme(MLMTableTrial,'RT~ wBetaP2sLFP_dB + (1|idSub)',...
    'Distribution',myDistribution,'Link',linkFun)
R2 = lme2.Rsquared

%3. RT ~ 1 + k * Dur1LFP + 1 | SubID
lme3 = fitglme(MLMTableTrial,'RT~ BurstP2sDur1LFP + (1|idSub)',...
    'Distribution',myDistribution,'Link',linkFun)
R2 = lme3.Rsquared

%4. RT ~ 1 + k * Dur2LFP + 1 | SubID
lme4 = fitglme(MLMTableTrial,'RT~ BurstP2sDur2LFP + (1|idSub)',...
    'Distribution',myDistribution,'Link',linkFun)
R2 = lme4.Rsquared

%5. RT ~ 1 + k * NumLFP + 1 | SubID
lme5 = fitglme(MLMTableTrial,'RT~ BurstP2sNumLFP + (1|idSub)',...
    'Distribution',myDistribution,'Link',linkFun)
R2 = lme5.Rsquared

%6. RT ~ 1 + k1 * TorN + k2 * betaLFP + k3 * gammaLFP + k4 * alphaLFP  + 1 | SubID
lme6 = fitglme(MLMTableTrial,'RT~ idTorN+wBetaP2sLFP_dB+wGammaP2sLFP_dB+wAlphaP2sLFP_dB+(1|idSub) ',...
    'Distribution',myDistribution,'Link',linkFun)
R2 = lme6.Rsquared

%7. RT ~ 1 + k * betaEEG + 1 | SubID
lme7 = fitglme(MLMTableTrial,'RT~ wBetaP2sEEG_dB + (1|idSub)',...
    'Distribution',myDistribution,'Link',linkFun)
R2 = lme7.Rsquared

%8. RT ~ 1 + k1 * TorN + k2 * betaEEG + k3 * alphaEEG + 1 | SubID
lme8 = fitglme(MLMTableTrial,'RT~ idTorN+wBetaP2sEEG_dB + wAlphaP2sEEG_dB + (1|idSub) ',...
    'Distribution',myDistribution,'Link',linkFun)
R2 = lme8.Rsquared

%9. RT ~ 1 + k1 * TorN + k2 * betaLFP + k3 * gammaLFP + k4 * betaEEG + 1 | SubID
lme9 = fitglme(MLMTableTrial,'RT~ idTorN+wBetaP2sLFP_dB+wGammaP2sLFP_dB+wBetaP2sEEG_dB+(1|idSub) ',...
    'Distribution',myDistribution,'Link',linkFun)
R2 = lme9.Rsquared

function fun_plot_Dot_Bar(MLMTableTrial,myColor,variableName)
idxTremor = [0 0 1 1 1 1 0 0 0 0 0 0 1 1 0 0 0 0 1 1 1];
setSub = unique(MLMTableTrial.idSub);
nSub = length(setSub);
myFontsize = 12;
figure
for iVar = 1:length(variableName)
    idxVar = find(strcmp(MLMTableTrial.Properties.VariableNames,variableName{iVar}));
    value_T = [];
    value_N = [];
    pValue = [];
    for iSub=1:nSub
%         idxT = find(MLMTableTrial.idSub==setSub(iSub) & MLMTableTrial.idTorN==1 & MLMTableTrial.idDay==1);
%         idxN = find(MLMTableTrial.idSub==setSub(iSub) & MLMTableTrial.idTorN==0 & MLMTableTrial.idDay==1);
        idxT = find(MLMTableTrial.idSub==setSub(iSub) & MLMTableTrial.idTorN==1);
        idxN = find(MLMTableTrial.idSub==setSub(iSub) & MLMTableTrial.idTorN==0);
        value_T(1,iSub) = nanmean(MLMTableTrial{:,:}(idxT,idxVar),1);
        value_T(2,iSub) = nanstd(MLMTableTrial{:,:}(idxT,idxVar),1)./sqrt(length(idxT));
        value_N(1,iSub) = nanmean(MLMTableTrial{:,:}(idxN,idxVar),1);
        value_N(2,iSub) = nanstd(MLMTableTrial{:,:}(idxN,idxVar),1)./sqrt(length(idxN));
        [~,pValue(iSub)] = ttest2(MLMTableTrial{:,:}(idxN,idxVar),MLMTableTrial{:,:}(idxT,idxVar));
    end
    rtDiff = squeeze(value_T(1,:)-value_N(1,:));
    [~,idsort] = sort(rtDiff);
    rtDiff = rtDiff(idsort);
    value_T = (squeeze(value_T(:,idsort)));
    value_N = (squeeze(value_N(:,idsort)));
    subplot(3,4,4*(iVar-1)+1)
    hold on
    axis(gca, 'square');
    main_fig        = findobj(gca,'Type','axes');
    axpos           = get(main_fig, 'Position');
    % make the main figure a bit smaller
    axpos(3)        = axpos(3) * 0.85;
    axpos(4)        = axpos(4) * 0.85;
    set(main_fig, 'Position', axpos);
    h_inset         = copyobj(main_fig, main_fig.Parent);
    idx1 = 0;
    idx2 = nSub;
    for iSub=1:nSub
        h = ploterr(squeeze(value_N(1,iSub)), squeeze(value_T(1,iSub)), ...
            squeeze(value_N(2,iSub)), ...
            squeeze(value_T(2,iSub)), ...
            '.', 'abshhxy', 0);
        % set colors according to yet another scale
        if(rtDiff(iSub)<0)
            idx1 = idx1 + 1;
            set(h(1), 'color', myColor.colors4(3, : ), 'markersize', 20);
            set(h(2), 'color', myColor.colors4(3, : ), 'linewidth', 0.5);
            set(h(3), 'color', myColor.colors4(3, : ), 'linewidth', 0.5);
        else
            idx2 = idx2 - 1;
            set(h(1), 'color', myColor.colors4(4, : ), 'markersize', 20);
            set(h(2), 'color', myColor.colors4(4, : ), 'linewidth', 0.5);
            set(h(3), 'color', myColor.colors4(4, : ), 'linewidth', 0.5);
        end
    end
    xlims   = get(main_fig, 'xlim');
    ylims   = get(main_fig, 'ylim');
    xlabel('No-training');
    ylabel('Training');
    minAx   = min([xlims ylims]);
    maxAx   = max([xlims ylims]);
    xlim(main_fig, [minAx maxAx]);
    ylim(main_fig, [minAx maxAx]);
    set(main_fig, 'xtick', get(main_fig, 'ytick'),'FontSize',myFontsize);
    l = refline(1,0);
    set(l, 'color', 'k', 'linestyle', '-', 'linewidth', 0.5);
    hold off
    xlim(main_fig, [minAx maxAx]);
    ylim(main_fig, [minAx maxAx]);
    set(main_fig, 'xtick', get(main_fig, 'ytick'));
    diff_N_T = squeeze(value_N(1,:)-value_T(1,:));
    hHist = bar(h_inset,[-0.25 0.25],[length(find(diff_N_T<0)) 0],0.5,'FaceColor',myColor.colors4(4,:),'EdgeColor','k');
    hold on
    hHist = bar(h_inset,[-0.25 0.25],[0 length(find(diff_N_T>0))],0.5,'FaceColor',myColor.colors4(3,:),'EdgeColor','k');
    hold off
    insetSize = axpos(3) * 0.7;
    rightTopX = axpos(1) + axpos(3) - 0.4*insetSize;
    rightTopY = axpos(2) + axpos(4) - 0.05*insetSize;
    set(h_inset,'view', [45 90], ...
        'Position', [rightTopX rightTopY insetSize insetSize], ...
        'box', 'off', 'ytick', [], 'ycolor', 'w', 'fontsize', 3, ...
        'xlim',[-0.75 0.75],'xtick',[],...
        'FontSize',11);
    
    subplot(3,4,4*(iVar-1)+2)
    Y = [nanmean(value_N(1,:)) nanmean(value_T(1,:))];
    SEM = [std(value_N(1,:))/sqrt(nSub) std(value_T(1,:))/sqrt(nSub)];
    hBar1 = bar([Y(1) 0],'EdgeColor','k','FaceColor',myColor.colors4(1, : ));
    hold on
    hBar2 = bar([0 Y(2)],'EdgeColor','k','FaceColor',myColor.colors4(2, : ));
    errorbar([1 2],Y,SEM,'o','color',[0.5 0.5 0.5],'linewidth',2)
    hBar1.BarWidth = 0.4;
    hBar2.BarWidth = 0.4;
    [h,pNormal] = adtest(value_N(1,:)-value_T(1,:));
    if(pNormal>0.05)
        [h,p,ci,stats] = ttest(value_T(1,:),value_N(1,:));
    else
       [p,h,stats] = signrank(value_T(1,:),value_N(1,:))
    end
    if(p<=0.05)
        H=sigstar({[1,2]},p);
        set(H,'Color','b')
    end
    xticks([1:2])
    xlim([0.5 2.5])
    xticklabels({'N','T'});
    ylabel(variableName{iVar})
    if(pNormal>0.05)
        title(['t' num2str(stats.df) '=' num2str(stats.tstat) ',p=' num2str(p) ', p2=' num2str(pNormal)]);
    else
        if(isfield(stats,'zval'))
            title(['z=' num2str(stats.zval) ',p=' num2str(p) ', p2=' num2str(pNormal)]);
        else
            title(['p=' num2str(p) ', p2=' num2str(pNormal)]);
        end
    end
    set(gca,'FontSize',myFontsize);
%     plot([1 2],[value_N(1,:)' value_T(1,:)'],'Color','k','LineStyle','-.');
%     plot(1,value_N(1,:)','Marker','^','MarkerSize',6,'MarkerEdgeColor','k');
%     plot(2,value_T(1,:)','Marker','^','MarkerSize',6,'MarkerEdgeColor','k','MarkerFaceColor',[187 187 187]/255);
    hold off
    axis square
    pValue
end

function fun_corrMovBeta_TrainPerf(MLMTableTrial)
corrType = 'Pearson'
nameYLFP = {'wBeta4sLFP_Nor'};
subSet = unique(MLMTableTrial.idSub);
nSub = length(subSet);
idxYLFP = [];
for iVar =1:length(nameYLFP)
    idxYLFP(iVar) = find(strcmp(MLMTableTrial.Properties.VariableNames,nameYLFP{iVar}));
end
valueYLFP = []
for iSub = 1:nSub
    valueYLFP(iSub,:) = nanmean(MLMTableTrial{:,:}(MLMTableTrial.idSub==subSet(iSub) & MLMTableTrial.idTorN==1,idxYLFP))-...
        nanmean(MLMTableTrial{:,:}(MLMTableTrial.idSub==subSet(iSub) & MLMTableTrial.idTorN==0,idxYLFP));
end
load('betaReduced.mat')
BetaRedLFP = []
for iSub = 1:length(betaReduced)
    BetaRedLFP = [BetaRedLFP;betaReduced{iSub}*100];
end
% Correlation LFP
nRow = size(BetaRedLFP,2);
nCol = length(idxYLFP);
figure
for iRow = 1:nRow
    for iPlot = 1:nCol
        subplot(nRow,nCol,(iRow-1)*nCol+iPlot)
        pointX = BetaRedLFP(:,iRow);
        pointY = -valueYLFP(:,iPlot);
        scatter(pointX,pointY,40,'MarkerFaceColor','r','MarkerEdgeColor','r',...
            'MarkerFaceAlpha',.25,'MarkerEdgeAlpha',.25);
        hold on
        [r,p] = corr(pointX,pointY,'rows','complete','Type',corrType);
        fitobject = fit(pointX,pointY,'poly1');
        h = plot(fitobject,pointX,pointY)
        hold off
        set(h(2),'LineWidth',3,'Color','k');
        delete(h(1))
        xlabel('MovBetaReduction');
        ylabel(['Diff in ' nameYLFP{iPlot}])
        title(['r = ' num2str(r,4) ', p = ' num2str(p)]);
        axis square
        set(gca,'box','on');
    end
end

function fun_day1day2Learning(MLMTableSession)

MLMTableSessionDiff = MLMTableSession(MLMTableSession.idTorN==1,:);
MLMTableSessionDiff{:,:}(:,8:end) = MLMTableSession{:,:}(MLMTableSession.idTorN==1,8:end) - MLMTableSession{:,:}(MLMTableSession.idTorN==0,8:end);
% MLMTableSessionDiff{:,:}(:,8:end) = MLMTableSession{:,:}(MLMTableSession.idTorN==1,8:end);

varName = {'Ball','wBeta4sLFP_Nor','Burst4sDur1LFP'};
idxVar = [];
for iVar = 1:length(varName)
    idxVar = [idxVar find(strcmp(MLMTableSessionDiff.Properties.VariableNames,varName{iVar}))];
end
subSet = unique(MLMTableSessionDiff.idSub(MLMTableSessionDiff.idTwoDays==1));
valDay1 = []; valDay2 = [];
for iSub = 1:length(subSet)
    idxRange = find(MLMTableSessionDiff.idSub==subSet(iSub) & MLMTableSessionDiff.idDay==1);
    valDay1(iSub,:) = mean(MLMTableSessionDiff{:,:}(idxRange,idxVar))
    idxRange = find(MLMTableSessionDiff.idSub==subSet(iSub) & MLMTableSessionDiff.idDay==2);
    valDay2(iSub,:) = mean(MLMTableSessionDiff{:,:}(idxRange,idxVar))
end
figure
% myColor = cbrewer('div', 'RdBu', 20);
% myColor = flipud(myColor)
% myColor = myColor([1:4 17:20],:);
myColor(1,:) = [180 241 116];
myColor(2,:) = [246 184 11];
myColor(3,:) = [0 0 0];
myColor(4,:) = [35 222 148];
myColor(5,:) = [255 0 0];
myColor(6,:) = [159 89 131];
myColor(7,:) = [246 184 253];
myColor(8,:) = [0 129 255];
myColor = myColor/255;
for iPlot = 1:length(varName)
    subplot(3,6,iPlot)
    hold on
    for iLine = 1:size(valDay1,1)
        line([1 2],[valDay1(iLine,iPlot) valDay2(iLine,iPlot)],'Color',myColor(iLine,:),'Marker','o','LineWidth',3);
    end
    xlim([0 3])
    title(varName{iPlot})
    xticks([1 2])
    xticklabels({'Day1','Day2'});
    if(iPlot==1);legend({'H11','H12','H13','H14','H17','H18','H19','H20'});end
    set(gca,'box','on','FontSize',15,'LineWidth',2);
    
    
    subplot(3,6,6+iPlot)
    hBar1 = bar([mean(valDay1(:,iPlot)) 0],'EdgeColor','k','FaceColor',[1 1 1],'LineWidth',3)
    hold on
    hBar2 = bar([0 mean(valDay2(:,iPlot))],'EdgeColor','k','FaceColor',[1 1 1],'LineWidth',3)
    errorbar([1 2],[mean(valDay1(:,iPlot)) mean(valDay2(:,iPlot))],...
        [std(valDay1(:,iPlot))/sqrt(length(valDay1(:,iPlot))) std(valDay2(:,iPlot))/sqrt(length(valDay2(:,iPlot)))],...
        'o','color',[0.5 0.5 0.5],'linewidth',2)
    hBar1.BarWidth = 0.4;
    hBar2.BarWidth = 0.4;
    [p,h] = signrank(valDay1(:,iPlot),valDay2(:,iPlot));
    if(p<=0.05)
        H=sigstar({[1,2]},p);
        set(H,'Color','b')
    end
    xlim([0 3])
    if(iPlot==1);ylabel('Diff Between Training and No Training');end
    xticklabels({'Day1','Day2'});
    set(gca,'FontSize',15,'LineWidth',2);
end

varName = {'wPowTremorACC_Nor'};
idxVar = [];
for iVar = 1:length(varName)
    idxVar = [idxVar find(strcmp(MLMTableSessionDiff.Properties.VariableNames,varName{iVar}))];
end
subSet = unique(MLMTableSessionDiff.idSub(MLMTableSessionDiff.idTwoDays==1 & MLMTableSessionDiff.idTremor==1));
valDay1 = []; valDay2 = [];
for iSub = 1:length(subSet)
    idxRange = find(MLMTableSessionDiff.idSub==subSet(iSub) & MLMTableSessionDiff.idDay==1);
    valDay1(iSub,:) = mean(MLMTableSessionDiff{:,:}(idxRange,idxVar))
    idxRange = find(MLMTableSessionDiff.idSub==subSet(iSub) & MLMTableSessionDiff.idDay==2);
    valDay2(iSub,:) = mean(MLMTableSessionDiff{:,:}(idxRange,idxVar))
end


subplot(3,6,iPlot+1)
myColorTremor = myColor([3 4 7 8],:);
hold on
for iLine = 1:size(valDay1,1)
    line([1 2],[valDay1(iLine,1) valDay2(iLine,1)],'Color',myColorTremor(iLine,:),'Marker','o','LineWidth',3);
end
xlim([0 3])
title(varName{1})
xticks([1 2])
xticklabels({'Day1','Day2'});
set(gca,'box','on','FontSize',15,'LineWidth',2);


subplot(3,6,6+iPlot+1)
hBar1 = bar([mean(valDay1(:,1)) 0],'EdgeColor','k','FaceColor',[1 1 1],'LineWidth',3)
hold on
hBar2 = bar([0 mean(valDay2(:,1))],'EdgeColor','k','FaceColor',[1 1 1],'LineWidth',3)
errorbar([1 2],[mean(valDay1(:,1)) mean(valDay2(:,1))],...
    [std(valDay1(:,1))/sqrt(length(valDay1(:,1))) std(valDay2(:,1))/sqrt(length(valDay2(:,1)))],...
    'o','color',[0.5 0.5 0.5],'linewidth',2)
hBar1.BarWidth = 0.4;
hBar2.BarWidth = 0.4;
[p,h] = signrank(valDay1(:,1),valDay2(:,1));
if(p<=0.05)
    H=sigstar({[1,2]},p);
    set(H,'Color','b')
end
xlim([0 3])
xticklabels({'Day1','Day2'});
set(gca,'FontSize',15,'LineWidth',2);

function fun_MovRelPow()
load('MovRelPow.mat')
x = linspace(-2+1/2048,3,size(avgPower2,2));
y = 1:1:100;
figure
subplot(2,2,4)
beta = [15 15 20 20 20 20 21 20 15 15 15 25 20 24 15 19 22 21 19 19 18];
nSmooth = 1;
for i=1:size(avgPower2,1)
    temp = smooth(avgPower2(i,:),nSmooth)';
    avgPower3(i,:) = temp;
end
bl_power = nanmean(avgPower3(:,idxBL),2);
avgPower3 = 100*(bsxfun(@rdivide,avgPower3,bl_power)-1);
imagesc(x,y,avgPower3);
hold on
rectangle('Position',[0,mean(beta)-2,3,5],...
    'Curvature',[0.8,0.4],...
    'LineWidth',2,'LineStyle','--')
line([0 0],[1 95],'LineStyle','--','Color','r');
hold off
colorbar
caxis([-20 30]);
set(gca,'YDir','normal') 
title('Training','FontSize',20);
xlabel('Time (s)','FontSize',20);
ylabel('Frequency (Hz)','FontSize',20);
xlim([-1 2])
ylim([1 95])

function fun_compareEMG(MLMTableTrial)
load('raw_EMG_avg.mat')
exEMG = find((MLMTableTrial.idSub==8 & MLMTableTrial.idTrial==6) | (MLMTableTrial.idSub==12 & MLMTableTrial.idDay==1 & MLMTableTrial.idTrial==50));
fs = 2048;
% nSub = length(EMG_trial)
trial_T = [];
trial_N = [];
avg_T = [];
avg_N = [];
% valid_Trial = ~isnan(sum(MLMTableTrial{:,[1:21 23:end]},2));
valid_Trial = ~isnan(sum(MLMTableTrial{:,:},2));
valid_Trial(exEMG) = 0;
timeAll = linspace(-2,4,fs*6);

% for iSub = 1:nSub
%     % Training
%     idxTrial = find(MLMTableTrial.idSub==iSub & MLMTableTrial.idDay==1 & MLMTableTrial.idTorN==1);
%     if(sum(MLMTableTrial.idTremor(idxTrial))>0) % Exclude subjects with consistent tremor during recording.
%         continue; 
%     end
%     figure
%     subplot(2,2,1)
%     idxTrial = idxTrial(valid_Trial(idxTrial)~=0)-find(MLMTableTrial.idSub==iSub & MLMTableTrial.idTrial==1 & MLMTableTrial.idDay==1)+1; % Exlude in-valid trials
% %     emg_T = EMG_trial{1,iSub}(idxTrial,:)-mean(EMG_trial{1,iSub}(idxTrial,1:4096),2);
%     emg_T = EMG_trial{1,iSub}(idxTrial,:);
%     trial_T = [trial_T;emg_T];
%     avg_T = [avg_T; mean(emg_T)];
%     plot(timeAll,emg_T','-.');
%     hold on
%     plot(timeAll,mean(emg_T),'-','LineWidth',2,'Color','k');
%     hold off
%     title(['Sub' num2str(iSub) ': Training'])
%     xlabel('Time (s)')
%     
%     % Training
%     subplot(2,2,2)
%     idxTrial = find(MLMTableTrial.idSub==iSub & MLMTableTrial.idDay==1 & MLMTableTrial.idTorN==0);
%     idxTrial = idxTrial(valid_Trial(idxTrial)~=0)-find(MLMTableTrial.idSub==iSub & MLMTableTrial.idTrial==1 & MLMTableTrial.idDay==1)+1; % Exlude in-valid trials
% %     emg_N = EMG_trial{1,iSub}(idxTrial,:)-mean(EMG_trial{1,iSub}(idxTrial,1:4096),2);
%     emg_N = EMG_trial{1,iSub}(idxTrial,:);
%     trial_N = [trial_N;emg_N];
%     avg_N = [avg_N; mean(emg_N)];
%     plot(timeAll,emg_N','-.');
%     hold on
%     plot(timeAll,mean(emg_N),'-','LineWidth',2,'Color','k');
%     hold off
%     title(['Sub' num2str(iSub) ': No training'])
%     xlabel('Time (s)')
%     
%     subplot(2,1,2)
%     bar(mean(EMG_trial{1,iSub},2))
%     xlabel('Trials')
% end
% figure
% sup1 = subplot(4,3,1)
% plot(timeAll,trial_T','-.');
% hold on
% plot(timeAll,mean(trial_T),'-','LineWidth',2,'Color','k');
% hold off
% title('Training')
% xlabel('Time (s)')
% set(gca,'box','off')
% 
% sup2 = subplot(4,3,4)
% plot(timeAll,trial_N','-.');
% hold on
% plot(timeAll,mean(trial_N),'-','LineWidth',2,'Color','k');
% hold off
% title('No training')
% xlabel('Time (s)')
% set(gca,'box','off')
% linkaxes([sup1 sup2],'y')
% 
% subplot(3,3,7)
% idxNFB = find(timeAll>=0 & timeAll<=4);
% avgEMG_R = mean(trial_N(:,idxNFB),2);
% ems_R = std(avgEMG_R)/sqrt(length(avgEMG_R));
% avgEMG_T = mean(trial_T(:,idxNFB),2);
% ems_T = std(avgEMG_T)/sqrt(length(avgEMG_T));
% hBar1 = bar([mean(avgEMG_R) 0],'EdgeColor','k','FaceColor',[146 197 226]./255)
% hold on
% hBar2 = bar([0 mean(avgEMG_T)],'EdgeColor','k','FaceColor',[244 166 128]./255)
% errorbar([1 2],[mean(avgEMG_R) mean(avgEMG_T)],[ems_R ems_T],...
%     'o','color',[0.5 0.5 0.5],'linewidth',2)
% hBar1.BarWidth = 0.4;
% hBar2.BarWidth = 0.4;
% [h,p,ci,stats] = ttest2(avgEMG_R,avgEMG_T);
% if(p<=0.05)
%     H=sigstar({[1,2]},p);
%     set(H,'Color','b')
% end
% xticks([1:2])
% xlim([0.5 2.5])
% % ylabel('Power [dB]')
% xticklabels({'N','T'});
% set(gca,'FontSize',20);
% title(['t' num2str(stats.df) '=' num2str(stats.tstat) ',p=' num2str(p)]);
% axis square

% save('raw_EMG_avg.mat','avg_T','avg_N');
load('raw_EMG_avg.mat');
for i=1:12
    myColor(i,:) = [rand() rand() rand()];
end
sup3 = subplot(4,3,3)
hold on
for iAvg = 1:size(avg_T,1)
    hLine = plot(timeAll,avg_T(iAvg,:),'-','Color',myColor(iAvg,:));
    hLine.Color(4) = 0.5;
end
plot(timeAll,mean(avg_T),'-','Color','k','LineWidth',2)
hold off

sup4 = subplot(4,3,6)
hold on
for iAvg = 1:size(avg_N,1)
    hLine = plot(timeAll,avg_N(iAvg,:),'-','Color',myColor(iAvg,:));
    hLine.Color(4) = 0.5;
end
plot(timeAll,mean(avg_N),'-','Color','k','LineWidth',2)
hold off
linkaxes([sup3 sup4],'y')

subplot(3,3,9)
idxNFB = find(timeAll>=0 & timeAll<=4);
avgEMG_R = mean(avg_N(:,idxNFB),2);
ems_R = std(avgEMG_R)/sqrt(length(avgEMG_R));
avgEMG_T = mean(avg_T(:,idxNFB),2);
ems_T = std(avgEMG_T)/sqrt(length(avgEMG_T));
hBar1 = bar([mean(avgEMG_R) 0],'EdgeColor','k','FaceColor',[146 197 226]./255)
hold on
hBar2 = bar([0 mean(avgEMG_T)],'EdgeColor','k','FaceColor',[244 166 128]./255)
errorbar([1 2],[mean(avgEMG_R) mean(avgEMG_T)],[ems_R ems_T],...
    'o','color',[0.5 0.5 0.5],'linewidth',2)
hBar1.BarWidth = 0.4;
hBar2.BarWidth = 0.4;
[h,p,ci,stats] = ttest(avgEMG_R,avgEMG_T);
if(p<=0.05)
    H=sigstar({[1,2]},p);
    set(H,'Color','b')
end
xticks([1:2])
xlim([0.5 2.5])
% ylabel('Power [dB]')
xticklabels({'N','T'});
set(gca,'FontSize',20);
title(['t' num2str(stats.df) '=' num2str(stats.tstat) ',p=' num2str(p)]);
axis square

function fun_plotPosTorN()
load('PowTorN.mat')

figure
x = freqs2use;
y1 = smooth(nanmean(psd_N2,1),3);
e1 = std(squeeze(nanmean(psd_N2,1)))/sqrt(21);
y2 = smooth(nanmean(psd_T2,1),3);
e2 = std(squeeze(nanmean(psd_T2,1)))/sqrt(21);

for i=1:95
    [p(i)] = signrank(psd_N2(:,i),psd_T2(:,i));
end
subplot(2,4,[5:8]);
[l1,p1] = boundedline(x, y1, e1,'alpha');
set(l1,'Color',[0 114 189]/255,'LineWidth',4);
set(p1,'FaceColor',[0 114 189]/255);
hold on
[l2,p2] = boundedline(x, y2, e2,'alpha');
set(l2,'Color',[217 83 25]/255,'LineWidth',4);
set(p2,'FaceColor',[217 83 25]/255);

hold off
axis tight;
xlim([3 95])
xlabel('Frequency (Hz)','FontSize',20);
ylabel('Power (% against baseline)','FontSize',20);
legend([l1 l2],{'Notraining','Training'},'FontSize',20)
set(gca,'box','on','LineWidth',2,'FontSize',20);

function fun_RestMovePSD()
load('RestMovPSD.mat')
beta = [15 15 20 20 20 20 21 20 15 15 15 25 20 24 15 19 22 21 19 19 18];
figure
subplot(2,2,3)
x = F;
y1 = 10*log10(smooth(nanmean(restPowLFP,1),4)');
e1 = smooth(std(10*log10(restPowLFP),1)/sqrt(21),4)';
y2 = 10*log10(smooth(nanmean(movePowLFP,1),4)');
e2 = smooth(std(10*log10(movePowLFP),1)/sqrt(21),4)';
[l1,p1] = boundedline(x, y1, e1,'alpha');
set(l1,'Color',[0 0 0]/255,'LineWidth',4);
set(p1,'FaceColor',[0 0 0]/255);
hold on
[l2,p2] = boundedline(x, y2, e2,'alpha');
x = [mean(beta)-2 mean(beta)+2 mean(beta)+2 mean(beta)-2];
y = [-20 -20 0 0];
patch(x,y,[0.5 0.5 0.5],'facealpha',0.1)
hold off
set(l2,'Color',[255 0 0]/255,'LineWidth',4);
set(p2,'FaceColor',[255 0 0]/255);
xlim([5 40])
ylim([-18 6])
xlabel('Frequency (Hz)','FontSize',20);
ylabel('dB/Hz','FontSize',20);
title('LFP beta power spectral','FontSize',20)
legend([l1 l2],{'Resting','Movement-related'},'FontSize',16);
set(gca, 'box', 'on','LineWidth',2,'FontSize',20);

function tableOutput = fun_mlmSession(tableInput)
nSub = length(unique(tableInput.idSub));
tableOutput = [];
for iSub = 1:nSub
    idxRange1 = find(tableInput.idSub==iSub);
    tableTemp1 = tableInput(idxRange1,:);
    nDay = length(unique(tableTemp1.idDay));
    tableDay = [];
    for iDay = 1:nDay
        idxRange2 = find(tableTemp1.idDay==iDay);
        tableTemp2 = tableTemp1(idxRange2,:);
        nSession = length(unique(tableTemp2.idSession));
        tableSession = [];
        for iSession = 1:nSession
            idxRangeT = find(tableTemp2.idSession==iSession & tableTemp2.idTorN==1);
            idxRangeN = find(tableTemp2.idSession==iSession & tableTemp2.idTorN==0);
            tableItem = [];
            for iItem = 1:length(tableTemp2.Properties.VariableNames)
                tableItem = [tableItem table([nanmean(tableTemp2{idxRangeT,iItem},1);nanmean(tableTemp2{idxRangeN,iItem},1)],...
                    'VariableNames',tableTemp2.Properties.VariableNames(iItem))];
            end
            tableSession = [tableSession;tableItem];
        end
        tableDay = [tableDay;tableSession];
    end
    tableOutput = [tableOutput;tableDay];
    fprintf('Sub:%d\n',iSub);
end

function fun_plot_raw_RT(MLMTableTrial)
load('raw_FRC.mat')
exFRC = find((MLMTableTrial.RT<0.2 | MLMTableTrial.RT>2));
Hemisphere = [1:21];
nFile = length(raw_FRC)
valid_Trial = ~isnan(sum(MLMTableTrial{:,:},2));
valid_Trial(exFRC) = 0;
index = find(timeTrial>=-0.5 & timeTrial<=2);
X = timeTrial(index);
idCue = find(X<=0);
idCue = idCue(end);
trial_T = [];
trial_N = [];
for iFile = 21
    figure
    sub1 = subplot(2,2,1)
    idxHem = find(MLMTableTrial.idSub==Hemisphere(iFile));
    idxTrial = find(MLMTableTrial.idSub==Hemisphere(iFile)& MLMTableTrial.idTorN==1);
    idxTrial = idxTrial(valid_Trial(idxTrial)~=0)-idxHem(1)+1; % Exlude in-valid trials
    dataTemp = raw_FRC{1,iFile}(idxTrial,index);
    frc_T = dataTemp - mean(dataTemp(:,1:idCue),2);
    trial_T = [trial_T;frc_T];
    plot(X,frc_T','-.');
    hold on
    plot(X,mean(frc_T),'-','LineWidth',2,'Color','k');
    hold off
    title(['Hem' num2str(iFile) ': Training'])
    xlabel('Time (s)')
    axis square
    
    sub2 = subplot(2,2,2)
    idxTrial = find(MLMTableTrial.idSub==Hemisphere(iFile) & MLMTableTrial.idTorN==0);
    idxTrial = idxTrial(valid_Trial(idxTrial)~=0)-idxHem(1)+1; % Exlude in-valid trials
    dataTemp = raw_FRC{1,iFile}(idxTrial,index);
    frc_N = dataTemp - mean(dataTemp(:,1:idCue),2);
    trial_N = [trial_N;frc_N];
    plot(X,frc_N','-.');
    hold on
    plot(X,mean(frc_N),'-','LineWidth',2,'Color','k');
    hold off
    title(['Hem' num2str(iFile) ': No Training'])
    xlabel('Time (s)')
    axis square
    
    
    sub3 = subplot(2,2,3)
    plot(X,frc_N','-.','Color',[66 149 193 200]/255);
    hold on
    plot(X,frc_T','-.','Color',[244 166 128 200]/255);
    plot(X,mean(frc_N),'-','LineWidth',3,'Color',[0 0 0 200]/255);
    plot(X,mean(frc_T),'-','LineWidth',3,'Color',[255 0 0 125]/255);
    hold off
    title(['Hem' num2str(iFile) ': Training and No Training'])
    xlabel('Time (s)')
    axis square
    set(gca,'FontSize',12,'LineWidth',1)
    linkaxes([sub1 sub2 sub3],'y')
end

function fun_plot_MatchedBeta(MLMTableTrial,myColor,variableName)
ratio1 =  .75; % a ratio to match beta power
ratio2 =  .25; % a ratio to match beta power

idxTremor = [0 0 1 1 1 1 0 0 0 0 0 0 1 1 0 0 0 0 1 1 1];
setSub = unique(MLMTableTrial.idSub);
nSub = length(setSub);
myFontsize = 12;
figure

idxBeta = find(strcmp(MLMTableTrial.Properties.VariableNames,variableName{1}));
beta_T = [];
beta_N = [];
idxRT = find(strcmp(MLMTableTrial.Properties.VariableNames,variableName{2}));
RT_T = [];
RT_N = [];
for iSub=1:nSub
    % select trials with higher beta in training condition
    idxT = find(MLMTableTrial.idSub==setSub(iSub) & MLMTableTrial.idTorN==1);
    [sort_beta, sort_trial] = sort(MLMTableTrial{:,:}(idxT,idxBeta),'descend');
    sort_beta(isnan(sort_beta)) = [];
    sort_trial(isnan(sort_beta)) = [];
    nValidTrial = floor(ratio1*length(sort_beta));
    beta_T(1,iSub) = nanmean(sort_beta(1:nValidTrial),1);
    beta_T(2,iSub) = nanstd(sort_beta(1:nValidTrial),1)./sqrt(nValidTrial);
    RT_T(1,iSub) = nanmean(MLMTableTrial{:,:}(sort_trial(1:nValidTrial),idxRT),1);
    RT_T(2,iSub) = nanstd(MLMTableTrial{:,:}(sort_trial(1:nValidTrial),idxRT))./sqrt(nValidTrial);
    
    % select trials with lower beta in no training condition
    idxN = find(MLMTableTrial.idSub==setSub(iSub) & MLMTableTrial.idTorN==0);
    [sort_beta, sort_trial] = sort(MLMTableTrial{:,:}(idxN,idxBeta),'ascend');
    sort_beta(isnan(sort_beta)) = [];
    sort_trial(isnan(sort_beta)) = [];
    nValidTrial = floor(ratio1*length(sort_beta));
    beta_N(1,iSub) = nanmean(sort_beta,1);
    beta_N(2,iSub) = nanstd(sort_beta,1)./sqrt(nValidTrial);
    RT_N(1,iSub) = nanmean(MLMTableTrial{:,:}(sort_trial(1:nValidTrial),idxRT),1);
    RT_N(2,iSub) = nanstd(MLMTableTrial{:,:}(sort_trial(1:nValidTrial),idxRT))./sqrt(nValidTrial);
end
% plot beta
fun_plot_individual_dot(beta_T,beta_N,[3 4 1],myColor);
fun_plot_line_bar(beta_T,beta_N,[3 4 2],myColor);
ylabel(variableName{1})

% plot RT
fun_plot_individual_dot(RT_T,RT_N,[3 4 5],myColor);
fun_plot_line_bar(RT_T,RT_N,[3 4 6],myColor);
ylabel(variableName{2})

function fun_plot_individual_dot(beta_T,beta_N,paramPlot,myColor)
nSub = size(beta_T,2);
subplot(paramPlot(1),paramPlot(2),paramPlot(3))
rtDiff = squeeze(beta_T(1,:)-beta_N(1,:));
[~,idsort] = sort(rtDiff);
rtDiff = rtDiff(idsort);
value_T = (squeeze(beta_T(:,idsort)));
value_N = (squeeze(beta_N(:,idsort)));
hold on
axis(gca, 'square');
main_fig        = findobj(gca,'Type','axes');
axpos           = get(main_fig, 'Position');
% make the main figure a bit smaller
axpos(3)        = axpos(3) * 0.85;
axpos(4)        = axpos(4) * 0.85;
set(main_fig, 'Position', axpos);
h_inset         = copyobj(main_fig, main_fig.Parent);
idx1 = 0;
idx2 = nSub;
for iSub=1:nSub
    h = ploterr(squeeze(value_N(1,iSub)), squeeze(value_T(1,iSub)), ...
        squeeze(value_N(2,iSub)), ...
        squeeze(value_T(2,iSub)), ...
        '.', 'abshhxy', 0);
    % set colors according to yet another scale
    if(rtDiff(iSub)<0)
        idx1 = idx1 + 1;
        set(h(1), 'color', myColor.colors4(3,:), 'markersize', 20);
        set(h(2), 'color', myColor.colors4(3,:), 'linewidth', 0.5);
        set(h(3), 'color', myColor.colors4(3,:), 'linewidth', 0.5);
    else
        idx2 = idx2 - 1;
        set(h(1), 'color', myColor.colors4(4,:), 'markersize', 20);
        set(h(2), 'color', myColor.colors4(4,:), 'linewidth', 0.5);
        set(h(3), 'color', myColor.colors4(4,:), 'linewidth', 0.5);
    end
end
xlims   = get(main_fig, 'xlim');
ylims   = get(main_fig, 'ylim');
xlabel('No-training');
ylabel('Training');
minAx   = min([xlims ylims]);
maxAx   = max([xlims ylims]);
xlim(main_fig, [minAx maxAx]);
ylim(main_fig, [minAx maxAx]);
set(main_fig, 'xtick', get(main_fig, 'ytick'));
l = refline(1,0);
set(l, 'color', 'k', 'linestyle', '-', 'linewidth', 0.5);
hold off
xlim(main_fig, [minAx maxAx]);
ylim(main_fig, [minAx maxAx]);
set(main_fig, 'xtick', get(main_fig, 'ytick'));
% hHist = histogram(h_inset, squeeze(value_N(1,:)-value_T(1,:)),'BinEdges',[-10 0 10],'BinWidth',0.5,'edgecolor', [0 0 0], ...
%     'facecolor', [1 1 1],'LineWidth',1);
diff_N_T = squeeze(value_N(1,:)-value_T(1,:));
hHist = bar(h_inset,[-0.25 0.25],[length(find(diff_N_T<0)) 0],0.5,'FaceColor',myColor.colors4(4,:),'EdgeColor','k')
hold on
hHist = bar(h_inset,[-0.25 0.25],[0 length(find(diff_N_T>0))],0.5,'FaceColor',myColor.colors4(3,:),'EdgeColor','k')
hold off
insetSize = axpos(3) * 0.7;
rightTopX = axpos(1) + axpos(3) - 0.4*insetSize;
rightTopY = axpos(2) + axpos(4) - 0.05*insetSize;
set(h_inset,'view', [45 90], ...
    'Position', [rightTopX rightTopY insetSize insetSize], ...
    'box', 'off', 'ytick', [], 'ycolor', 'w', 'fontsize', 3, ...
    'xlim',[-0.75 0.75],'xtick',[],...
    'FontSize',11);

function fun_plot_line_bar(beta_T,beta_N,paramPlot,myColor)
nSub = size(beta_T,2);
myFontsize = 12;
subplot(paramPlot(1),paramPlot(2),paramPlot(3))
rtDiff = squeeze(beta_T(1,:)-beta_N(1,:));
[~,idsort] = sort(rtDiff);
rtDiff = rtDiff(idsort);
value_T = (squeeze(beta_T(:,idsort)));
value_N = (squeeze(beta_N(:,idsort)));
Y = [nanmean(value_N(1,:)) nanmean(value_T(1,:))];
SEM = [std(value_N(1,:))/sqrt(nSub) std(value_T(1,:))/sqrt(nSub)]
hBar1 = bar([Y(1) 0],'EdgeColor','k','FaceColor',myColor.colors4(1,:))
hold on
hBar2 = bar([0 Y(2)],'EdgeColor','k','FaceColor',myColor.colors4(2,:))
errorbar([1 2],Y,SEM,'o','color',[0.5 0.5 0.5],'linewidth',2)
hBar1.BarWidth = 0.4;
hBar2.BarWidth = 0.4;
[h,pNormal] = adtest(value_N(1,:)-value_T(1,:))
if(pNormal>0.05)
    [h,p,ci,stats] = ttest(value_T(1,:),value_N(1,:));
else
    [p,h,stats] = signrank(value_T(1,:),value_N(1,:))
end
if(p<=0.05)
    H=sigstar({[1,2]},p);
    set(H,'Color','b')
end
xticks([1:2])
xlim([0.5 2.5])
xticklabels({'N','T'});
if(pNormal>0.05)
    title(['t20=' num2str(stats.tstat) ',p=' num2str(p) ', p2=' num2str(pNormal)]);
else
    if(isfield(stats,'zval'))
        title(['z=' num2str(stats.zval) ',p=' num2str(p) ', p2=' num2str(pNormal)]);
    else
        title(['p=' num2str(p) ', p2=' num2str(pNormal)]);
    end
end
set(gca,'FontSize',myFontsize);
% plot([1 2],[beta_N(1,:)' beta_T(1,:)'],'Color','k','LineStyle','-.');
% plot(1,beta_N(1,:)','Marker','^','MarkerSize',6,'MarkerEdgeColor','k');
% plot(2,beta_T(1,:)','Marker','^','MarkerSize',6,'MarkerEdgeColor','k','MarkerFaceColor',[187 187 187]/255);
hold off
axis square

function fun_plot_BurstProfile(MLMTableTrial)
load('burstIndividual.mat');
burstProfileT = [];
burstProfileN = [];
subSet = unique(MLMTableTrial.idSub)
for iFile = 1:length(subSet)
    idxTrial = find(MLMTableTrial.idSub==subSet(iFile));
    idxT = find(MLMTableTrial.idTorN(idxTrial)==1);
    valueTemp = reshape(burstIndividual2{1,iFile}(idxT,:),[],1);
    valueTemp(valueTemp==0) = [];
    burstT = valueTemp;
    
    idxN = find(MLMTableTrial.idTorN(idxTrial)==0);
    valueTemp = reshape(burstIndividual2{1,iFile}(idxN,:),[],1);
    valueTemp(valueTemp==0) = [];
    burstN = valueTemp;
    profileSub = zeros(9,2);
    for iSeg = 1:8
        profileSub(iSeg,1) = length(find(burstN>=iSeg*100 & burstN<(iSeg+1)*100));
        profileSub(iSeg,2) = length(find(burstT>=iSeg*100 & burstT<(iSeg+1)*100));
    end
    profileSub(iSeg+1,1) = length(find(burstN>=(iSeg+1)*100));
    profileSub(iSeg+1,2) = length(find(burstT>=(iSeg+1)*100));
    burstProfileN(iFile,:) = profileSub(:,1);
    burstProfileT(iFile,:) = profileSub(:,2);
end

figure
burstProfileN4 = [sum(burstProfileN(:,1:3),2) sum(burstProfileN(:,4:7),2) sum(burstProfileN(:,8:end),2)];
burstProfileT4 = [sum(burstProfileT(:,1:3),2) sum(burstProfileT(:,4:7),2) sum(burstProfileT(:,8:end),2)];
subplot(2,4,4)
hold on
for iSeg = 1:size(burstProfileN4,2)
    bar(iSeg-0.2, mean(burstProfileN4(:,iSeg)),'BarWidth',0.25,'FaceColor',[146 197 226]/255);
    bar(iSeg+0.2, mean(burstProfileT4(:,iSeg)),'BarWidth',0.25,'FaceColor',[244 166 128]/255);
    SEM1 = std(burstProfileN4(:,iSeg))/sqrt(length(subSet));
    SEM2 = std(burstProfileT4(:,iSeg))/sqrt(length(subSet));
    errorbar([iSeg-0.2 iSeg+0.2],[mean(burstProfileN4(:,iSeg)) mean(burstProfileT4(:,iSeg))],[SEM1 SEM2],...
        'o','color',[0.5 0.5 0.5],'linewidth',2)
    [h,p,ci,stats] = ttest(burstProfileN4(:,iSeg),burstProfileT4(:,iSeg));
    p = p*3;
    if(p<0.05)
        H=sigstar({[iSeg-0.2,iSeg+0.2]},p);
        set(H,'Color','b')
    end
end
hold off
xlim([0.3 3.7])
xticks([1:3])
xticklabels({'<0.4','0.4-0.8','>0.8'});
xlabel('Time windows (s)')
ylabel('Amount of bursts for all trials')
legend({'No training','Training'})

function fun_Tremor_LFPtheta(MLMTableTrial)
subSet = unique(MLMTableTrial.idSub)
myDistribution = 'Normal';
linkFun = 'identity';
fitMethod = 'MPL';

myColor = distinguishable_colors(length(subSet)+1);
figure
subplot(3,4,2)
hold on
fitVal = [];reaVal = [];
for iSub=1:length(subSet)
    idxSubT = find(MLMTableTrial.idSub==subSet(iSub) & MLMTableTrial.idTorN==1);
    idxSubN = find(MLMTableTrial.idSub==subSet(iSub) & MLMTableTrial.idTorN==0);
    fitVal(end+1,:) = [10*log10(nanmean(10.^(MLMTableTrial.wPowTremorACC_dB(idxSubT)/10))) 10*log10(nanmean(10.^(MLMTableTrial.wPowTremorACC_dB(idxSubN)/10)))];
    reaVal(end+1,:) = [10*log10(nanmean(10.^(MLMTableTrial.wTremor4sLFP_dB(idxSubT)/10))) 10*log10(nanmean(10.^(MLMTableTrial.wTremor4sLFP_dB(idxSubN)/10)))];
    scatter(fitVal(end,1),reaVal(end,1),80,'MarkerEdgeColor','none','MarkerFaceColor',[146 197 226]/255,'MarkerFaceAlpha',0.8)
    scatter(fitVal(end,2),reaVal(end,2),80,'MarkerEdgeColor','none','MarkerFaceColor',[244 166 128]/255,'MarkerFaceAlpha',0.8)
end
% legend({'H1','H2','H3','H4','H5','H6','H7','H8','H9'})
fitVal = reshape(fitVal,[],1);
reaVal = reshape(reaVal,[],1);
NonNaN = intersect (find(~isnan(fitVal)),find(~isnan(reaVal)));
fitobject = fit(fitVal,reaVal,'poly1');
h = plot(fitobject,fitVal,reaVal)
set(h(2),'LineWidth',3,'Color','k');
delete(h(1))
hold off
xlabel("Real tremor power (dB)");
ylabel("Fitted tremor power (dB)");
% xlim([min([min(fitVal) min(reaVal)]) max([max(fitVal) max(reaVal)])]);
% ylim([min([min(fitVal) min(reaVal)]) max([max(fitVal) max(reaVal)])]);
[rho,pval] = corr(fitVal,reaVal,'Rows','complete')
title(['R = ' num2str(rho) ', p = ' num2str(pval)]);
set(gca,'FontSize',12);
axis square

function fun_day1day2Baseline()

fileDay1 = {'P6L1','P6R1','P7L1','P7R1','P10L1','P10R1','P11L1','P11R1'};
fileDay2 = {'P6L2','P6R2','P7L2','P7R2','P10L2','P10R2','P11L2','P11R2'};
betaDay = [];
for iFile = 1:length(fileDay1)
    cali1 = load(['D:\MatlabSystem\data\Raw_80\Cali\Trials2s_Session1T4\betaPow\' fileDay1{iFile} '.mat']);
    betaDay(iFile,1) = cali1.powerBetaLFP;
    cali2 = load(['D:\MatlabSystem\data\Raw_80\Cali\Trials2s_Session1T4\betaPow\' fileDay2{iFile} '.mat']);
    betaDay(iFile,2) = cali2.powerBetaLFP;
end
% betaDay = 10*log10(betaDay);
betaDay = 100*(betaDay./mean(betaDay,2)-1)
figure
myColor(1,:) = [180 241 116];
myColor(2,:) = [246 184 11];
myColor(3,:) = [0 0 0];
myColor(4,:) = [35 222 148];
myColor(5,:) = [255 0 0];
myColor(6,:) = [159 89 131];
myColor(7,:) = [246 184 253];
myColor(8,:) = [0 129 255];
myColor = myColor/255;

subplot(3,6,1)
hold on
for iLine = 1:size(betaDay,1)
    line([1 2],[betaDay(iLine,1) betaDay(iLine,2)],'Color',myColor(iLine,:),'Marker','o','LineWidth',3);
end
xlim([0 3])
title('Baseline beta power (%)')
xticks([1 2])
xticklabels({'Day1','Day2'});
legend({'H11','H12','H13','H14','H17','H18','H19','H20'});
set(gca,'box','on','FontSize',15,'LineWidth',2);


subplot(3,6,7)
hBar1 = bar([mean(betaDay(:,1)) 0],'EdgeColor','k','FaceColor',[1 1 1],'LineWidth',3)
hold on
hBar2 = bar([0 mean(betaDay(:,2))],'EdgeColor','k','FaceColor',[1 1 1],'LineWidth',3)
errorbar([1 2],[mean(betaDay(:,1)) mean(betaDay(:,2))],...
    [std(betaDay(:,1))/sqrt(length(betaDay(:,1))) std(betaDay(:,2))/sqrt(length(betaDay(:,2)))],...
    'o','color',[0.5 0.5 0.5],'linewidth',2)
hBar1.BarWidth = 0.4;
hBar2.BarWidth = 0.4;
[p,h] = signrank(betaDay(:,1),betaDay(:,2));
if(p<=0.05)
    H=sigstar({[1,2]},p);
    set(H,'Color','b')
end
xlim([0 3])
ylabel('Diff Between Training and No Training');
xticklabels({'Day1','Day2'});
set(gca,'FontSize',15,'LineWidth',2);
