%--------------------------------------
% BAMS (Bayesian Adaptive Multiscale Smoother)
%--------------------------------------
    % ver 1.1
    % This is a demo of BAMS method
    % Authors: Brani Vidakovic and Fabrizio Ruggeri
    % Blue Devil's Wavelet Group, February 2000
    % Needs: Wavelab 0.8X
    %--------------------------------------
%=== Initialization ===================
    clear all; 
    error = [ ]; 
    %======================================
  choice=menu('Which Signal?', 'Blocks', 'Bumps', ...
        'Doppler', 'HeaviSine', 'Exit BAMS');
  switch(choice);
      case(1), signal_type = 'Blocks'; wavelet = 'Haar'; par=[];
      case(2), signal_type = 'Bumps'; wavelet = 'Daubechies'; par=6;
      case(3), signal_type = 'Doppler'; wavelet = 'Symmlet'; par=8;
      case(4), signal_type = 'HeaviSine'; wavelet = 'Symmlet'; par=8;
      case(5), close all;
  end

  N = input('Input length of the sample, must be a power of 2 [1024]:');
    if isempty(N)
      N=1024;
    end
  SNR = input('Input Signal-to-Noise Ratio [7]:');
    if isempty(SNR)
      SNR=7;
    end
  KK = input('Input the number of simulational runs [10]:');
    if isempty(KK)
      KK=10;
    end
    smodata = 0.*linspace(1,N,N);
    %======================================
for k  = 1:KK,
    seed=k;        %set rn generator seed as an index; print the cycle
    disp(['run # ' num2str(k) '/'  num2str(KK)]);
    %----------------------
    coarsest=floor(log2( log (N))) + 1; %coarsest level
    randn('state', seed );              %initialize random number generator
    %--------------------------------------------------------------
    signal = MakeSignal(signal_type, N);
    sigma = std(signal);
    signal = signal .* SNR/sigma;
    %--------------------------------------------------------------
    nois = randn(N,1,1);    
    data = signal + nois';
    sdata = std(data);
    %==============================================================
    xaxis=(1:N)/N;
    figure(1)
    clf
    set(0, 'DefaultAxesFontSize', 15);
    subplot(321)
        plot(xaxis, signal)
        title('signal')
    subplot(322)
        plot(xaxis, data)
        title('noisy signal')
    %-------------------------------
    filt = MakeONFilter( wavelet, par );
    data_wd = FWT_PO (data, coarsest, filt);
    data_smowd = data_wd;
    finest_lev=log2(N)-1;
    %----------setting the parameters: mu----------------------
    niz=data_wd(dyad(finest_lev));
    nizo=sort(niz); mm=length(nizo);
    low=floor(1/4 * mm); high=floor(3/4 * mm);
    pseudos = abs(nizo(high)-nizo(low))/1.5; %hoaglin, mosteller, tukey: book
    mu = 1/pseudos^2;                        %Robust EDA, Wiley 1983
    %----------setting the parameters: tau and epsilon (levelwise)-------
    shape=ones(1, finest_lev);
    tauu= sqrt((std(data).^2 - 1/mu));
    for i = finest_lev:-1:coarsest,
          eps(i) = 1 - 1/(i - coarsest +1)^1.5;
          tau(i) = tauu;
    %----------------------------------------------------------
           a = data_wd(dyad(i));
           b = 0 .* a;
           nn=length(a);
                  for  j = 1:nn,
                    [zz,  b(j)] = bayesrule(a(j), mu, tau(i), eps(i));
                  end
           data_smowd(dyad(i))= b; 
    end
    smo_data = IWT_PO( data_smowd, coarsest, filt);
    smodata = smodata + smo_data;  %needed for calculating the bias
%----------------------------------------------------------
    %   more  plots and calculations
%----------------------------------------------------------
    subplot(323)
           plot(xaxis, smo_data)
           title('denoised signal')
    subplot(324)
           plot(xaxis, signal - smo_data, 'g-')
           title('residuals')
           x = -10:0.1:10; y =0.*x;  w=0.*x;
           mm=length(x);
           for i = finest_lev:-1:coarsest,
                   for m = 1:mm,
                   [y(m), w(m)] = bayesrule(x(m), mu, tau(i), eps(i));
           end
    subplot(325)
           plot (x,y);
           axis([-10 10 0 0.8])
           hold on
           title('marginal distributions')
    subplot(326)
           plot (x,w);
           hold on
           plot (x, x, 'r --');
           axis([-5 5 -5 5])
           title('bayes shrinkage rules')
    end
        error = [error  1/N * sum( (signal - smo_data).^2 )];
    end
    smodata=smodata/KK;
%==== OUTPUTS on the screen [BIAS-SQUARED, VAR, MSE, STDERR]=======
    biassq = 1/N * sum((signal - smodata).^2);
    variance = mean(error) - biassq;
    mse= mean(error);
    std_se = std(error);
    disp('=====================');
    disp(['Bias squared=' num2str(biassq)]);
    disp(['Variance=' num2str(variance)]);
    disp(['MSE=' num2str(mse)]);
    disp(['STDERR=' num2str(std_se)]);
%-------------------------------------------------------------------
%
% Copyright (c) 2000  B. Vidakovic and F. Ruggeri
%