clear; clc;

cedpath = 'C:\CEDMATLAB\CEDS64ML';
addpath(cedpath);
CEDS64LoadLib(cedpath);

sessionname = 'LAMP0_002';

mvt = 1;
fps = 100;

%% position
optictrack = csvread([sessionname,'.csv'],7,1);
idx = find(optictrack(:,2)==0);
optictrack(idx,2:8) = nan(size(idx,1),7);
optictrack = fillmissing(optictrack,'linear');
t = optictrack(:,1);
x = optictrack(:,6);
y = optictrack(:,8);
vel = sqrt((x(2:end) - x(1:end-1)).^2 + (y(2:end) - y(1:end-1)).^2)*fps;
vel = smoothdata(vel,'movmean',fps*mvt);

%% Open a CED file
% CED recording samplingrate sr1; tick time samplingrate sr2
sr1 = 20000; sr2 = 80000; sr = 100;
fhand1 = CEDS64Open([sessionname,'.smrx']);
maxTimeTicks = CEDS64MaxTime(fhand1) + 1;

% stimulation
ic = 1;
[~, fVals, ~] = CEDS64ReadWaveF(fhand1, ic, maxTimeTicks, 0, maxTimeTicks);
driver = fVals;

% signal
ic = 4;
[fRead, fVals, fTime] = CEDS64ReadWaveF(fhand1, ic, maxTimeTicks, 0, maxTimeTicks);
signal = fVals;

%optitrack sync
ic = 3;
[iRead, vi64T] = CEDS64ReadEvents(fhand1, ic, maxTimeTicks, 0, maxTimeTicks);
dSecs = CEDS64TicksToSecs(fhand1, vi64T );
start = ceil(dSecs(1)); stop = floor(dSecs(end));

% time stamps
tim = CEDS64TicksToSecs(fhand1, [1:maxTimeTicks]');

%% processing
tim = downsample(tim(start*sr2+1:stop*sr2),sr2/sr);
signal = downsample(signal(start*sr1+1:stop*sr1),sr1/sr);
vel = interp1(t(1:end-1)+dSecs(1),vel,tim);
vel = fillmissing(vel,'linear');

Ach.timestamps = tim;
Ach.samplingRate = sr;
Ach.data = signal;
Ach = bz_Filter(Ach,'passband',[0,30],'filter','butter');

AchSignal = smoothdata(detrend(Ach.data),'movmean',sr*mvt);
AchSignal = fillmissing(AchSignal,'linear');

figure();
plot(tim,zscore(vel)); hold on;
plot(tim,zscore(AchSignal)); xlim([start,stop]);
cor = corrcoef(AchSignal,vel);
title(sprintf('Average windows: %d sec.; Corrcoef: %.3f', mvt,cor(1,2)));

savefig(gcf, sprintf('%s_%dsec_smooth',sessionname,mvt));

unloadlibrary ceds64int;