% Detect EP trigger points for chirp waves (local maximum)
% Step1: load .ns6 file
% Step2: load proper Chirp File.mat (should be correct in fs, frequencies and duration)
% Step3: output EPindex as the correct EP detect points. 

clear all; close all;
%% input parameters
Th=1000;
PlotFigure=1; % if ==1, plot figure (will be slow)
ChSig1=[136]; % EP channel
ChRef1=[0]; % if ==0, no reference.

EPChannel1=1; % if==1, EPChannelRaw=Signal1; ==2, EPchannelRaw=Signal2;

%%%%%%%%% Put input parameters into structure %%%%%%%%%%%%%
save temp.mat;
a_InputParameters=load('temp.mat');
delete temp.mat;
%% load data
[GetFileName GetFilePath]=uigetfile('*.ns*'); % get filename and path
NSx=openNSx([GetFilePath GetFileName]);

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;
%% Re-asign channelID to legacy variables %%%%%%%%%%
 ChannelSignal1=find(NSx.MetaTags.ChannelID==ChSig1);
 if ChRef1==0; ChannelRef1=0; 
 else ChannelRef1=find(NSx.MetaTags.ChannelID==ChRef1);
 end
%% Assign channel %%%
    if ChannelRef1==0;
        Signal1=data(:,ChannelSignal1);
    else
        Signal1=data(:,ChannelSignal1)-data(:,ChannelRef1);
    end
    clear data;


%% load ChirpWaves
[C,D]=uigetfile('Chirp*.mat');
load([D,C],'ChirpWave');
a_InputParameters.ChirpFileSource=C;
clear C D;

%% Add threshold and find first chirp points
EPChannelRaw=Signal1;
clear Signal1;
A=find(EPChannelRaw<=Th);
EPChannelRaw(A)=0;
clear A;



=min(find(EPChannelRaw>0));
SigChirp(1:length(EPChannelRaw))=0;
SigChirp(X:X+length(ChirpWave)-1)=ChirpWave;
[A,EPindex]=findpeaks(SigChirp);


save('EPindex_ChirpWave','NSx','a_InputParameters','ChirpWave','SigChirp','EPindex','-v7.3');


%% plot
Tvec=0:1:(length(EPChannelRaw)-1);
%Tvec=Tvec/fs;

if PlotFigure==1;
   figure;
    H1=plot(Tvec,EPChannelRaw);

    figure;
    H2=plot(Tvec,[EPChannelRaw';SigChirp*3000]);

end