function [eegTimes ] = PsychLSLAudioAnalysis(allDataStreams)
%Analysis function test test synchrony between LSL events and audio data.
%Written by Matthew Grivich, SCCN/UCSD.


    %These parameters will need bo be changed to match the data.
    eegStreamNumber = 2;
    eegChannelNumber = 1;
    audioStreamNumber = 1;
    audioChannelNumber = 1;
    
    snapMarkerStreamNumber = 3;
    stdThreshold = 1.5;
    eegThreshold = 1.5;
    desiredEvent = 'tone';
    %desiredEvent = 'tone-after';
    
    preTime = .010;
    
    eeg =  allDataStreams.item{eegStreamNumber}.data(:,eegChannelNumber); 
 
    eeg = eeg-mean(eeg);
    eegTime = allDataStreams.item{eegStreamNumber}.timeStamp;
    
    
    eventTimes = allDataStreams.item{snapMarkerStreamNumber}.timeStamp;
    eventVals = allDataStreams.item{snapMarkerStreamNumber}.event.label; 

    audio =  allDataStreams.item{audioStreamNumber}.data(:,audioChannelNumber);
    audioTime = allDataStreams.item{audioStreamNumber}.timeStamp;
    
    %eegTime = eegTime - audioTime(450);
    %eventTimes = eventTimes - audioTime(450);
    %audio = audio(500:end);
    %audioTime = -audioTime(450) + audioTime(500:end);

    audioStd = std(audio);
    eegStd = std(eeg);
    
    figure;
    hold on;
    plot(audioTime(1:10:end), audio(1:10:end)/audioStd, 'color', 'black');
    plot(eegTime, eeg/eegStd+10, 'color', 'blue');
    leg = legend('Audio Capture', 'EEG');
    set(leg,'Location', 'NorthEast');
    set(leg,'Interpreter', 'None');
    
    xlabel('Time (s)');
    ylabel('Amplitude / Standard Deviation');
 

    audioIndex = 1;
    eegIndex = 1;
    
    eegLatencies = zeros(1,length(eventTimes));
    audioLatencies = zeros(1,length(eventTimes));
    eegTimes = zeros(1,length(eventTimes));
    eventCount = 0;
    %for each event, find eeg and audio latency (difference between event
    %marker and onset of audio.
    for i=1:length(eventTimes)-10
    %    if((eventTimes(i) > 655 && eventTimes(i) < 656) || (eventTimes(i) > 1529 && eventTimes(i) < 1531)) 
    %        continue;
    %    end
        if strcmp(eventVals(i,:),desiredEvent) == 1 || strcmp(eventVals(i),desiredEvent) == 1
            
         
            eventCount = eventCount + 1;

            line([eventTimes(i),eventTimes(i)], get(gca,'Ylim'), 'color','r')
            while (eegTime(eegIndex) < eventTimes(i) - preTime || eeg(eegIndex)/eegStd < eegThreshold) && eegIndex <= length(eeg)
                eegIndex = eegIndex + 1;
            end
            
            %interpolate
            eegLatencies(eventCount) = eegTime(eegIndex-1) + (eegTime(eegIndex) - eegTime(eegIndex-1)) * (eegThreshold*eegStd - eeg(eegIndex-1)) / (eeg(eegIndex)-eeg(eegIndex-1)) - eventTimes(i);
            
            eegTimes(eventCount) = eegTime(eegIndex);
            
            while (audioTime(audioIndex) < eventTimes(i) -preTime || audio(audioIndex)/audioStd < stdThreshold) && audioIndex <= length(audio)
                audioIndex = audioIndex + 1;
            end

            %interpolate
            audioLatencies(eventCount) = audioTime(audioIndex-1) + (audioTime(audioIndex) - audioTime(audioIndex-1)) * (stdThreshold*audioStd - audio(audioIndex-1)) / (audio(audioIndex)-audio(audioIndex-1)) - eventTimes(i);
           
            
        end
        
        
    end

    eegLatencies = eegLatencies(1:eventCount);
    eegTimes = eegTimes(1:eventCount);
    audioLatencies = audioLatencies(1:eventCount);
    
    line(get(gca,'Xlim'),[eegThreshold + 10,eegThreshold + 10],'color','red');
    line(get(gca,'Xlim'),[stdThreshold,stdThreshold],'color','red');
    hold off
    
    %Make various analysis plots
    figure
    hold all
    scatter(eegLatencies, audioLatencies);
    xlabel('EEG Latency (s)');
    ylabel('Audio Capture Latency (s)');

    hold off
    
    histFig( eegLatencies*1000, 'EEG Latency (ms)' );
   
    histFig( audioLatencies*1000, 'Audio Capture Latency (ms)' );   

    figure
    plot(eegTimes,eegLatencies - audioLatencies);
    ylabel('EEG delay - audio delay (s)');
    xlabel('EEG time (s)');
    
    figure
    plot(eegTimes,eegLatencies);
    ylabel('EEG delay (s)');
    xlabel('EEG time (s)');
    
    figure
    plot(eegTimes,audioLatencies);
    ylabel('audio delay (s)');
    xlabel('EEG time (s)');
    
end

