function StatMat=CalcEcgPeak
%   This software is protected by german copyright and international treaties.             
%   Copyright 2004 Markus Junghöfer & Peter Peyk. All Rights Reserved.                     
%                                                                                          
%   THIS SOFTWARE AND DOCUMENTATION IS PROVIDED "AS IS," AND COPYRIGHT HOLDERS MAKE        
%   NO REPRESENTATIONS OR WARRANTIES, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO,    
%   WARRANTIES OF MERCHANTABILITY OR FITNESS FOR ANY PARTICULAR PURPOSE OR THAT THE USE    
%   OF THE SOFTWARE OR DOCUMENTATION WILL NOT INFRINGE ANY THIRD PARTY PATENTS, COPYRIGHTS,
%   TRADEMARKS OR OTHER RIGHTS. COPYRIGHT HOLDERS WILL NOT BE LIABLE FOR ANY DIRECT,       
%   INDIRECT, SPECIAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF ANY USE OF THE SOFTWARE      
%   OR DOCUMENTATION.                                                                      
SampRate=687.17;
NPointsPerPic=round(333.333./(1000./SampRate));

% FilePath1='K:\MEG\FPP\Flucht\14\1.1\FPP14.1.1.ecg.pck';
% FilePath2='K:\MEG\FPP\Flucht\14\2.3\FPP14.2.3.ecg.pck';
% FilePath3='K:\MEG\FPP\Flucht\14\3.5\FPP14.3.5.ecg.pck';

FilePath1=[];
FilePath2=[];
FilePath3=[];

DefDataFilePath=SetDefPath('r','*.1.ecg.pck');
[File1,Path1,FilePath1]=ReadFilePath(FilePath1,DefDataFilePath,'Choose *.1.ecg.pck file');
if File1==0; return; end
SetDefPath('w',Path1);
a1=ReadAscii(FilePath1);

DefDataFilePath=SetDefPath('r','*.3.ecg.pck');
[File2,Path2,FilePath2]=ReadFilePath(FilePath2,DefDataFilePath,'Choose *.3.ecg.pck file');
if File2==0; return; end
SetDefPath('w',Path2);
a2=ReadAscii(FilePath2);

DefDataFilePath=SetDefPath('r','*.5.ecg.pck');
[File3,Path3,FilePath3]=ReadFilePath(FilePath3,DefDataFilePath,'Choose *.5.ecg.pck file');
if File3==0; return; end
SetDefPath('w',Path3);
a3=ReadAscii(FilePath3);

CondFilePath1=[FilePath1(1:length(FilePath1)-7),'CON'];
[CondVec1,CondFile1,CondPath1,CondFilePath1,EndPos1]=ReadData(CondFilePath1,1,[],'ascii','ascii');
[CondVec1,CondFile1,CondPath1,CondFilePath1,EndPos1]=ReadData(CondFilePath1,1,[],'ascii','ascii',EndPos1);
CondFilePath2=[FilePath2(1:length(FilePath2)-7),'CON'];
[CondVec2,CondFile2,CondPath2,CondFilePath2,EndPos2]=ReadData(CondFilePath2,1,[],'ascii','ascii');
[CondVec2,CondFile2,CondPath2,CondFilePath2,EndPos2]=ReadData(CondFilePath2,1,[],'ascii','ascii',EndPos2);
CondFilePath3=[FilePath3(1:length(FilePath3)-7),'CON'];
[CondVec3,CondFile3,CondPath3,CondFilePath3,EndPos3]=ReadData(CondFilePath3,1,[],'ascii','ascii');
[CondVec3,CondFile3,CondPath3,CondFilePath3,EndPos1]=ReadData(CondFilePath3,1,[],'ascii','ascii',EndPos3);

%=========================================================================================
hFig1=figure('Units','normal','Position',[0.075 0 .45  .95],'Name','ECG Peaks','NumberTitle','off');
hFig2=figure('Units','normal','Position',[.53 0 .47  .95],'Name','ECG Peaks','NumberTitle','off');

StdFak=2;
d1=diff(a1);
d2=diff(a2);
d3=diff(a3);
AbsMin=min([min(d1),min(d2),min(d3)]);
AbsMax=max([max(d1),max(d2),max(d3)]);
StatMat=zeros(3,9);
ColorStr='rbg';
TitleStr=char('Aversive vs. Neutral');
TitleStr=char(TitleStr,'Pleasant vs. Neutral');
TitleStr=char(TitleStr,'Aversive vs. Pleasant');
for Ind=1:3
    eval(['a0=a',int2str(Ind),';d0=d',int2str(Ind),';CondVec0=CondVec',int2str(Ind),';'])
    
    DiffCondVec0=diff(CondVec0);
    BreakInd=find(DiffCondVec0>15000)
    if isempty(BreakInd)
        error('No break?')
    end

    d0Std=std(d0);
    d0Mean=mean(d0);
    UpLim=d0Mean+StdFak.*d0Std;
    LoLim=d0Mean-StdFak.*d0Std;
    i0=find(d0<=UpLim & d0>=LoLim);
    j0=find(d0>UpLim | d0<LoLim);
    
    Nd0=length(d0);
    ni0=length(i0);
    nj0=length(j0);
    rn0=100.*(nj0./Nd0);
    
    
    Maxd0=max(d0);
    Mind0=min(d0);
    
    figure(hFig1);
    subplot(3,1,Ind)
    
    plot(a0(2:Nd0+1),d0,ColorStr(Ind))
    hold on
    plot([a0(2) a0(Nd0+1)],[UpLim UpLim],[ColorStr(Ind),'--'])
    plot([a0(2) a0(Nd0+1)],[LoLim LoLim],[ColorStr(Ind),'--'])
    plot([CondVec0(1) CondVec0(1)],[.9.*AbsMin 1.1.*AbsMax],'k')
    plot([CondVec0(BreakInd) CondVec0(BreakInd)],[.9.*AbsMin 1.1.*AbsMax],'k')
    plot([CondVec0(BreakInd+1) CondVec0(BreakInd+1)],[.9.*AbsMin 1.1.*AbsMax],'k')
    plot([CondVec0(end) CondVec0(end)],[.9.*AbsMin 1.1.*AbsMax],'k')
    if ~isempty(j0)
        plot(a0(j0+1),d0(j0),[ColorStr(Ind),'*'])
    end
    TrigVec1=[1:25:BreakInd];
    TrigVec2=[BreakInd+1:25:length(CondVec0)];
    for TrigInd=1:length(TrigVec1)
        plot([CondVec0(TrigVec1(TrigInd))+1 CondVec0(TrigVec1(TrigInd))+1],[.9.*AbsMin 1.1.*AbsMax],'k--')
    end
    for TrigInd=1:length(TrigVec2)
        plot([CondVec0(TrigVec2(TrigInd))+1 CondVec0(TrigVec2(TrigInd))+1],[.9.*AbsMin 1.1.*AbsMax],'k--')
    end
    axis([2 a0(Nd0+1) .9.*AbsMin 1.1.*AbsMax])
    title(TitleStr(Ind,:))
    xlabel('Time []')
    ylabel('R-wave distance [ms]')
    figure(hFig2)
    subplot(3,2,(Ind-1).*2+1)
    
    MinSeg1=CondVec0(1);
    MaxSeg1=CondVec0(BreakInd)+NPointsPerPic;
    xi=linspace(MinSeg1,MaxSeg1,1000);
    x1=a0(i0);
    y1=d0(i0);
    Seg1Ind=find(x1>=MinSeg1 & x1<=MaxSeg1);
    NInd1=x1(Seg1Ind(end))-x1(Seg1Ind(1))+1;
    
    Int1=interp1(x1(Seg1Ind),y1(Seg1Ind),xi,'linear','extrap');
    MeanInt1=mean(Int1);
    StdInt1=std(Int1);
    w1=linspace(MinSeg1,MaxSeg1,1000);
    plot(w1,Int1,ColorStr(Ind))
    hold on
    plot([w1(1) w1(end)],[MeanInt1 MeanInt1],['--',ColorStr(Ind)])
    axis([MinSeg1 MaxSeg1 AbsMin AbsMax])
    title(['Mean: ',int2str(MeanInt1)])
    xlabel('Time []')
    ylabel('R-wave distance [ms]')
    
    subplot(3,2,(Ind-1).*2+2)
    MinSeg2=CondVec0(BreakInd+1);
    MaxSeg2=CondVec0(end)+NPointsPerPic;
    xi=linspace(MinSeg2,MaxSeg2,1000);
    x2=a0(i0);
    
    y2=d0(i0);
    Seg2Ind=find(x2>=MinSeg2 & x2<=MaxSeg2);
    NInd2=x2(Seg2Ind(end))-x2(Seg2Ind(1))+1;
    
    Int2=interp1(x2(Seg2Ind),y2(Seg2Ind),xi,'linear','extrap');
    MeanInt2=mean(Int2);
    StdInt2=std(Int2);
    w2=linspace(MinSeg2,MaxSeg2,1000);
    plot(w2,Int2,ColorStr(Ind))
    hold on
    plot([w2(1) w2(end)],[MeanInt2 MeanInt2],['--',ColorStr(Ind)])
    axis([MinSeg2 MaxSeg2 AbsMin AbsMax])
    title(['Mean: ',int2str(MeanInt2)])
    xlabel('Time []')
    ylabel('R-wave distance [ms]')
    
    MeanInt12=mean([Int1,Int2]);
    StdInt12=std([Int1,Int2]);
    StatMat(Ind,1)=NInd1./1000;
    StatMat(Ind,2)=NInd2./1000;
    StatMat(Ind,3)=(NInd1+NInd2)./1000;
    StatMat(Ind,4)=MeanInt1;
    StatMat(Ind,5)=MeanInt2;
    StatMat(Ind,6)=MeanInt12;
    StatMat(Ind,7)=StdInt1;
    StatMat(Ind,8)=StdInt2;
    StatMat(Ind,9)=StdInt12;   
    StatMat(Ind,10)=ni0;
    StatMat(Ind,11)=nj0;
    StatMat(Ind,12)=rn0;
    
end
SaveAscii(StatMat,[FilePath1(1:length(FilePath1)-8),'.msna'])
%=========================================================================================
%=========================================================================================
