function [f,TrialIndVec] = NeuroScan2Ecg(FileMat)


%   NeuroScan2Ecg

%   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.                                                                      

%global a1 b1 PeakDiffMat
global PeakDiffMat
a1=[]
b1=[]
PeakDiffMat=[]

MeanPeakMat=[]
StdPeakMat=[]
MedPeakMat=[]
%==================================================================
%
%	NeuroScan2Ecg
%
%
%	Markus Junghoefer	[2002]
%
%
%	Function definition
%
%	function NsStruct=NeuroScan2Ecg(NsFilePath)
% 
% 	if nargin<1; FilePath=[]; end
%

%=================================================================
if nargin<1;FileMat=[]; end

%FileMat='Rudi:ApplicationsOS9:Matlab5.2:Plot2d3d:Plot2dTestData:EEG:BinData:Neuroscan:06fastpic2.1.cnt'
DefFilePath=SetDefPath('r','*.cnt');
[NFiles,FileMat,VersionString,BatchFilePath]=ReadFileNames(FileMat,DefFilePath,'Choose cnt file(s) or a batch file:');
if NFiles==0; return; end
[BatchFile,BatchPath]=SepFilePath(BatchFilePath)
SetDefPath('w',BatchPath);

for FileInd=1:NFiles
	[NsFile,NsPath,NsFilePath]=GetFileNameOfMat(FileMat,FileInd);
	
	f=ReadNeuroScanStruct(NsFilePath);
	NsHeader=f.header;
	NsEvents=f.events;
	%==========================================
	NTotTrials=NsHeader.numevents;
	NNsChan=NsHeader.nchannels;
	SampRate=NsHeader.rate
	
	EcgChan=64;
	CellVec=[1 2 3 4 5];
	NBadStimtype=0;
	NOKStimtype=0;
	TrialStatusVec=zeros(NTotTrials,1);
	for i=1:NTotTrials
		StimtypeTmp=NsEvents(i).stimtype;
		if StimtypeTmp>0 & StimtypeTmp<6
			TrialStatusVec(i)=1;
		end
		StimTypeVec(i)=NsEvents(i).stimtype;
		StimOffsetVec(i)=NsEvents(i).offset;	
	end
	UseTrialIndVec=find(TrialStatusVec);
	NonUseTrialIndVec=find(~TrialStatusVec)
	NTotTrials=length(UseTrialIndVec);
	NCells=length(CellVec);
	for i=1:NCells
		NTrialsVec(i)=length(find(StimTypeVec==CellVec(i)));
	end	
	NTrialsVec=NTrialsVec
	StimTypeVec=StimTypeVec(UseTrialIndVec);
	StimOffsetVec=StimOffsetVec(UseTrialIndVec);
	StimTypeMat=zeros(NCells,10);
	StimOffsetMat=zeros(NCells,10);
	for CellInd=1:NCells
		Tmp1=find(StimTypeVec==CellInd);
		DiffTmp1=diff(Tmp1);
		Tmp2=find(DiffTmp1>1);
		StimTypeMat(CellInd,:)=Tmp1([1,Tmp2+1]);
		StimOffsetMat(CellInd,:)=StimOffsetVec(Tmp1([1,Tmp2+1]));
	end	
	StimTypeMat
	StimOffsetMat
	%=========Ecg Header======================
	EcgFilePath=[NsFilePath(1:length(NsFilePath)-4),'.ecg'];
	%EcgFid=fopen(EcgFilePath,'w','b');
	
	fprintf(1,'\n\n\n');
	fprintf(1,'Total cells:          %g\n',NCells);
	fprintf(1,'Total trials:         %g\n',NTotTrials);
	%==========================================
	Nyquist=SampRate./2;
	LineFrequ=50;
	StopFrequ=[LineFrequ-1 LineFrequ+1];
	LowFrequPass=40;
    LowFrequStop=50;
	HighFrequPass=0.5;
    HighFrequStop=0.1;
    PassdB=3;
    StopdB=45;
	LowN=6;
	HighN=3;
	StopN=6;
    [NHigh,WnHigh] = buttord(HighFrequPass*2/SampRate,HighFrequStop*2/SampRate,PassdB,StopdB)
    [BHigh,AHigh]  = butter(NHigh,WnHigh,'high')
    [NLow,WnLow]   = buttord(LowFrequPass*2/SampRate,LowFrequStop*2/SampRate,PassdB,StopdB)
	[BLow,ALow]    = butter(NLow,WnLow,'low');
	[BStop,AStop]  = butter(StopN,StopFrequ./Nyquist,'stop');
	NNsPoints=0;
	Pre=0
	Post=3150
	PreSeek=Pre.*NNsChan.*2;
	PostSeek=Post.*NNsChan.*2;
	NPoints=(StimOffsetVec(end)-StimOffsetVec(1))./(NNsChan.*2)+Pre+Post+1;
	StimPointMat=(StimOffsetMat-StimOffsetVec(1))./(NNsChan.*2)+Pre+1;
	if isempty(a1)
		a1=zeros(NPoints,1);
		NsFid=fopen(NsFilePath,'r','l')
		FSeekStatus=fseek(NsFid,StimOffsetVec(1)-PreSeek,'bof');
		if FSeekStatus==-1
			error('fseek error in NeuroScan2Ecg')
		end
		for PoindInd=1:NPoints
			Tmp=fread(NsFid,NNsChan,'int16');
			a1(PoindInd)=Tmp(EcgChan,:);
		end
		fclose(NsFid);
	end

	%=========Filter data =========
	if isempty(b1)
		b1=filtfilt(BLow,ALow,a1);
    	b1=filtfilt(BHigh,AHigh,b1);
	end
	%=========================================================================================
	hRawEcg=figure('Units','normal','Position',[0 .65 1 .3],'Name','Raw ECG','NumberTitle','off','Color',[1 1 1]);
	FiltEcgFigTitle='Filtered raw ECG';
	hPeak=figure('Units','normal','Position',[0 .33 1 .3],'Name',FiltEcgFigTitle,'NumberTitle','off','Color',[1 1 1]);
	hBlockEcg=figure('Units','normal','Position',[0 0 1 .3],'Name','Block Ecg','NumberTitle','off','Color',[1 1 1]);
	hDistPeak=figure('Units','normal','Position',[.4 0 .3 .4],'Name','R-distance','NumberTitle','off','Color',[1 1 1]);
    hMedDistPeak=figure('Units','normal','Position',[.7 0 .3 .4],'Name','Median R-distance','NumberTitle','off','Color',[1 1 1]);
	Minb1=min(b1);
	Maxb1=max(b1);
	figure(hRawEcg)
	plot(b1)
	hold on
	ColVec='kbcrm'
	for CellInd=1:NCells
		for i=1:10
			plot([StimPointMat(CellInd,i) StimPointMat(CellInd,i)],[Minb1 Maxb1],ColVec(CellInd))
		end
	end
	axis([1 1000 min(b1(1:1000)) max(b1(1:1000))])
    
    %=========Check for data reversal=========
    Reversed=questdlg('Data reversed ?','Data reversed ?','No','Yes','No');
    switch Reversed,
        case 'Yes', 
            set(hRawEcg,'Name','Reversed Raw ECG');
            set(hPeak,'Name','Filtered and Reversed Raw ECG');
            a1=-a1;
            b1=-b1;
    end % switch
    %=========================================
    set(hRawEcg,'Name','Please choose an ECG pattern (interval)')

	[x1,y1,z1]=ginput(2);
	Pat=b1(x1(1):x1(2));
	NPat=length(Pat);
    axis([x1(1) x1(2) min(Pat) max(Pat)])
    pause(2)
    set(hRawEcg,'Name','Raw Ecg')
	axis([1 NPoints Minb1 Maxb1])
    PreBlock=-125;
	BlockLength=3150;
	PostBlock=0;
    NPeakDiffOld=0;
    NBlock=10;
    PeakDiffMat=zeros(2,NCells,NBlock,30);
	for CellInd=1:NCells
		for BlockInd=1:NBlock
			
			figure(hRawEcg)
			plot(StimPointMat(CellInd,BlockInd),Maxb1,['*',ColVec(CellInd)])
			
			c1=b1(StimPointMat(CellInd,BlockInd)-PreBlock:StimPointMat(CellInd,BlockInd)+BlockLength+PostBlock);
			Nc1=length(c1);
			Minc1=1.1.*min(c1);
			Maxc1=1.1.*max(c1);
			
			figure(hBlockEcg)	
			plot(c1,ColVec(CellInd))
			hold on
			plot([PreBlock+1 PreBlock+1],[Minc1 Maxc1],'k')
			plot([PreBlock+1+BlockLength PreBlock+1+BlockLength],[Minc1 Maxc1],'k')
			axis([1 Nc1 Minc1 Maxc1])
			hold off
			
            c1=xcorr(c1,Pat);
            v=version; 
            if str2num(v(1))>5
			    c1=c1(Nc1:end);
            else
                c1=c1(1:Nc1+NPat);
            end

			figure(hPeak)
			Nc1=length(c1);
			Minc1=1.1.*min(c1);
			Maxc1=1.1.*max(c1);
			
			while 1
				plot(c1,ColVec(CellInd))
				hold on
				plot([PreBlock+1 PreBlock+1],[Minc1 Maxc1],'k')
				plot([PreBlock+1+BlockLength PreBlock+1+BlockLength],[Minc1 Maxc1],'k')
				axis([1 Nc1 Minc1 Maxc1])
				title('Left mouse button => Peak above ginput of point; Right mouse button => Peak above ginput of line');  
				n1=[];
				[x1,y1,z1]=ginput(1);
				if z1==1
					m1=find(c1>y1);
					n1(1)=m1(1);
					j=1;
					for i=2:length(m1)
						if m1(i)-m1(i-1)==1
							if c1(m1(i))>c1(n1(j));
								n1(j)=m1(i);					
							end
						else
							j=j+1;
							n1(j)=m1(i);
						end
					end
					plot(n1,c1(n1),'r*')
					hold off
				else
					[x2,y2,z2]=ginput(1);
					m=(y2-y1)./(x2-x1);
					mx0=y1-m.*x1;
					m1=mx0+m.*[1:length(c1)]';
					plot([1:NPoints],m1,'k--')
					m1=find(c1-m1>0);
					n1(1)=m1(1);
					j=1;
					for i=2:length(m1)
						if m1(i)-m1(i-1)==1
							if c1(m1(i))>c1(n1(j));
								n1(j)=m1(i);					
							end
						else
							j=j+1;
							n1(j)=m1(i);
						end
					end
					plot(n1,c1(n1),'r*')
					hold off
				end	
                PeakTime=n1./SampRate.*1000
                NPeakTime=length(PeakTime);
				PeakDiff=diff(n1)./SampRate.*1000
                NPeakDiff=length(PeakDiff);
                if str2num(v(1))<6
                    PeakDiff=fliplr(PeakDiff);
                end
                PeakDiffMat(1,CellInd,BlockInd,1:NPeakTime)=PeakTime;
                PeakDiffMat(2,CellInd,BlockInd,1:NPeakDiff)=PeakDiff;
                
				MeanPeak=mean(PeakDiff)
				StdPeak=std(PeakDiff)
				MedPeak=median(PeakDiff)
				MeanPeakMat(CellInd,BlockInd)=MeanPeak;
				StdPeakMat(CellInd,BlockInd)=StdPeak;
				MedPeakMat(CellInd,BlockInd)=MedPeak;
				BlockMeanPeak=mean(MeanPeakMat(1:CellInd,1:BlockInd),2)
				BlockStdPeak=std(StdPeakMat(1:CellInd,1:BlockInd)')
				BlockMedPeak=median(MedPeakMat(1:CellInd,1:BlockInd)')
				title('Left mouse button => Redo; Right mouse button => Accept'); 
				set(hPeak,'Name',[FiltEcgFigTitle,'; Cell: ',int2str(CellInd),'; Block: ',int2str(BlockInd),'; NPeak: ',int2str(length(PeakDiff)),'; Mean: ',int2str(MeanPeak),'; Block Mean: ',int2str(BlockMeanPeak(CellInd))])
                pause(.01)
                %break;
 				[x3,y3,z3]=ginput(1);
 				if z3>1
 					break;
                end	
			end
			figure(hDistPeak)
			plot(PeakDiff,ColVec(CellInd))
			hold on
            figure(hMedDistPeak)
            PeakDiffMatTmp=squeeze(PeakDiffMat(2,:,:,:));
			MedDiffTmp=squeeze(median(PeakDiffMatTmp,2));
            SizeMedDiffTmp=size(MedDiffTmp);
            b=median(MedDiffTmp(:,2:min(13,SizeMedDiffTmp(2))),2);
            for CellIndTmp=1:NCells
                plot(MedDiffTmp(CellIndTmp,2:min(13,SizeMedDiffTmp(2))),ColVec(CellIndTmp))
                hold on
            end
            pause(.1)
		end
	end
end
% SaveNDimData(MeanPeakMat,[NsFilePath(1:length(NsFilePath)-4),'.mea'])
% SaveNDimData(StdPeakMat,[NsFilePath(1:length(NsFilePath)-4),'.std'])
% SaveNDimData(MedPeakMat,[NsFilePath(1:length(NsFilePath)-4),'.med'])
SaveNDimData(PeakDiffMat,[NsFilePath(1:length(NsFilePath)-4),'.pmat'])

%PeakDiffMat=ReadNDimData([],[],'*.pmat')
%SaveNDimData(PeakDiffMat,[FileMat(1:length(FileMat)-4),'.pmat'])

Tmp=mean(StdPeakMat,2)
Tmp=mean(StdPeakMat(:,2:10),2)

PeakDiffMatTmp=squeeze(PeakDiffMat(2,:,:,:));
a=squeeze(median(PeakDiffMatTmp,2))
b=median(a(:,2:13),2)
hMedMedDistPeak=figure('Units','normal','Position',[.1 0 .3 .4],'Name','Median','NumberTitle','off','Color',[1 1 1]);
for CellInd=1:NCells
    plot(a(CellInd,2:13),ColVec(CellInd))
    hold on
end
close(hRawEcg)
close(hPeak)
close(hBlockEcg)
pause(2)
close(hDistPeak)
close(hMedDistPeak)
close(hMedMedDistPeak)
return
