%% Legacy file: Apub_NS2_Extract2_v3.m
% Version 13: set legend off on Subplot 1,2 and label all details in all legends
% Version 12: Able to do sequential ExtEventIndex based on SaveFileOrder=0
%    if SaveFileOrder==0, ExtEventIndex=ExtEventIndex0
%    if SaveFileOrder==1, ExtEventIndex=ExtEventIndex1
%     and so on
% Version 11: Add self normalization (to its own window period),& corresponding a_ExtractSummary and Excel file
% Version 10: Make legend order (v9) optional (% LegendOrder==1, has legend order) 
% Version 9: Add legend order (number) to each events in the graphic legend
% Version 8: Allow save additional mat and fig files. (e.g.Mean_PSD_new2e1.mat)
% Version 7: Support Extractwindow from A_OnOffExtract_byTFPlot
% Version 6: Support reading EventListRaw2 (which is manually changed variable stored in EventList.mat)
% Version 5: support reading EventList from NicoletEEG extract, also correct ExtractWindows in a_ExtraPara to correct time
% Version 4: use Title1 to more clearly identify filter/timeshift/timecut issues
clear all; close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Input Parameters %
CallExcel=0;
LegendOrder=1; % if ==1, legend order (+), else negative

SaveFileOrder=0; % If ==0, generate "Mean_PSD_new2.mat"; if any positive number, e.g. if the number is 1, generate "Mean_PSD_new2e1.mat" 
LoadEventList=0; % if==0, use ExtractWindows assigned below, no SubBoundary
                 % if 1, load event list from EventList.mat, % input ExtEventIndex below
                 % if 2, load EventListRaw2 (mannualy corrected version)
                 % if 11, load EventListTFPlot.mat, and find EventWindow (automatic) %
                 % if 12, load EventListTFPlot.mat, and find EventWindow2 (manually assinged)
                 % if 111, include both ExtractWindows (list) and loaded EventWindow
  ExtEventIndex0=[4,5,6,8,9,10,12,13,14,15,16,18,19,20,21,22,24,25,26,27,38,28,39,41,44,47,50,54,57,59,60,63,66,70,73,76,77]; %ALL % Input raw index to extract from EventListRaw table, can be multiple 
  ExtEventIndex1=[4,5,9,13,14,15,19,20,21,25,26,27,39,47,57,76,77];  %e1: Rest,kinetic, posture, 500g % Input raw index to extract from EventListRaw table, can be multiple
  ExtEventIndex2=[4,5,6,10,16,22,38,41,50,59,63,70]; %e2: right 
  ExtEventIndex3=[4,5,8,12,18,24,28,44,54,60,66,73]; %e3: left 
  ExtEventIndex4=[4,5,6,10,41,50,59,63,70]; %e2: right, no 4,5,6 Hz 
  ExtEventIndex5=[4,5,8,12,44,54,60,66,73]; %e3: left, no 4,5,6 Hz
  %  ExtEventIndex4=[3,4,8,12,19,26,32,5,7,9,11,34,38]; %e4: optional 
%  ExtEventIndex5=[3,4,7,11,18,25,31,38]; %e5: optional 
%  ExtEventIndex6=[3,4,7,11,18,25,31,38]; %e6: optional 
  
%  ExtEventIndex=[5,12,17,21,27]; % Input raw index to extract from EventListRaw table, can be multiple 
SubBoundary=[1,1]; % Adjust boundary for ExtractWindows. Subtract column 1 (1 sec) for event begin, column 2 (1 sec) for event end)
  
ExtractFreqRange=[5,40]; % Define the frequency range to extract the data
ExtractWindows =[1,180;283,392;673,782;1063,1172;1453,1562;1843,1952;2233,2342;2623,2732]; % Define the extract windows (can be multiple)
%
nSig1PlotLimit=[0,60,0,15];
nSig2PlotLimit=[0,60,0,3];

SumWindow=[40,50]; % Define the frequency range for normalize the data 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% load excel parameters %
if CallExcel==1;
A=xlsread('Input parameters3.xlsx','G2:O2');
ExtractFreqRange=A(8:9); % Define the frequency range to extract the data
ExtractWindows =[A(1:2);A(3:4);A(5:6)]; % Define the extract windows (can be multiple)
else end;

save temp.mat;
a_ExtractPara=load('temp.mat');
delete temp.mat;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load Coherence_PSD_Profiles&Results.mat PSDY1 PSDY2 Coh FreqCoh a_InputParameters;
FreqCohRes=1/a_InputParameters.CohSec; % Defined Frequency resolution
Title1=a_InputParameters.Title1;
%%%%%%%%%%%%%%%%%%%%% Consistency check-up %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if FreqCohRes~=FreqCoh(2);
    errordlg('frequency resolution mismatch');
    return
end
%% Since Version 12, Sequential ExtEventIndex %%
    Temp01=sprintf('ExtEventIndex%d',SaveFileOrder);
    eval(['ExtEventIndex=',Temp01,';']);
    clear Temp01;
%% Load EventList data for events and legends
 % if LoadEventList not assinged, will use input parameters
    
if LoadEventList==1
    load('EventList.mat');
    Events=EventListRaw(ExtEventIndex,:);
    ExtractWindows=cell2mat(Events{:,4:5});
    ExtractWindows(:,1)=ExtractWindows(:,1)+SubBoundary(1);
    ExtractWindows(:,2)=ExtractWindows(:,2)-SubBoundary(2);
    a_ExtractPara.ExtractWindows=ExtractWindows;
end
if LoadEventList==2
    load('EventList.mat');
    Events=EventListRaw2(ExtEventIndex,:);
    ExtractWindows=cell2mat(Events{:,4:5});
    ExtractWindows(:,1)=ExtractWindows(:,1)+SubBoundary(1);
    ExtractWindows(:,2)=ExtractWindows(:,2)-SubBoundary(2);
    a_ExtractPara.ExtractWindows=ExtractWindows;
end

if LoadEventList==11
    load('EventListTFPlot.mat','EventWindow');
    ExtractWindows=EventWindow;
    ExtractWindows(:,1)=ExtractWindows(:,1)+SubBoundary(1);
    ExtractWindows(:,2)=ExtractWindows(:,2)-SubBoundary(2);
    a_ExtractPara.ExtractWindows=ExtractWindows;
end

if LoadEventList==111
    load('EventListTFPlot.mat','EventWindow');
    ExtractWindows=[ExtractWindows;EventWindow];
    ExtractWindows(:,1)=ExtractWindows(:,1)+SubBoundary(1);
    ExtractWindows(:,2)=ExtractWindows(:,2)-SubBoundary(2);
    a_ExtractPara.ExtractWindows=ExtractWindows;
end

% Asign legend
    AB=cell(size(ExtractWindows,1),1);
    K=cell(size(ExtractWindows,1),1);
    for i=size(ExtractWindows,1):-1:1
        AB{i,1}=sprintf('%d-%d %d',ExtractWindows(i,1),ExtractWindows(i,2));
      if LegendOrder==1;
        K{i,1}=sprintf('%d- ',i);
      else
        K{i,1}=[]; % create empty cell
      end
    end
    
if LoadEventList==1 | LoadEventList==2 % '|' means or
    Legend=strcat(K,Events{:,1},{' '},AB);
else 
    Legend=strcat(K,{'   '},AB);
end
clear AB K;

ExtractWindows(:,3) = ExtractWindows(:,2)-a_InputParameters.timeframe; % Define the Endtime by subtract the timeframe (for internal consistency of previous version)
[WindowNumber,b]=size(ExtractWindows); % Define the extraction size
clear b;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
meanPSDY1(:,1)=mean(PSDY1(:,ExtractWindows(1,1)+1:ExtractWindows(1,3)),2);
meanPSDY2(:,1)=mean(PSDY2(:,ExtractWindows(1,1)+1:ExtractWindows(1,3)),2);
meanCoh(:,1)=mean(Coh(:,ExtractWindows(1,1)+1:ExtractWindows(1,3)),2);

SumWindowIndex=SumWindow/FreqCohRes+1; % calculating the normalizing value from the SumWindow)
sumPSDY1=sum(meanPSDY1(SumWindowIndex(1):SumWindowIndex(2),1));
sumPSDY2=sum(meanPSDY2(SumWindowIndex(1):SumWindowIndex(2),1));
%sumCoh=sum(meanCoh(SumWindowIndex(1):SumWindowIndex(2),1));

if WindowNumber>=2;
    for i=WindowNumber:-1:1;
     meanPSDY1(:,i)=mean(PSDY1(:,ExtractWindows(i,1)+1:ExtractWindows(i,3)),2);
     meanPSDY2(:,i)=mean(PSDY2(:,ExtractWindows(i,1)+1:ExtractWindows(i,3)),2);
     meanCoh(:,i)=mean(Coh(:,ExtractWindows(i,1)+1:ExtractWindows(i,3)),2);
     sumPSDY1Self(i)=sum(meanPSDY1(SumWindowIndex(1):SumWindowIndex(2),i));
     sumPSDY2Self(i)=sum(meanPSDY2(SumWindowIndex(1):SumWindowIndex(2),i));
     nmeanPSDY1Self(:,i)=meanPSDY1(:,i)/sumPSDY1Self(i);
     nmeanPSDY2Self(:,i)=meanPSDY2(:,i)/sumPSDY2Self(i);    
    end
end
nmeanPSDY1=meanPSDY1/sumPSDY1;
nmeanPSDY2=meanPSDY2/sumPSDY2;

%%%%%%%%%%%%%%%%%% Find max and freqeuncy %%%%%%%%%%%%%%%%%%%%%%%%
ExtractFreqIndex=ExtractFreqRange/FreqCohRes+1;
[maxPSDY1,indexPSDY1]=max(meanPSDY1(ExtractFreqIndex(1):ExtractFreqIndex(2),:));
nmaxPSDY1=maxPSDY1/sumPSDY1;
[maxPSDY2,indexPSDY2]=max(meanPSDY2(ExtractFreqIndex(1):ExtractFreqIndex(2),:));
nmaxPSDY2=maxPSDY2/sumPSDY2;
[maxCoh,indexCoh]=max(meanCoh(ExtractFreqIndex(1):ExtractFreqIndex(2),:));
FreqmaxPSDY1=(ExtractFreqIndex(1)+indexPSDY1-2)*FreqCohRes;
FreqmaxPSDY2=(ExtractFreqIndex(1)+indexPSDY2-2)*FreqCohRes;
FreqmaxCoh=(ExtractFreqIndex(1)+indexCoh-2)*FreqCohRes;
nmaxPSDY1Self=maxPSDY1./sumPSDY1Self; % find normalized PSDY1 peak with its own normalization method.
nmaxPSDY2Self=maxPSDY2./sumPSDY2Self; % find normalized PSDY2 peak with its own normalization method.


%%%%%%%%%%%%%%%%%% Find max according to baseline freqeuncy %%%%%%%%%%%%%%%%%%%%%%%%
nmaxPSDY1BaseFreq=meanPSDY1(ExtractFreqIndex(1)+indexPSDY1(1)-1,:)/sumPSDY1;
nmaxPSDY2BaseFreq=meanPSDY2(ExtractFreqIndex(1)+indexPSDY2(1)-1,:)/sumPSDY2;
maxCohBaseFreq=meanCoh(ExtractFreqIndex(1)+indexCoh(1)-1,:);
nmaxPSDY1BaseFreqSelf=meanPSDY1(ExtractFreqIndex(1)+indexPSDY1(1)-1,:)./sumPSDY1Self; % Normalized according to its own time window
nmaxPSDY2BaseFreqSelf=meanPSDY2(ExtractFreqIndex(1)+indexPSDY2(1)-1,:)./sumPSDY2Self; % Normalized according to its own time window

%a_ExtractSummary=cell(14,WindowNumber+1);
 a_ExtractSummary=cell(28,WindowNumber+1);
%a_ExtractSummary(:,1)={'Window','nPSDY1','nPSDY2','Coh','Y1Frq','Y2Freq','CohFreq','','nPSDY1_bfreq','nPSDY2_bfreq','Coh_bfreq','Y1Frq','Y2Freq','CohFreq'};
 a_ExtractSummary(:,1)={'Window','nPSDY1','nPSDY2','Coh','Y1Frq','Y2Freq','CohFreq','','nPSDY1_bfreq','nPSDY2_bfreq','Coh_bfreq','Y1Frq','Y2Freq','CohFreq','','nSelfPSDY1','nSelfPSDY2','Coh','Y1Frq','Y2Freq','CohFreq','','nSelfPSDY1self_bfreq','nSelfPSDY2self_bfreq','Coh_bfreq','Y1Frq','Y2Freq','CohFreq'};

for i=WindowNumber:-1:1
%    a_ExtractSummary(:,i+1)={Legend{i,1},nmaxPSDY1(i), nmaxPSDY2(i), maxCoh(i),FreqmaxPSDY1(i),FreqmaxPSDY2(i),FreqmaxCoh(i),'',nmaxPSDY1BaseFreq(i), nmaxPSDY2BaseFreq(i), maxCohBaseFreq(i),FreqmaxPSDY1(1),FreqmaxPSDY2(1),FreqmaxCoh(1)};
     a_ExtractSummary(:,i+1)={Legend{i,1},nmaxPSDY1(i), nmaxPSDY2(i), maxCoh(i),FreqmaxPSDY1(i),FreqmaxPSDY2(i),FreqmaxCoh(i),'',nmaxPSDY1BaseFreq(i), nmaxPSDY2BaseFreq(i), maxCohBaseFreq(i),FreqmaxPSDY1(1),FreqmaxPSDY2(1),FreqmaxCoh(1),'',nmaxPSDY1Self(i),nmaxPSDY2Self(i),maxCoh(i),FreqmaxPSDY1(i),FreqmaxPSDY2(i),FreqmaxCoh(i),'',nmaxPSDY1BaseFreqSelf(i),nmaxPSDY2BaseFreqSelf(i),maxCohBaseFreq(i),FreqmaxPSDY1(1),FreqmaxPSDY2(1),FreqmaxCoh(1)};

end

clear Coh PSDY1 PSDY2 i;
if SaveFileOrder==0
   save Mean_PSD_new2.mat;
else 
    FF=sprintf('Mean_PSD_new2e%d.mat',SaveFileOrder);
    save(FF);
end
%%
%%%%%%  Plot %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
subplot(3,1,1);
h11=plot(FreqCoh ,meanPSDY1);
title(['PSD (Sig1) ' Title1]);
xlabel('Frequency (Hz)'); 
ylabel('PSD (dB/Hz)');
axis([nSig1PlotLimit(1:3),max(maxPSDY1)*1.2]);
legend(Legend);
legend('off');
%set(h11(1),'Color',[0,0,1]); % Set the line 1 by RGB color [0 0 1] (blue)

subplot(3,1,2);
plot(FreqCoh ,meanPSDY2);
title('PSD (Sig2)'); 
xlabel('Frequency (Hz)'); 
ylabel('PSD (dB/Hz)');
axis([nSig2PlotLimit(1:3),max(maxPSDY2)*1.2]);
legend(Legend);
legend('off');

subplot(3,1,3);
plot(FreqCoh ,meanCoh);
title('PSD (Coh)'); 
xlabel('Frequency (Hz)'); 
ylabel('Coh');
xlim([nSig2PlotLimit(1),nSig2PlotLimit(2)]);
legend(Legend);

if SaveFileOrder==0
   saveas(1,'PSD2.fig');
   saveas(1,'PSD2.jpg');
else 
    FF1=sprintf('PSD2e%d.fig',SaveFileOrder);
    saveas(1,FF1);
    FF2=sprintf('PSD2e%d.jpg',SaveFileOrder);
    saveas(1,FF2);
end
%%
figure(2);
subplot(3,1,1);
h21=plot(FreqCoh ,nmeanPSDY1);
title(['nPSD (Sig1) ' Title1]);
xlabel('Frequency (Hz)'); 
ylabel('PSD (dB/Hz)');
axis(nSig1PlotLimit);
legend(Legend);
legend('off');
set(h21(1),'Color',[0,0,1]); % Set the line 1 by RGB color [0 0 1] (blue)

subplot(3,1,2);
h22=plot(FreqCoh ,nmeanPSDY2);
title('nPSD (Sig2)'); 
xlabel('Frequency (Hz)'); 
ylabel('PSD (dB/Hz)');
axis(nSig2PlotLimit);
set(h22(1),'Color',[0,0,1]); % Set the line 1 by RGB color [0 0 1] (blue)
legend(Legend);
legend('off');

subplot(3,1,3);
h23=plot(FreqCoh ,meanCoh);
title('nPSD (Coh)'); 
xlabel('Frequency (Hz)'); 
ylabel('Coh');
xlim([nSig2PlotLimit(1),nSig2PlotLimit(2)]);
set(h23(1),'Color',[0,0,1]); % Set the line 1 by RGB color [0 0 1] (blue)
legend(Legend);

if SaveFileOrder==0
   saveas(2,'nPSD2.fig');
   saveas(2,'nPSD2.jpg');
else 
    FF1=sprintf('nPSD2e%d.fig',SaveFileOrder);
    saveas(2,FF1);
    FF2=sprintf('nPSD2e%d.jpg',SaveFileOrder);
    saveas(2,FF2);
end
%%
figure(3);
subplot(3,1,1);
h21=plot(FreqCoh ,nmeanPSDY1Self);
title(['nPSDSelf (Sig1) ' Title1]);
xlabel('Frequency (Hz)'); 
ylabel('PSD (dB/Hz)');
axis(nSig1PlotLimit);
legend(Legend);
legend('off');
set(h21(1),'Color',[0,0,1]); % Set the line 1 by RGB color [0 0 1] (blue)

subplot(3,1,2);
h22=plot(FreqCoh ,nmeanPSDY2Self);
title('nPSDSelf (Sig2)'); 
xlabel('Frequency (Hz)'); 
ylabel('PSD (dB/Hz)');
axis(nSig2PlotLimit);
set(h22(1),'Color',[0,0,1]); % Set the line 1 by RGB color [0 0 1] (blue)
legend(Legend);
legend('off');

subplot(3,1,3);
h23=plot(FreqCoh ,meanCoh);
title('nPSD (Coh)'); 
xlabel('Frequency (Hz)'); 
ylabel('Coh');
xlim([nSig2PlotLimit(1),nSig2PlotLimit(2)]);
set(h23(1),'Color',[0,0,1]); % Set the line 1 by RGB color [0 0 1] (blue)
legend(Legend);

if SaveFileOrder==0
   saveas(3,'nPSD2_self.fig');
   saveas(3,'nPSD2_self.jpg');
else 
    FF1=sprintf('nPSD2e%d_self.fig',SaveFileOrder);
    saveas(3,FF1);
    FF2=sprintf('nPSD2e%d_self.jpg',SaveFileOrder);
    saveas(3,FF2);
end


%%
%%%%%%%%%% Export Data to Excel %%%%%%%%%%%%%%%%%%%%%%%5
FF3=sprintf('PSDresultse%d.xls',SaveFileOrder);
xlswrite(FF3,a_ExtractSummary,'Sheet1','A4');
xlswrite(FF3,{'FileName',a_InputParameters.SourceFile},'Sheet1','A1');
xlswrite(FF3,{'Title',Title1},'Sheet1','A2');
xlswrite(FF3,meanPSDY1,'meanPSDY1','A2');
xlswrite(FF3,meanPSDY2,'meanPSDY2','A2');
xlswrite(FF3,meanCoh,'meanCoh','A2');
xlswrite(FF3,nmeanPSDY1,'nmeanPSDY1','A2');
xlswrite(FF3,nmeanPSDY2,'nmeanPSDY2','A2');
xlswrite(FF3,nmeanPSDY1Self,'nmeanPSDY1Self','A2');
xlswrite(FF3,nmeanPSDY2Self,'nmeanPSDY2Self','A2');