
function TESTPsychLSLScreenTimingAnalysis(allDataStreams)
%Analysis function to test synchrony between LSL events and screen flash
%events.
%Written by Matthew Grivich SCCN/UCSD.

    %These parameters will need to be changed to match the data.
    eegStreamNumber = 1;
    eegChannelNumber = 1;
    stdThreshold = 0;
 
    engineStreamNumber = 2;
    markerStreamNumber = 2;
    engineEvent = 'flip';
    onEvent = 'white';
    offEvent = 'black';
    
    preTime = .10;
 
    
    
    
    eeg =  allDataStreams.item{eegStreamNumber}.data(:,eegChannelNumber); 
    eegTime = allDataStreams.item{eegStreamNumber}.timeStamp;
    eventLatencies = allDataStreams.item{markerStreamNumber}.timeStamp; 
    eventVals = char(allDataStreams.item{markerStreamNumber}.event.label);
    eventLatencies = eventLatencies(1:end-1);


    engineLatencies = allDataStreams.item{engineStreamNumber}.timeStamp;
    engineVals = char(allDataStreams.item{engineStreamNumber}.event.label);

    
    samplingRate = allDataStreams.item{eegStreamNumber}.samplingRate;
    
    %filter out high frequency oscillations common in monitor light levels
    [z,p,k] =  butter(9,40/(samplingRate/2),'low'); % 9th order,  40 hz / (samplingrate/2), low pass

    [sos,g] = zp2sos(z,p,k);
    eeg = filtfilt(sos,g,eeg);
 
    onTimes = zeros(1, length(eventLatencies));
    offTimes = zeros(1, length(eventLatencies));
    onCount = 0;
    offCount = 0;
     

    for i=1:length(eventLatencies)-1
         if strcmp(strtrim(eventVals(i,:)),onEvent)
            onCount = onCount+1;
            onTimes(onCount) = eventLatencies(i);

         end       
       

         if strcmp(strtrim(eventVals(i,:)),offEvent)
            offCount = offCount+1;
            offTimes(offCount) = eventLatencies(i);

         end
    end    
    
    onTimes = onTimes(1:onCount);
    offTimes = offTimes(1:offCount);
    eegStd = std(eeg);
    eeg = (eeg-mean(eeg))/eegStd;
    figure;
    hold all;
    plot(eegTime, eeg, 'color', 'blue');
    leg = legend('EEG');
    set(leg,'Location', 'NorthEast');
    set(leg,'Interpreter', 'None');
    
    xlabel('Time (s)');
    ylabel('Volts');
    ylabel('Amplitude / Standard Deviation');
    
    eegIndex = 1;
   engineCount = 0;
   ylim = get(gca,'Ylim');
    for i=1:1000 %length(engineLatencies)
        if strcmp(strtrim(engineVals(i,:)),engineEvent)
        %if strcmp(engineVals(i,:),engineEvent)
            engineCount = engineCount+1;
            line([engineLatencies(i),engineLatencies(i)], ylim, 'color','g')
        end
    end    

    onCount = 0;
    offCount = 0;
    onLatencies = zeros(1,length(onTimes));
    offLatencies = zeros(1,length(offTimes));
    
    %for each event, find screen latency (difference between event marker
    %and onset of light transition)
    for i=1:length(onTimes) 
       
            line([onTimes(i),onTimes(i)], ylim, 'color','r')
 
            while eegIndex < length(eeg) && (eegTime(eegIndex) < onTimes(i) - preTime || eeg(eegIndex) < stdThreshold) 
               
                eegIndex = eegIndex + 1;
            end
    %        onLatencies(i) = eegTime(eegIndex) - onTimes(i);
            %interpolate 
            onLatencies(i) = eegTime(eegIndex-1) + (eegTime(eegIndex) - eegTime(eegIndex-1)) * (stdThreshold*eegStd - eeg(eegIndex-1)) / (eeg(eegIndex)-eeg(eegIndex-1)) - onTimes(i);
   
             onCount = onCount+1;
        
    end
    eegIndex = 1;
     for i=1:length(offTimes) 
       
            line([offTimes(i),offTimes(i)], ylim, 'color','r')

            while eegIndex < length(eeg) && (eegTime(eegIndex) < offTimes(i) - preTime || eeg(eegIndex) > stdThreshold) 
                eegIndex = eegIndex + 1;
            end
            %offLatencies(i) = eegTime(eegIndex) - offTimes(i);
            %interpolate
            offLatencies(i) = eegTime(eegIndex-1) + (eegTime(eegIndex) - eegTime(eegIndex-1)) * (stdThreshold*eegStd - eeg(eegIndex-1)) / (eeg(eegIndex)-eeg(eegIndex-1)) - offTimes(i);
   
            offCount = offCount+1;
    end
    
    line(get(gca,'Xlim'),[stdThreshold,stdThreshold],'color','red');
    hold off
   

    histFig( onLatencies*1000, 'On Latency (ms)' );

    histFig(offLatencies*1000, 'Off Latency (ms)');
    
    figure
    plot(onTimes,1000*(onLatencies));
    ylabel('On delay (ms)');
    xlabel('On time (s)');
    
    figure
    plot(offTimes,1000*(offLatencies));
    ylabel('Off delay (ms)');
    xlabel('On time (s)');
    
end

