close all; clear all; clc;

T=150*1e-3;  %[s] duration of post-stimulus time window
Fs=1e3;      %[Hz] sampling frequency
Ncs=4;       % Number of channels
Nts=T*Fs;    % Number of timepoints

ts=linspace(0,T,Nts); %[s] time points after stimulus

ft = @(p,t)(p(5)*exp(p(1)*p(3)+0.5*p(2)^2*p(3)^2-p(3)*t)...
    .*normcdf(t,p(1)+p(2)^2*p(3),p(2))+p(4)*normcdf(t,p(1),p(2)));

% Generate some synthetic data with random latency ranging 20÷70 [ms]
mu=zeros(1,Ncs); lt=zeros(1,Ncs); xt=zeros(Ncs,Nts);
for ch=1:Ncs
mu(ch)=(20+randi(50,1));
xt(ch,:)=feval(ft,[mu(ch) 10 3 0.1 6],ts*1e3);
lt(ch)=ts(find(xt(ch,:)>=1/3*max(xt(ch,:)),1,'first'));
end
yt=xt+.01*randn(Ncs,Nts);

% Get Latencies
[latTimes0, lsqCurveTimes0, lsqFittedCurve0, lsqFittedParam0, lsqCurveFitErr0] = computeLatencies(xt,ts*1e3,10);
[latTimes, lsqCurveTimes, lsqFittedCurve, lsqFittedParam, lsqCurveFitErr] = computeLatencies(yt,ts*1e3,10);

%% Plot results
% multi-line plot offset
ploffs0=repmat((1:Ncs)'*max(yt(:)),1,Nts);
ploffs=repmat((1:Ncs)'*max(yt(:)),1,length(lsqCurveTimes));

figure(); hold on
plot(-inf,-inf,'b-'); plot(-inf,-inf,'k-'); plot(-inf,-inf,'r-');
plot(-inf,-inf,'bo'); plot(-inf,-inf,'ko'); plot(-inf,-inf,'ro');
plot(-inf,-inf,'bx'); plot(-inf,-inf,'kx'); plot(-inf,-inf,'rx');
plot(ts*1e3,xt+ploffs0,'b-'); hold on;
plot(ts*1e3,yt+ploffs0,'k-'); hold on;
plot(lsqCurveTimes,lsqFittedCurve+ploffs,'r-'); hold on;

plot(lt*1e3,ploffs0(:,1),'bo');
plot(latTimes0,ploffs0(:,1),'ko');
plot(latTimes,ploffs0(:,1),'ro');

plot(mu,ploffs0(:,1),'bx');
plot(lsqFittedParam0(:,1),ploffs0(:,1),'kx');
plot(lsqFittedParam(:,1),ploffs0(:,1),'rx');
xlabel('time (ms)');
ylabel('channels'); set(gca,'YTick',[]);

hl=legend('x(t) sample signal','y(t)=x(t)+noise','f(t) fitted function','true latency',...
    'latency fit(x)','latency fit(y)','true mean','mean fit(x)','mean fit(y)');
hl.EdgeColor='none'; hl.Location='southeast';