% Matlab script to compare the inverse Dower transform with a real 
% record VCG assumes that you have an internet connection and have 
% installed WFDB libs under linux correctly. See www.physionet.org
%
% Gari D. Clifford April 2006  http://alum.mit.edu/www/gari/ecgtools


% download 10s of data form the PTBDB and save to file
!rdsamp -r ptbdb/patient001/s0010_re | head -10000 > allLeads.dat

% load the data
load allLeads.dat
% column order = t i ii iii avr avl avf v1 v2 v3 v4 v5 v6 x y z 

%find 8 lead ECG and real VCG
ecg = zeros(length(allLeads),8);
ecg(:,1:2) = allLeads(:,2:3);
ecg(:,3:8) = allLeads(:,8:13);

vcg = allLeads(:,14:16);


% Construct inverse dower transform matrix
T = [ -0.172 -0.074  0.122  0.231 0.239 0.194  0.156 -0.010 ;
       0.057 -0.019 -0.106 -0.022 0.041 0.048 -0.227  0.887 ;
      -0.229 -0.310 -0.246 -0.063 0.055 0.108  0.022  0.102 ];


% estimate VCG using inverse Dower Transform
vcg2 = T*ecg';

% compare
figure;plot(vcg2(1,:))
hold on;plot(vcg(:,1),'m')
legend('Inverse Dower X','Real X')
xlabel ('samples (ms)');
ylabel ('Amplitude');
title('Inverse Dower test on ptbdb/patient001/s0010-re')
axis([800 2000 -inf inf])
print('-dpng', 's00010_reVCGcf.png');

figure;plot(vcg2(1,:)/max(vcg2(1,:)))
hold on;plot(vcg(:,1)/max(vcg(:,1)),'m')
legend('Inverse Dower X','Real X')
xlabel ('samples (ms)');
ylabel ('Normalised Amplitude');
title('normalised Inverse Dower test on ptbdb/patient001/s0010-re')
axis([800 2000 -inf inf])
print('-dpng', 's00010_reVCGcfnorm.png');


figure;plot(vcg2')
title('Inverse Dower VCG ptbdb/patient001/s0010-re')
xlabel ('samples (ms)');
ylabel ('Amplitude');
legend('X','Y','Z')
axis([800 2000 -inf inf])
print('-dpng', 's00010_reVCGart.png');

figure;plot(vcg)
title('Recorded VCG ptbdb/patient001/s0010-re')
xlabel ('samples (ms)');
ylabel ('Amplitude');
legend('X','Y','Z')
axis([800 2000 -inf inf])
print('-dpng', 's00010_reVCGreal.png');


