% Demo of wavelet-based function estimation
clear all
close all
        %--- set plotting parameters -------
        lw = 2;  
        set(0, 'DefaultAxesFontSize', 15);
        fs = 15;
        msize = 6;

% (i) Make "Doppler" signal on [0,1]
t=linspace(0,1,1024);
sig = sqrt(t.*(1-t)).*sin((2*pi*1.05) ./(t+.05));
% and plot it
figure(1); plot(t, sig,'linewidth',lw); axis tight
print -deps wdemo1.eps
% (ii) Add noise of size 0.1. We are fixing
% the seed of random number generator for repeatability
% of example. We add the random noise to the signal
% and make a plot.
randn('seed',1)
sign = sig + 0.1 * randn(size(sig));
figure(2); plot(t, sign,'linewidth',lw); axis tight
print -deps wdemo2.eps
% (iii) Take the filter H, in this case this is SYMMLET 4
filt = [ -0.07576571478934 -0.02963552764595 ...
0.49761866763246 0.80373875180522 ...
0.29785779560554 -0.09921954357694 ...
-0.01260396726226 0.03222310060407];
% (iv) Transform the noisy signal in the wavelet domain.
% Choose L=8, eight detail levels in the decomposition.
sw = dwtr(sign, 8, filt);
% At this point you may view the sw. Is it disbalanced?
% Is it decorrelated?
%(v) Let's now threshold the small coefficients.
%The universal threshold is determined as
%lambda = sqrt(2 * log(1024)) * 0.1 = 0.3723
%
% Here we assumed $sigma=0.1$ is known. In real life
% this is not true and we estimate sigma. A robust estimator
% is 'MAD' estimator from the finest level of detail
% believed to be an image of noise.
finest = sw(513:1024);
sigma_est = 1/0.6745 * median(abs( finest - median(finest)));
lambda = sqrt(2 * log(1024)) * sigma_est;
% hard threshold in the wavelet domain
swt=sw .* (abs(sw) > lambda );
figure(3); plot([1:1024], swt ,'linewidth',lw) ; axis tight
print -deps wdemo3.eps
% (vi) Back-transform the thresholded object to the time
% domain. Of course, retain the same filter and value L.
sig_est = idwtr(swt, 8, filt);
figure(4); plot(t, sig_est,'linewidth',lw); hold on;  axis tight
print -deps wdemo4.eps