%% legacy file: Apub_NSx_CohPSD_v3.m
% Version b6: compatible with openephys
% Version 6: Automatic add subfolders
% Version 5: forced down sample to 1000 Hz; Solve memory and graphic problem
% Version 4: change from NS2 to NSx

clear all; close all;

%% Input Parameter %%%%%%%%%%%%%%%%%%
CallExcel=0; % Use spreadsheet instead of data here
DataSource=1; % 1:BlackRock; 2:OpenEphys

ChSig1=10;  
ChRef1=11;
ChSig2=129;
ChRef2=0;
timeframe=20; 
CohSec=1;   
BW=25.4; 
Timecut=0; % trim file
  TrimStartTime=0; 
  TrimStopTime=1300; 
Timeshift= 0; % in sec, shift (ChSig2-ChRef2), can be negative

Filter=1; % if =1, apply filter
 ChannelNotFilter=[137]; % assigned channel those do not apply filter
   Lowcut=0.3;
   Highcut=250;
%% Main function %%%%%%%%%%
if CallExcel==1;
A=xlsread('Input parameters3.xlsx','B2:E2');
ChSig1=A(1); % 指定 Signal1 的 Channel
ChSig2=A(3); % 指定 Siganl2 的 Channel
ChRef1=A(2) ; % 指定 Reference1 的 Channel 不要的話，設為０
ChRef2=A(4); % 指定 Reference2 的 Channel 不要的話，設為０
clear A;
else
end
%%%%%%%%% Put input parameters into structure %%%%%%%%%%%%%
save temp.mat;
a_InputParameters=load('temp.mat');
delete temp.mat;
%%%%%%%%%% Check error %%%%%%%%%%%%%%%%%%%%
if Highcut >250;
    error('Highcut is above 250 Hz');
end
%%%%%%%%  Load NSx data (BlackRock of OpenEphys) %%%%%%%%%%%%%%%%%%%%%%
if DataSource==1
[GetFileName GetFilePath]=uigetfile('*.ns*'); % get filename and path
NSx=openNSx([GetFilePath GetFileName]);
end
if DataSource==2
[NSx,EphysToMat]=Afun_OpenephysToNSx_v2(0,0,0);
GetFilePath=[pwd '\']; % pwd is matlab intrinsic code for current path
end
SourceFile=[NSx.MetaTags.Filename NSx.MetaTags.FileExt];
data=double(NSx.Data');
NSx=rmfield(NSx,'Data'); 
a_InputParameters.SourceFile=SourceFile;
fs=double(NSx.MetaTags.SamplingFreq);
a_InputParameters.fs=fs;

%%%%%%%% Create subfolders to store file %%%%%%%%
SubPath=sprintf('Ch%d-%d_Ch%d-%d_fil%d',ChSig1,ChRef1,ChSig2,ChRef2,Filter); % sprintf did not allow '\'
SubPath=[SubPath,'\']; 
FullFilePath=[GetFilePath,SubPath];
mkdir(FullFilePath); % make the folder (cannot write if the folder do not exist)
%%%%%%%% Re-asign channelID to legacy variables %%%%%%%%%%
ChannelSignal1=find(NSx.MetaTags.ChannelID==ChSig1);
ChannelSignal2=find(NSx.MetaTags.ChannelID==ChSig2);
if ChRef1==0; ChannelRef1=0; 
else ChannelRef1=find(NSx.MetaTags.ChannelID==ChRef1);
end
if ChRef2==0; ChannelRef2=0; 
else ChannelRef2=find(NSx.MetaTags.ChannelID==ChRef2);
end
Title1=[sprintf('%d-%d vs %d-%d ',ChSig1,ChRef1,ChSig2,ChRef2),SourceFile];

%%% Trim data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if Timecut~=0;
data=data(TrimStartTime*fs+1:TrimStopTime*fs,:);
Title1=['C*',Title1];
end;
   
%%% Assign channel %%%
if ChannelRef1==0;
    Signal1=data(:,ChannelSignal1);
else
    Signal1=data(:,ChannelSignal1)-data(:,ChannelRef1);
end
if ChannelRef2==0;
    Signal2=data(:,ChannelSignal2);
else
    Signal2=data(:,ChannelSignal2)-data(:,ChannelRef2);
end    
clear data;
%%% Filtering Data (optional) %%%%%%%%%%%

if Filter==1;
    Title1=['F*',Title1];
   if max(ismember([ChSig1,ChRef1],ChannelNotFilter))==0; % filter if ChSig1 and ChRef1 both not in the "nofilter channel" 
        Signal1=A_applyfilter(Signal1,fs,Lowcut,Highcut);
   else a_InputParameters.note(1)={'signal1 not filtered'};
   end
   if max(ismember([ChSig2,ChRef2],ChannelNotFilter))==0;
        Signal2=A_applyfilter(Signal2,fs,Lowcut,Highcut);
   else a_InputParameters.note(2)={'signal2 not filtered'};
   end
end
%%% Down sampling data %%%
if fs>1000;
    DownSamplePoints=1:fs/1000:length(Signal1);
    Signal1=(Signal1(DownSamplePoints));
    Signal2=(Signal2(DownSamplePoints));
    fs=1000;
end
%%%%% Calculate Coherence and PSD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% calculation
tic
N = fs*timeframe; 
freqStep = fs/N; 
N1s=fs;          
NCoh=N1s*CohSec; 
FreqCohRes=fs/NCoh;
EndNumber= fix(length(Signal1)/N1s)-timeframe-ceil(abs(Timeshift));
Shiftedframe= N1s*Timeshift;
if Shiftedframe<0
    AntiShiftframe=Shiftedframe*(-1);
    Shiftedframe=0;
else
    AntiShiftframe=0;
end

if Timeshift~=0
    Title1=['TS*',Title1];
end
    

for s=1:EndNumber
    
y1=Signal1((1+(s-1)*N1s+AntiShiftframe):(N+(s-1)*N1s)+AntiShiftframe);   
y2=Signal2((1+(s-1)*N1s+Shiftedframe):(N+(s-1)*N1s+Shiftedframe)); 

[Cy1y2,FreqCoh]=mscohere(y1,y2,hanning(NCoh),1/2*NCoh,NCoh,fs); % (y1,y2 coherece, with hamming window, take NCoh points, half points overlap, number of fft:NCoh, sampling frequency:fs
Coh(:,s)=Cy1y2; 
FreqLable(:,s)=FreqCoh;
TimeLable(1:(NCoh/2+1),s)=s-1;

[PxxY1,freqPsY1]=pwelch(y1,hanning(NCoh),1/2*NCoh,NCoh,fs); % PSD by Wilch method and Hanning window
[PxxY2,freqPsY2]=pwelch(y2,hanning(NCoh),1/2*NCoh,NCoh,fs);
PSDY1(:,s)=PxxY1;
PSDY2(:,s)=PxxY2;

end
a_InputParameters.Title1=Title1;
save([FullFilePath,'Coherence_PSD_Profiles&Results.mat'],'Coh','PSDY1','PSDY2','FreqCoh','a_InputParameters','NSx','TimeLable','FreqLable','Title1');
toc
%%  Plotting  %%%%%%%%%%

%figure(1);
%gca=pcolor(TimeLable,FreqLable,Coh);
%set(gca, 'LineStyle','none');
%colorbar('location','eastoutside');
%title(['Cohere ',Title1]); 
%xlabel('Time (seconds)'); 
%ylabel('Frequency (Hz)'); 
%saveas(1,'Coherence_0to250Hz.fig');

%figure(2);
%gca=pcolor(TimeLable,FreqLable,Coh);
%set(gca, 'LineStyle','none');
%colorbar('location','eastoutside');
%title(['Cohere ',Title1]); 
%xlabel('Time (seconds)'); 
%ylabel('Frequency (Hz)'); 
%axis([-inf,inf,0,100]);
%saveas(2,'Coherence_0to100Hz.fig');

%figure(3);
%gca=pcolor(TimeLable,FreqLable,Coh);
%set(gca, 'LineStyle','none');
%colorbar('location','eastoutside');
%title(['Cohere ',Title1]); 
%xlabel('Time (seconds)'); 
%ylabel('Frequency (Hz)'); 
%axis([-inf,inf,0,50]);
%saveas(3,'Coherence_0to50Hz.fig');

%figure(6);
%gca3=pcolor(TimeLable,FreqLable,PSDY1);
%set(gca3, 'LineStyle','none');
%colorbar('location','eastoutside');
%caxis([0 1000]);
%title(['PSD Sig1 ',Title1]); 
%xlabel('Time (seconds)'); 
%ylabel('Frequency (Hz)'); 
%axis([-inf,inf,0,100]);
%saveas(6,'PSD_Singal1_0to100Hz.fig');

%figure(8);
%gca3=pcolor(TimeLable,FreqLable,PSDY2);
%set(gca3, 'LineStyle','none');
%colorbar('location','eastoutside');
%caxis([0 1000]);
%title(['PSD Sig2 ',Title1]); 
%xlabel('Time (seconds)'); 
%ylabel('Frequency (Hz)');
%axis([-inf,inf,0,100]);
%saveas(8,'PSD_Signal2_0to100Hz.fig');
%close all;

figure(9);
subplot(3,1,1);
gca3=pcolor(TimeLable,FreqLable,PSDY1);
set(gca3, 'LineStyle','none');
colorbar('location','eastoutside');
caxis([0 1000]);
title([Title1]); 
ylabel('PSDY1(Hz)'); 
axis([-inf,inf,0,50]);

subplot(3,1,2);
gca3=pcolor(TimeLable,FreqLable,PSDY2);
set(gca3, 'LineStyle','none');
colorbar('location','eastoutside');
caxis([0 1000]);
ylabel('PSDY2(Hz)');
axis([-inf,inf,0,50]);

subplot(3,1,3);
gca=pcolor(TimeLable,FreqLable,Coh);
set(gca, 'LineStyle','none');
colorbar('location','eastoutside');
xlabel('Time (seconds)'); 
ylabel('Coh(Hz)'); 
axis([-inf,inf,0,50]);
saveas(9,[FullFilePath,'TFPLot_PSDY1Y2+Coh.fig']);
saveas(9,[FullFilePath,'TFPLot_PSDY1Y2+Coh.jpg']);

