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

NPoints=165000;
Mat=zeros(NPoints,3);
% [File,Path]=uigetfile('*.00A','Open EKG VPM file: ');
% FilePath=[Path,File];
FilePath='Paul:256Data:03:VPM:003.EKG'
                                                                      
[fid, message]=fopen(FilePath,'r','l');
fseek(fid,26,-1);
Mat(:,1)=fread(fid,NPoints,'int16');
fclose(fid);
FilePath=[FilePath(1:length(FilePath)-3),'EDA']
[fid, message]=fopen(FilePath,'r','l');
fseek(fid,26,-1);
Mat(:,2)=fread(fid,NPoints,'int16');
fclose(fid);
FilePath=[FilePath(1:length(FilePath)-3),'TRI']
[fid, message]=fopen(FilePath,'r','l');
fseek(fid,26,-1);
Mat(:,3)=fread(fid,NPoints,'int16');
fclose(fid);
Mat(:,1)=-1.*(Mat(:,1)-mean(Mat(1:1000,1)));
Mat(:,2)=Mat(:,2)-mean(Mat(1:1000,2));
Mat(:,3)=Mat(:,3)-mean(Mat(1:1000,3));
% for i=1:80
% 	plot(Mat((i-1).*5000+1:i.*5000,1))
% 	pause
% end
%===========================================
EcgVec=zeros(2000,1);
NEcg=0;
for i=2:NPoints
	if Mat(i,1)>700 & Mat(i-1,1)<700
		NEcg=NEcg+1;
		EcgVec(NEcg)=i;
	end
end
EcgVec=EcgVec(1:NEcg);
DiffECGVec=EcgVec(2:length(EcgVec))-EcgVec(1:length(EcgVec)-1);
% plot(DiffECGVec)
% %axis([100 200 225 325])
% pause
%===========================================

EDAVec=Mat(:,2);
% plot(EDAVec)
% %axis([100 200 225 325])
% pause


%===========================================
TrigVec=zeros(2000,1);
NTrig=0;
for i=2:NPoints
	if Mat(i,3)<-500 & Mat(i-1,3)>-500
		NTrig=NTrig+1;
		TrigVec(NTrig)=i;
	end
end
TrigVec=TrigVec(1:NTrig);
for i=2:37
	DiffVec(i-1)=TrigVec(i)-TrigVec(i-1);
end

if DiffVec(35)<82 | DiffVec(35)>85 | DiffVec(36)<164 | DiffVec(36)>170
	error('trigger problems')
end
StartTrig=37;

NegHigh=[2 9  13 20 23 30 33 39 42 49 52 59 63 70 72 80 83 89 92 98 ];
NegLow =[4 8  12 19 22 28 34 40 43 50 53 60 65 67 74 79 82 87 95 100];
Neutral=[1 6  11 16 21 26 31 36 41 46 51 56 61 66 71 76 81 86 91 96 ];
PosLow =[5 7  14 17 25 27 35 37 44 47 55 58 64 68 75 78 85 88 94 99 ];
PosHigh=[3 10 15 18 24 29 32 38 45 48 54 57 62 69 73 77 84 90 93 97 ];

NegHighInd=[];
NegLowInd=[];
NeutralInd=[];
PosLowInd=[];
PosHighInd=[];

NegHighEcgDiff=[];
NegLowEcgDiff=[];
NeutralEcgDiff=[];
PosLowEcgDiff=[];
PosHighEcgDiff=[];

NegHighEDA=[];
NegLowEDA=[];
NeutralEDA=[];
PosLowEDA=[];
PosHighEDA=[];

for BlockInd=1:50
	BlockInd=BlockInd
	BlockStart=TrigVec(StartTrig+(BlockInd-1).*18)+1;
	BlockEnd=TrigVec(StartTrig+BlockInd.*18);
	BlockLength=BlockEnd-BlockStart+1
	if BlockLength>3005 | BlockLength<2995
		error('wrong BlockLength')
	end
	TmpInd=[BlockStart:BlockEnd];
	EcgInd=find(EcgVec>BlockStart & EcgVec<BlockEnd);
	EcgDiff=4.*(EcgVec(EcgInd(2:length(EcgInd)))-EcgVec(EcgInd(1:length(EcgInd)-1)))';
	MeanEDA=mean(EDAVec(BlockStart+500:BlockEnd+500));
	if any(NegHigh==BlockInd)
		NegHighInd=[NegHighInd,TmpInd];
		NegHighEcgDiff=[NegHighEcgDiff,EcgDiff];
		NegHighEDA=[NegHighEDA,MeanEDA];
	elseif any(NegLow==BlockInd)
		NegLowInd=[NegLowInd,TmpInd];
		NegLowEcgDiff=[NegLowEcgDiff,EcgDiff];
		NegLowEDA=[NegLowEDA,MeanEDA];	
	elseif any(Neutral==BlockInd)
		NeutralInd=[NeutralInd,TmpInd];
		NeutralEcgDiff=[NeutralEcgDiff,EcgDiff];
		NeutralEDA=[NeutralEDA,MeanEDA];		
	elseif any(PosLow==BlockInd)
		PosLowInd=[PosLowInd,TmpInd];
		PosLowEcgDiff=[PosLowEcgDiff,EcgDiff];
		PosLowEDA=[PosLowEDA,MeanEDA];		
	elseif any(PosHigh==BlockInd)
		PosHighInd=[PosHighInd,TmpInd];
		PosHighEcgDiff=[PosHighEcgDiff,EcgDiff];
		PosHighEDA=[PosHighEDA,MeanEDA];		
	end
end
hEDAFig=figure('Units','normal','Position',[0 .5 .5 .5],'Name','EDA by Condition','NumberTitle','off');
plot(EDAVec)
hold on
plot(NegHighInd,EDAVec(NegHighInd),'ro')
plot(PosHighInd,EDAVec(PosHighInd),'bo')
plot(NeutralInd,EDAVec(NeutralInd),'ko')
plot(NegLowInd,EDAVec(NegLowInd),'mo')
plot(PosLowInd,EDAVec(PosLowInd),'go')
axis([0 20000 min(EDAVec) max(EDAVec)]) %length(EDAVec)

% hECGFig=figure('Units','normal','Position',[0 0 .5 .5],'Name','ECG by Condition','NumberTitle','off');
% plot(DiffECGVec)
% hold on
% plot(NegHighInd,DiffECGVec(NegHighInd),'ro')
% plot(PosHighInd,DiffECGVec(PosHighInd),'bo')
% plot(NeutralInd,DiffECGVec(NeutralInd),'ko')
% plot(NegLowInd,DiffECGVec(NegLowInd),'mo')
% plot(PosLowInd,DiffECGVec(PosLowInd),'go')
% axis([0 20000 min(DiffECGVec) max(DiffECGVec)]) %length(DiffECGVec)

MeanNegHighEcgDiff=mean(NegHighEcgDiff)
MeanNegLowEcgDiff=mean(NegLowEcgDiff)
MeanNeutralEcgDiff=mean(NeutralEcgDiff)
MeanPosLowEcgDiff=mean(PosLowEcgDiff)
MeanPosHighEcgDiff=mean(PosHighEcgDiff)

MedNegHighEcgDiff=median(NegHighEcgDiff)
MedNegLowEcgDiff=median(NegLowEcgDiff)
MedNeutralEcgDiff=median(NeutralEcgDiff)
MedPosLowEcgDiff=median(PosLowEcgDiff)
MedPosHighEcgDiff=median(PosHighEcgDiff)

MeanNegHighEDA=mean(NegHighEDA)
MeanNegLowEDA=mean(NegLowEDA)
MeanNeutralEDA=mean(NeutralEDA)
MeanPosLowEDA=mean(PosLowEDA)
MeanPosHighEDA=mean(PosHighEDA)

MedNegHighEDA=median(NegHighEDA)
MedNegLowEDA=median(NegLowEDA)
MedNeutralEDA=median(NeutralEDA)
MedPosLowEDA=median(PosLowEDA)
MedPosHighEDA=median(PosHighEDA)

hEDAEcgFig=figure('Units','normal','Position',[.5 .5 .5 .5],'Name','ECG & EDA','NumberTitle','off');

subplot(2,2,1)
plot(NegHighEcgDiff,'r')
hold on
plot(MedNegHighEcgDiff.*ones(size(NegHighEcgDiff)),'r--')
plot(NegLowEcgDiff,'m')
plot(MedNegLowEcgDiff.*ones(size(NegLowEcgDiff)),'m--')
plot(NeutralEcgDiff,'w')
plot(MedNeutralEcgDiff.*ones(size(NeutralEcgDiff)),'k--')
plot(PosLowEcgDiff,'g')
plot(MedPosLowEcgDiff.*ones(size(PosLowEcgDiff)),'g--')
plot(PosHighEcgDiff,'b')
plot(MedPosHighEcgDiff.*ones(size(PosHighEcgDiff)),'b--')
axis([1 100 800 1400])

subplot(2,2,3)
plot([MedPosHighEcgDiff MedPosLowEcgDiff MedNeutralEcgDiff MedNegLowEcgDiff MedNegHighEcgDiff])

subplot(2,2,2)
plot(NegHighEDA,'r')
hold on
plot(MedNegHighEDA.*ones(size(NegHighEDA)),'r--')
plot(NegLowEDA,'m')
plot(MedNegLowEDA.*ones(size(NegLowEDA)),'m--')
plot(NeutralEDA,'w')
plot(MedNeutralEDA.*ones(size(NeutralEDA)),'k--')
plot(PosLowEDA,'g')
plot(MedPosLowEDA.*ones(size(PosLowEDA)),'g--')
plot(PosHighEDA,'b')
plot(MedPosHighEDA.*ones(size(PosHighEDA)),'b--')
axis([1 10 350 500])

subplot(2,2,4)
plot([MedPosHighEDA MedPosLowEDA MedNeutralEDA MedNegLowEDA MedNegHighEDA])




		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
		
