cd D:\blink\blink_magic\eyedata
 clear all
 load gaze_all
 
 for i=1:20 
     i
     dat=gaze{i}(:,4);
     
     dat(isnan(dat))=0;
     fs=1000;
     windowSize=20;%20ms smoothing
     pupil=filtfilt(ones(1,windowSize)/windowSize,1,dat);
     plot(pupil)
    
     % detect close time
    thre=input('threshold');
    close all
    
    cl=zeros(length(pupil),1);
    cl(find(pupil<thre))=1;
    
    % detect on and offset
    on=find(diff(cl)==1);
    off=find(diff(cl)==-1);

    if length(on) > length(off)
        on(length(on))=[];
    else
    end
    
    if on(1)>off(1)
        off(1)=[];
    else
    end
    dur=off-on;

    thre1=1000; % maximum blink duration
    thre2=20;%minimum blink duration
    
    del=[find(dur<thre2)';find(dur>thre1)'];
    bk_onset=on-20; % 10ms ahead
    bk_offset=off+20;bk_dur=dur;
    
    bk_onset(del)=[];
    bk_offset(del)=[];
    bk_dur(del)=[];
    
    df=diff(bk_onset);
    del2=find(df<100);
    bk_onset(del2)=[];

    figure
    plot(pupil);hold on
    plot(bk_onset, pupil(bk_onset),'ro');
    pause
    close all
    
    bk_all{i}=bk_onset;
    clear pupil bk_onset dat

 end
 
 save gaze_all bk_all -append

% bk duration :300ms on(1), off(0)
bk_hist=zeros(124500,20);
for i=1:20
   for j=1:length(bk_all{i})
       bk_hist(bk_all{i}(j):bk_all{i}(j)+300,i)=1;
   end
end

all_bk_hist=sum(bk_hist,2);
sr_bk_hist=zeros(124500,1000);
for k=1:1000
    dat=zeros(124500,20);
    for i=1:20
        %shuffling bk onset time
        df1=[bk_all{i}(1);diff(bk_all{i});124001-bk_all{i}(length(bk_all{i}))];
        rdf=df1(randperm(length(df1)));
        sr_bk=cumsum(rdf);
        sr_bk(length(sr_bk))=[];
       for j=1:length(sr_bk)
           dat(sr_bk(j):sr_bk(j)+300,i)=1;
       end
    end
    sr_bk_hist(:,k)=sum(dat,2);
    %sr_histc(:,k)=histc(sr_bk_hist(1:124000,k),bin);
end

% z-score based on the 1000 surrogate data
th_av=mean(sr_bk_hist,2);
th_sd=std(sr_bk_hist');

th_zscore=zeros(124500,1);
th_p=zeros(124500,1);
for i=1:124500
        [h,p,ci,zval] =ztest(all_bk_hist(i),th_av(i),th_sd(i));
        th_zscore(i,1)=zval;
        th_p(i,1)=p;
end

% z -value transformation
z_bk=(all_bk_hist-mean(sr_bk_hist,2))./std(sr_bk_hist')';
method=[20.7,20.73;30.9,31.0; 36.9, 37.1; 38.86, 39.06; 51.60, 52.00; 59.83, 60.00; 61.86, 62.00];
effect=[24.33,26.66; 33.23, 35.60; 40.90, 49.63; 55.26, 59.90; 63.00, 77.16; 117.63, 121.10];
st_final=[28.04;31.55;36.28;51.36;61.21;77.94;93.07;104.14;107.22;112.41;123.61];
figure(2)
plot(1/1000:1/1000:124.5,all_bk_hist,'linewidth',2);hold on;box off;axis([0 124 0 10]);xlabel('time(sec)');ylabel('blink frequency');
hold on
plot(st_final, all_bk_hist(st_final*1000),'ro','linewidth',2); % more than 6
for i=1:length(effect)
line(effect(i,:),[5.5 5.5],'linewidth',4)
end

for i=1:length(method)
line(method(i,:),[7 7],'linewidth',4)
end


% correlogram for each pair
 clear all
 cd D:\blink\blink_magic\eyedata
 load gaze_all

 clear z_hist real_hist
bin=-1950:300:1950;
k=1;
for i=[1:20] 
    for j=[1:20]
        if i==j
        else
        dat1=bk_all{i};
       % each pair　ごとに、1000回サロゲートを作って、z-valueで評価する
        dat2=bk_all{j};
        clear dfp
        for s=1:length(dat1)
            dfp(:,s)=dat2-dat1(s);
        end
        rh=sum(histc(dfp,bin),2)/length(dat1);
        real_hist(:,k)=rh;
        
        % surrogate data of dat2
        for t=1:1000
            clear sr_bk df1 rdf sdfp 
        df1=[dat2(1);diff(dat2);124001-dat2(length(dat2))];
        rdf=df1(randperm(length(df1)));
        sr_bk=cumsum(rdf);
        sr_bk(length(sr_bk))=[];
        
             for s=1:length(dat1)
                sdfp(:,s)=sr_bk-dat1(s);
            end
        srh(:,t)=sum(histc(sdfp,bin),2)/length(dat1);
        end
        srh_m=mean(srh,2);
        srh_sd=std(srh')';
        z_hist(:,k)=(rh-srh_m)./srh_sd;
        k=k+1;
         end
       
    end
end

z_hist(size(z_hist,1),:)=[];
tb=-1800:300:1800;
mz=mean(z_hist,2);
sez=std(z_hist')'/sqrt(190);
figure(1)
errorbar(tb,mz,sez,'-ro','linewidth',3);
set(gca,'XTick',[-1800:300:1800],'XtickLabel',[-1.8:0.3:1.8]);
xlabel('Time(s)');ylabel('z-value for frequency of blink onset asynchrony');
%plot(tb,mean(z_hist,2),'-ro');hold on

% 
% for i=1:size(z_hist,1)
%     [h(i) p(i)]=ttest(z_hist(i,:),0);
% end

save gaze_all h p z_hist real_hist tb mz sez -append