%% 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

%% Function version of  Apub_NSx_CohPSD_vb6.m
% 2022.05.30: path: directory of the target file/ file: filename/ para: input parameters

function F_NSx_CohPSD(path,file,para)
%% Input Parameter %%%%%%%%%%%%%%%%%%
DataSource=para.DataSource; % 1:BlackRock

ChSig1=para.ChSig1;  
ChRef1=para.ChRef1;
ChSig2=para.ChSig2;
ChRef2=para.ChRef2;
timeframe=para.timeframe; 
CohSec=para.CohSec;   
BW=para.BW;

Timecut=para.Timecut; % trim file
  TrimStartTime=para.TrimStartTime; 
  TrimStopTime=para.TrimStopTime; 
Timeshift=para.Timeshift; % in sec, shift (ChSig2-ChRef2), can be negative

Filter=para.Filter; % if =1, apply filter
 ChannelNotFilter=para.ChannelNotFilter; % assigned channel those do not apply filter
   Lowcut=para.Lowcut;
   Highcut=para.Highcut;
%%%%%%%%%% Check error %%%%%%%%%%%%%%%%%%%%
if Highcut >250;
    error('Highcut is above 250 Hz');
end
%%%%%%%%  Load NSx data (BlackRock of OpenEphys) %%%%%%%%%%%%%%%%%%%%%%
GetFilePath=[path,'\'];
NSx=openNSx([GetFilePath,file]);
SourceFile=[NSx.MetaTags.Filename NSx.MetaTags.FileExt];
data=double(NSx.Data');
NSx=rmfield(NSx,'Data'); 
para.SourceFile=SourceFile;
fs=double(NSx.MetaTags.SamplingFreq);
para.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 para.note(1)={'signal1 not filtered'};
   end
   if max(ismember([ChSig2,ChRef2],ChannelNotFilter))==0;
        Signal2=A_applyfilter(Signal2,fs,Lowcut,Highcut);
   else para.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
para.Title1=Title1;
save([FullFilePath,'Coherence_PSD_Profiles&Results.mat'],'Coh','PSDY1','PSDY2','FreqCoh','para','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,60]);

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,60]);

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,60]);
saveas(9,[FullFilePath,'TFPLot_PSDY1Y2+Coh.fig']);
saveas(9,[FullFilePath,'TFPLot_PSDY1Y2+Coh.jpg']);
toc
end