% Demo of wavelet-based function estimation
clear all
close all
% (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)


% (ii) Add noise of size 0.1. Make sure the noise is fixed 
% by fixing the seed
randn('seed',1)
sign = sig + 0.1 * randn(size(sig));

% (iii) You may plot the noisy signal here.

figure(2)
plot(t, sign)

% (iv) 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]; 

% (v) Transfer the signal in the wavelet domain.
% Choose L=8, eight levels of decomposition

sw = dwtr(sign, 8, filt);

% At this point you may view the sw. Is it disbalanced?
% Is it decorrelated?

%(vi) Let's now threshold the small coefficients.
%The universal threshold is determined as
%lambda = sqrt(2 * log(1024)) * 0.1
%   ans = 0.3723

swt=sw .* (abs(sw) > 0.3723 );

% (vii) Return now thresholded object back to the time
% domain. Of course with the same filter and L.

a=idwtr(swt,8, filt);

% (viii) Check if you made a good estimate...

figure(3)
plot(t, a, '-')