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

% a1=ReadAscii('Rudi:Data:MEG:FPP:Flucht:45:1.3:FPP45.1.3.ecg');
% a2=ReadAscii('Rudi:Data:MEG:FPP:Flucht:45:2.1:FPP45.2.1.ecg');
% a3=ReadAscii('Rudi:Data:MEG:FPP:Flucht:45:3.5:FPP45.3.5.ecg');
% a1=ReadAscii('Rudi:Data:MEG:FPP:Flucht:30:1.3:FPP30.1.3.ecg');
% a2=ReadAscii('Rudi:Data:MEG:FPP:Flucht:30:2.1:FPP30.2.1.ecg');
                                                                     
% a3=ReadAscii('Rudi:Data:MEG:FPP:Flucht:30:3.5:FPP30.3.5.ecg');
% a1=ReadAscii('Rudi:Data:MEG:FPP:Flucht:44:2.3:FPP44.2.3.ecg');
% a2=ReadAscii('Rudi:Data:MEG:FPP:Flucht:44:1.1:FPP44.1.1.ecg');
% a3=ReadAscii('Rudi:Data:MEG:FPP:Flucht:44:3.5:FPP44.3.5.ecg');
% a1=ReadAscii('Data:Data:MEG:FP:FPP:Flucht:30:1.3:FPP30.1.3.ecg');
% a2=ReadAscii('Data:Data:MEG:FP:FPP:Flucht:30:2.1:FPP30.2.1.ecg');
% a3=ReadAscii('Data:Data:MEG:FP:FPP:Flucht:30:3.5:FPP30.3.5.ecg');
a1=ReadAscii([],'Choose *.3.ecg file',[],[],'*.3.ecg.raw.p');
a2=ReadAscii([],'Choose *.1.ecg file',[],[],'*.1.ecg.raw.p');
a3=ReadAscii([],'Choose *.5.ecg file',[],[],'*.5.ecg.raw.p');
%=========================================================================================
% SampRate=1;
% ch_LowPassFilt=0;
% ch_HighPassFilt=1;
% HighDefFreqsMat=[SampRate 0.01 0.02];
% HighDefRipMat=[3 45];
% LowDefFreqsMat=[SampRate 10 20]; 
% LowDefRipMat=[3 45];
% 
% [ch_LowPassFilt,LowB,LowA,ch_HighPassFilt,HighB,HighA,LowFreqsMat,HighFreqsMat] = ...
% GetHighLowFiltCoeff(SampRate,HighDefFreqsMat,HighDefRipMat,LowDefFreqsMat,LowDefRipMat,...
% ch_LowPassFilt,ch_HighPassFilt);
% %=========================================================================================
%  FiltMat=zeros(2,length(HighB))
%  HighA
%  HighB
%  FiltMat(1,:)=HighA;
%  FiltMat(2,:)=HighB;
%  SaveData(FiltMat,'Paul:Matlab:Plot2d3d:Plot2dUtil:FilterCoeff:High01.hf',1,[],'float64','float64');
% return
%=========================================================================================
FiltFilePath=which('High01.hf');
FiltMat=ReadData(FiltFilePath,1,[],'float64','float64');
HighA=FiltMat(1,:)
HighB=FiltMat(2,:)
%=========================================================================================

StdFak=2;
d1=[];
for i=2:length(a1)
	d1(i-1)=a1(i)-a1(i-1);
end
e1=(filtfilt(HighB,HighA,d1));
f1=e1;
StdTmp=std(e1);
MeanTmp=mean(e1);
i1=find(e1<=MeanTmp+StdFak.*StdTmp & e1>=MeanTmp-StdFak.*StdTmp);
j1=find(e1>MeanTmp+StdFak.*StdTmp | e1<MeanTmp-StdFak.*StdTmp);
figure
subplot(3,1,1)
plot(a1(2:length(a1)),d1,'r')
hold on
plot(a1(2:length(a1)),e1,'b')
if ~isempty(j1)
	plot(a1(j1),d1(j1),'r*')
	plot(a1(j1),e1(j1),'b*')
end
%=========================================================================================
d2=[];
for i=2:length(a2)
	d2(i-1)=a2(i)-a2(i-1);
end
e2=(filtfilt(HighB,HighA,d2));
f2=e2;
StdTmp=std(e2);
MeanTmp=mean(e2);
i2=find(e2<=MeanTmp+StdFak.*StdTmp & e2>=MeanTmp-StdFak.*StdTmp);
j2=find(e2>MeanTmp+StdFak.*StdTmp | e2<MeanTmp-StdFak.*StdTmp);
subplot(3,1,2)
plot(a2(2:length(a2)),d2,'r')
hold on
plot(a2(2:length(a2)),e2,'b')
if ~isempty(j2)
	plot(a2(j2),d2(j2),'r*')
	plot(a2(j2),e2(j2),'b*')
end
%=========================================================================================
d3=[];
for i=2:length(a3)
	d3(i-1)=a3(i)-a3(i-1);
end
e3=(filtfilt(HighB,HighA,d3));
f3=e3;
StdTmp=std(e3);
MeanTmp=mean(e3);
i3=find(e3<=MeanTmp+StdFak.*StdTmp & e3>=MeanTmp-StdFak.*StdTmp);
j3=find(e3>MeanTmp+StdFak.*StdTmp | e3<MeanTmp-StdFak.*StdTmp);
subplot(3,1,3)
plot(a3(2:length(a3)),d3,'r')
hold on
plot(a3(2:length(a3)),e3,'b')
if ~isempty(j3)
	plot(a3(j3),d3(j3),'r*')
	plot(a3(j3),e3(j3),'b*')
end
%=========================================================================================
figure
plot(a1(i1),1./d1(i1),'r')
hold on
plot(a2(i2),1./d2(i2),'b')
plot(a3(i3),1./d3(i3),'g')
%=========================================================================================
x1=linspace(max([a1(1),a2(1),a3(1)]),min([a1(length(a1)),a2(length(a2)),a3(length(a2))]),1000);
Int1=interp1(a1(i1),1./d1(i1),x1);
Int2=interp1(a2(i2),1./d2(i2),x1);
Int3=interp1(a3(i3),1./d3(i3),x1);
%=========================================================================================
figure
plot(Int1,'r')
hold on
plot(Int2,'b')
plot(Int3,'g')
