%calculate the offset between Logitech video and Spike2 audio
%this code calculates the cross-correlation between the high-pass filtered
%webcam audio and the low-pass filtered spike2 audio to determine the
%offset between these two data streams

% (1) before running this code, scroll through audio and find a time with low
%frequency noise to use for alignment
%(2) load video, grab audio from webcam video, and high pass filter
%(3) load your high quality audio recording, low pass filter, and
%downsample to match sample rate of webcam audio
%(4) run a cross-correlation on the two audio signals
%(5) get the time lag of maximum correlation

close all; clear all;
% filename='VM8 final opto.wmv';
filename='Video 1760.wmv'; %won't work on mp4 files!!!
filename2='50to150.mat';
channel='opto_Ch1';

spike2SampleRate=250000;
%video is about 3s delayed relative to Spike2
vidTimeToGrab=[50 150]; %specified in seconds


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% RUN AFTER SETTING PARAMS!!!

%open webcam video and audio
% [y,Fs]=audioread(filename);
[videoGrabbed, audioGrabbed] = mmread(filename,[],vidTimeToGrab);
soundData=audioGrabbed.data;
soundData=soundData(:,1);
soundRate=audioGrabbed.rate; %48000 Hz
soundDuration=length(soundData)/soundRate;
x=1:length(soundData);
x=x-x(1);
x=x/soundRate; %make a time base
%high pass filter soundData, over 300 Hz
N=2;
Wn=8000/48000;
% Wn=500/48000; %Wn is the cutoff frequency, normalized from 0 to 1, 
%1 corresponds to half the sample rate of the signal
[B,A] = butter(N,Wn,'high');
soundDataFiltered= filter(B,A,soundData);
% soundDataFiltered=soundData;

%load spike2 audio
spike2=load(filename2);
spike2=spike2.(channel);
spike2=spike2.values;
spike2Duration=length(spike2)/spike2SampleRate; %in seconds
xoriginal=1:length(spike2);
xoriginal=xoriginal-xoriginal(1);
xoriginal=xoriginal/spike2SampleRate;

%low pass filter spike 2 data
N2=2;
Wn2=(soundRate/2)/spike2SampleRate; %Wn is the cutoff frequency, normalized from 0 to 1, 
%1 corresponds to half the sample rate of the signal
[B2,A2] = butter(N2,Wn2,'low');
spike2Filtered= filter(B2,A2,spike2);
% spike2Filtered=spike2;
xlim([0 10]);

%downsample spike2 to run xcorr against logitech audio
endPoint=xoriginal(end);
newNumberPoints=floor(endPoint*soundRate);
newEndPoint=newNumberPoints/soundRate;
xnew=0:(1/soundRate):newEndPoint;
spike2Interp=interp1(xoriginal,spike2Filtered,xnew);

%look at original aligment 
figure(1);
plot(x,soundData); hold on;
plot(xoriginal,spike2/200,'r');
xlim([0 3]);
title('original logitech vs. spike2');

% %compare hpf logitech and lpf spike2
figure(2); plot(x,soundDataFiltered);
hold on;
plot(xoriginal,spike2/500,'r');
xlim([0 5]);
title('hpf logitech vs. lpf spike2');

%run a cross-correlation to get time lag of maximum similarity
r=xcorr(spike2Interp,soundDataFiltered);
%length of r will be the length of the longer of the 2 vectors
%times 2 and minus 1 ((length*2)-1);
xvector=1:length(r);
toSubtract=(length(r)+1)/2;
xvector=xvector-toSubtract;
xvector=xvector/soundRate;
figure(3); plot(xvector,r);
title('cross-correlation');
peakR=find(r==max(r));
timeOfPeak=xvector(peakR); %in seconds
display(timeOfPeak);

%if time of peak is negative, need to adjust time chunk grabbed from spike2
%file in the same direction (for example, time lag=-0.3, vid chunk grabbed
%from 3 to 6s, spike2 grabbed should be 2.7 to 5.7s)
%peak is negative if video started FIRST

%adjust spike2 times
h=figure(4);
offset=-timeOfPeak;
xAdjusted=xnew+offset;
plot(x,soundDataFiltered*20); %logitech audio
hold on; %spike2 audio
plot(xAdjusted,spike2Interp/10,'r'); %spike2 audio
title('adjusted data, with filtering and downsample');
saveas(h,'aligned spike2 and logitech audio.fig');

savename='spike 2 offset.mat';
save(savename,'xvector','r','peakR','timeOfPeak');