% Main program for analyzing Sirotem data

orient landscape;

autop = 0;		% print graphs if this = 1;

cjo = -1;		% controls # of filter passes, let = -1
			%   to get most possible for basis


% Which data sets should I analyze?

datasets = [7 8 9 0];
			% 7 = siros257.dat
			% 8 = siros258.dat
			% 9 = siros259.dat
			% 0 = siros260.dat


% Which wavelets should I use?

wvstr = 'S8';		% H1 = Haar wavelets
			% D4 = Daubechies 4
			% S8 = Symmlet 8


% ------------------------- Start Program

% Obtain lo-pass wavelet coeffs (filter for transformation) 
  clo = MakeFilter(wvstr);

% Obtain appropriate value of jo
  jo = Checkjo(cjo,clo);


% ---------------------- Perform for each specified data set

for irow = 1:length(datasets);
  id = datasets(irow);
  figure(1);  clf reset;
  figure(2);  clf reset;
  figure(3);  clf reset;
  figure(4);  clf reset;
  figure(5);  clf reset;

  % Obtain data and cut to a power of two
  [y,datastr] = LoadSiros(id);
  y = Cut2Pow2(y);
  n = length(y);

  % Obtain (stationary) wavelet coeffs of data
  [mcoefs,fcoefs] = Dswt(clo,y,jo);

  % Estimate variance of data
  sighat = medSDhat(mcoefs(:,size(mcoefs,2)));
  sigstr = num2str(sighat);

  % Denoising threshold value -- DJ's
  bthresh = sqrt(2*log(n))*sighat 

  % Obtain all the discrete wavelet transforms (different shifts)
  DWTs = DWTshifts(fcoefs(:,1),mcoefs);

  % Obtain corresponding order of decimation
  s = log2(n) - jo ;
  decim = DecimateOrder(size(DWTs,2),s);


  % Sc2Targ = sum(Y.^2)		% for checking DWT algorithm
  % Sc2DWTs = sum(DWTs.^2)


  % Plot original data, mcoeffs, and fcoeffs  (figure 1 and 2)
  SpWvCoefs(autop,y,mcoefs,fcoefs,datastr,wvstr)

  % plot p of the sorted coeffs^2 associated with each shifted 
  % DWT and find the # greater than a hard threshold value  (figure 3)	
  [nbig,p] = SpDWTc2(DWTs,bthresh,wvstr,datastr,sigstr);


% ------------------------Test different thresholds

  Tmult = 1;
  while 1
    % Choose size of threshold or go to next data set
    [flag,thresh] = ThreshSize(bthresh,Tmult);
    if flag == 1;
      break;
    end;
    d = 1;

    % Update the threshold on the sorted coeffs^2 plot  (figure 3)
    % and plot Cumulative Sum of Squares/Total Sum of Squares
    % vs. the shift #   (figure 4)
    [nbig,p] = SpDWTc2(DWTs,thresh,wvstr,datastr,sigstr);
    SpCumSc2(n,DWTs,mcoefs,y,wvstr,datastr);

    % Find "best" and "worst" dwt according to the 
    % sorted coeffs^2 plot
    [nbest,nworst] = SCompYhats(DWTs,nbig)
    best = DWTs(:,nbest);
    worst = DWTs(:,nworst);

    % Hilight the "best" and "worst" dwts on figure 3 and figure 4
    HiLiteDWT(n,p,best,nbest,thresh,'g');
    HiLiteDWT(n,p,worst,nworst,thresh,'b');

    % Reconstruct "best" and "worst" dwts based on "large" and "small"
    % coeffs and plot on figure(5)
    SpTarget(n,y,wvstr,datastr);
    SpRlargeC(best,decim,nbest,thresh,clo,jo,'g');
    SpRlargeC(worst,decim,nworst,thresh,clo,jo,'b');
    SpRsmallC(best,decim,nbest,thresh,clo,jo,'g');
    SpRsmallC(worst,decim,nworst,thresh,clo,jo,'b');


% ------------------ Examine a particular dwt more closely
  
    q2 = 1;
    while 1

      %Mouse selection of a dwt for further study
      [button,ndwt,dwt] = chooseDWT(jo,n,DWTs);
      if button == 3;
	ndwt = [ ];
        break
      end;

      % Re-draw all plots if necessary
      SpWvCoefs(autop,y,mcoefs,fcoefs,datastr,wvstr);
      if q2 ~= 1
        [nbig,p] = SpDWTc2(DWTs,thresh,wvstr,datastr,sigstr);
        SpCumSc2(n,DWTs,mcoefs,y,wvstr,datastr);
        HiLiteDWT(n,p,best,nbest,thresh,'g');
        HiLiteDWT(n,p,worst,nworst,thresh,'b');
        SpTarget(n,y,wvstr,datastr);
        SpRlargeC(best,decim,nbest,thresh,clo,jo,'g');
        SpRlargeC(worst,decim,nworst,thresh,clo,jo,'b');
        SpRsmallC(best,decim,nbest,thresh,clo,jo,'g');
        SpRsmallC(worst,decim,nworst,thresh,clo,jo,'b');
      end;

      % Hilight the chosen dwt in the Sorted Squared Coeffs plot
      % and in the Cumulative Squared Coefs Plot
      HiLiteDWT(n,p,dwt,ndwt,thresh,'r');

      % Plot the reconstructions based on "large" and "small" coeffs of
      % the chosen dwt
      SpRlargeC(dwt,decim,ndwt,thresh,clo,jo,'r');
      SpRsmallC(dwt,decim,ndwt,thresh,clo,jo,'r');

      % Connect the mcoef and fcoef plots
      nj = size(fcoefs,2);
      SConnectDWT(n,ndwt,mcoefs,fcoefs,nj);

      disp('Red is the selected DWT.')
      disp('Blue is the "worst" DWT (according to figure 3),')
      disp('  and green is the "best." ')
      disp('      ')

      % Display the number of coeffs larger than the threshold
      d = DispNbig(nbig,ndwt,nbest,nworst,d);

      q2 = q2 + 1 ;

    end; 	% while

    disp('Blue is the "worst" DWT (according to figure 3),')
    disp('  and green is the "best." ')
    disp('      ')

    % Display the number of coeffs larger than the threshold
    d = DispNbig(nbig,ndwt,nbest,nworst,d);

  end;		% while

end;		% for








