%
%  This function allows us to play with svd
%  The signal to noise ratio (SNR) is in dB
%
function  hwk2a(SNR);
%
  Z = [-0.2042 0.0559 0.3715 -0.2077;
       0.1347 0.2208 0.5294 0.1336;
       0.4803 0.3889 0.6862 0.4670;
       0.6078 0.6088 1.2159 0.6081];
%
% open an output file
%
  fid = fopen('hwk2a.dat','wt');
%
  fprintf(fid,' File hwk2a.m, SNR = %f dB \n\n', SNR); 
%
%  compute a forward solution
%
  X = [1;1;1;1];
  Y = Z * X;
  fprintf(fid,'The true forward solution (Y) is: \n');
  fprintf(fid,' %f \n',Y);
%
%  now compute the zero noise inverse
%
  x = inv(Z)*Y;  
  fprintf(fid,'\nThe estimated value of x with no noise is: \n');
  fprintf(fid,' %f \n',x);
%
%  now compute the singular value
%  decomposition and try again
%
  [U,S,V]=svd(Z);
  beta = U'*Y;
  beta_no_noise = beta;
  alpha = inv(S)*beta;
  x = V*alpha;
  fprintf(fid,'\nThe (no noise) svd estimate of x is \n');
  fprintf(fid,' %f \n',x);
%
%  set the random number generator to a new seed each time
%
   randn('state',sum(100*clock'))
%
%  now add noise to the observation
%
%  SNR = 20*log_10(RMS/sigma)
%
  RMS = norm(Y)/sqrt(length(Y));
  sigma = RMS/(10^(0.05*SNR));
  noise = sigma*randn(length(Y),1);
%
%  determine the new y with the noise
%
  Y_noise = Y + noise;
  fprintf(fid,'\nThe new observed data (Y) with added noise is: \n');
  fprintf(fid,' %f \n',Y_noise);
%
%  determine the inverse based on the new y
%
  x = inv(Z)*Y_noise;
  fprintf(fid,'\nThe estimated x based on the inverse of Z is: \n');
  fprintf(fid,' %f \n',x);
%
%  determine the corresponding forward computed value
%
  Y = Z*x;
  fprintf(fid,'\nThe corresponding Y is: \n');
  fprintf(fid,' %f \n', Y');
%
%  now let's look at the singular values
%
  sing_val = diag(S);
  fprintf(fid,'\nThe singular values are: \n');
  fprintf(fid,' %f \n',sing_val');
  beta = U'*Y_noise;
  fprintf(fid,'\ndelta beta/sigma is (absolute values) \n');
  fprintf(fid,' %f \n',abs(beta-beta_no_noise)./sing_val);
%
%  first remove the last coefficient
%
  beta(4) = 0;
  alpha = inv(S)*beta;
  x = V*alpha;
  Y = Z*x;
  fprintf(fid,'\nThe estimate of x removing the last singular mode is: \n');
  fprintf(fid,' %f \n',x);
  fprintf(fid,'\nThe corresponding Y is: \n');
  fprintf(fid,' %f \n',Y); 
%
%  now remove the last two coefficicents
%
  beta(3) = 0;
  alpha = inv(S)*beta;
  x = V*alpha;
  Y = Z*x;
  fprintf(fid,'\nThe estimate of x removing the last two singular modes is: \n');
  fprintf(fid,' %f \n',x);
  fprintf(fid,'\nThe corresponding Y is: \n');
  fprintf(fid,' %f \n',Y);
%    
%  all done for now
%
  fclose(fid);
  return;
