
% 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, 5, filt);

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

%(vi) Let's now apply Bayesian Shrinkage. 
% Assume that the likelihood of a detail wavelet coefficient is
% normal (theta, sigma^2) and that the prior on theta
% is also normal (0, tau^2). The Bayes rule is:
%tau^2/(tau^2 + sigma^2) d. Take $tau^2=0.01 and 
%sigma^2=0.1.
swt = sw;
swt(2^5+1:end) = swt(2^5+1:end).*0.01/0.11;

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

a=idwtr(swt,5, filt);

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

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