%
%  This program implements truncated singular value decomposition 
%  for solving the inverse problem of electrocardiography.
%  This version is for the endocardial data
%
%  Outputs: 
%      est = the estimated endocardial electrograms
%  Inputs:
%      Nsvd = the number of expansion modes to use
%
function est = probe_tsvd(Nsvd)
%  
%  read in the necessary data
%
%  the transfer matrix
%
  load ZHB.ASCII
  Z = ZHB;
%
%  TMAT selects those points on the endocardium that correspond to 
%  where we have sensors
%
  load TMAT
%
%  load in the endocardial data measured at the sensors
%
  load endo_data
  endo = endo_data';
%  
%  load in the probe data
%
  load probe_data
  probe = probe_data';
%
%  determine some of the parameters we need
%
  N_electrodes = 66;
  N = 201;
%
%  compute the svd of Z
%
  [U,S,V] = svd(Z);
  [m,p] = size(Z);
  temp_s = diag(S);
  s = temp_s(1:m);
  s2 = s.*s;
  tV = V(:,1:Nsvd);
%
%  plot the singular values
%
  plot([1:m],s,'x'); grid; xlabel('Mode'); ylabel('Singular Values');
%
  est = [];
%
%  compute the endocardial estimate at each time
%
  for time = 1:N
    up = U'*probe(:,time);
    top = s.*up;
    bot = s.*s;
    alpha = (top./bot);
    full_est = tV*alpha(1:Nsvd);
    est(:,time) = TMAT*full_est;
  end;
%
%  Now make the epicardial estimate plot
% 
  time = [1:N];
  ave_cc = plot_15(est,endo,time,Nsvd);
  fprintf(' The average correlation coefficient for endoardial electrograms is %f \n', ave_cc);
%  
return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Plot the estimated endocardial electrograms
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
function ave_cc = plot_15(ya,ytrue,time,Nsvd);
%
%  
  scale = 1.2;
%
%  plot the data
%
  figure;
  fs = 14;
%   whitebg(1,[0 0 0]);
%   set(1,'Color',[0 0 0]);
%   set(1,'InvertHardCopy','off');
%
  orient landscape;
%
  total_cc = [];
  for k=1:4
      p = k;
      for j=1:4
         ii = (k-1)*4+j;
         if (ii ~= 16)
           ymin = scale*min(ytrue(ii,:));
           ymax = scale*max(ytrue(ii,:));
           subplot(4,4,p);
           axis([time(1) time(end) ymin ymax]);
           hold on
           plot(time,ya(ii,:),'b-','Linewidth',1.5);
           plot(time,ytrue(ii,:),'r-','Linewidth',1.5);
           axis off;
           cc = corrcoef(ya(ii,:),ytrue(ii,:));
           total_cc = [total_cc cc(1,2)];
           title([num2str(cc(1,2),2)]);
           hold off
         end;
         p = p+4;
      end;
  end;  
  text(10,-12,['N_{svd}  = ' num2str(Nsvd)]);
  print -deps2 endo.eps
  ave_cc = mean(total_cc);
  return;
