%   CalcEcgRaw.m

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


%=========================================================================================
DefDataFilePath=SetDefPath('r','*.raw');
[File,Path,FilePath]=ReadFilePath([],DefDataFilePath);
if File==0; return; end
SetDefPath('w',Path);
                                                                    
a0=ReadData(FilePath,1,[],'float32','float32');
a1=a0;
%=========================================================================================
SampRate=678.17;
Nyquist=SampRate./2;
LineFrequ=50;
SampRate=1;
HighFiltStatus=1;
HighFrequ=2;
HighN=6;
LowFiltStatus=1;
LowFrequ=100;
LowN=6;
StopFiltStatus=1;
StopFrequ=[LineFrequ-1 LineFrequ+1];
StopN=6;
if HighFiltStatus
	fprintf(1,'Highpass filtering of ECG data: HighFrequ: %g Hz; Order: %g;\n\n',HighFrequ,HighN);
	[BHigh,AHigh]=butter(HighN,HighFrequ./Nyquist,'high');
	a1=filtfilt(BHigh,AHigh,a1);
else
	fprintf(1,'No highpass filtering of ECG data.\n\n');
end
if LowFiltStatus
	fprintf(1,'Lowpass filtering of ECG data: LowFrequ: %g Hz; Order: %g;\n\n',LowFrequ,LowN);
	[BLow,ALow]=butter(LowN,LowFrequ./Nyquist,'low');
	a1=filtfilt(BLow,ALow,a1);
else
	fprintf(1,'No lowpass filtering of ECG data.\n\n');
end
if StopFiltStatus
	fprintf(1,'Stopband filtering of ECG data: LineFrequ: %g Hz; Order: %g;\n\n',LineFrequ,StopN);
	[BStop,AStop]=butter(StopN,StopFrequ./Nyquist,'stop');
	a1=filtfilt(BStop,AStop,a1);
else
	fprintf(1,'No stopband filtering of ECG line noise.\n\n');
end
%=========================================================================================
NPoints=length(a1);
NPoints2Plot=2000;
%=========Quick test for data reversal=========
NInt=10;
NPointsVec=floor(linspace(1,NPoints,NInt));
ReverseStatus=0;
for i=1:NInt-1
	MaxVec(i)=max(a1(NPointsVec(i):NPointsVec(i+1)));
	MinVec(i)=min(a1(NPointsVec(i):NPointsVec(i+1)));
	MeanVec(i)=mean(a1(NPointsVec(i):NPointsVec(i+1)));
end
MedMax=median(MaxVec);
MedMin=median(MinVec);
MedMean=median(MeanVec);
if abs(MedMax-MedMean)<abs(MedMin-MedMean);
	ReverseStatus=1;
end
if ReverseStatus
	RawEcgFigTitle='Reversed raw ECG';
	FiltEcgFigTitle='Filtered and reversed raw ECG';
	a1=-a1;
else
	RawEcgFigTitle='Raw ECG';
	FiltEcgFigTitle='Filtered raw ECG';
end
%=========================================================================================
hRawEcg=figure('Units','normal','Position',[0 .6 1 .3],'Name',RawEcgFigTitle,'Color',[1 1 1],'NumberTitle','off');
plot(a1)
axis([1, NPoints,min(a1),max(a1)])

c1=a1;
hPeak=figure('Units','normal','Position',[0 .3 1 .3],'Name',FiltEcgFigTitle,'Color',[1 1 1],'NumberTitle','off');

while 1
	plot(c1)
	axis([1,NPoints,0,1.1.*max(c1)])
	hold on
	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:NPoints]';
		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	
	title('Left mouse button => Redo; Right mouse button => Accept');  
	[x3,y3,z3]=ginput(1);
	if z3>1
		break;
	end	
end
%====================================================================================
figure(hRawEcg)
plot(c1)
axis([1,NPoints,0,1.1.*max(c1)])
plot(a1)
hold on
plot(n1,a1(n1),'r*')
hold off
axis([1,NPoints,1.1.*min(a1),1.1.*max(a1)])
%====================================================================================
SaveAscii(n1,[FilePath(1:length(FilePath)-4),'.pck'])
