%% Denoising with Wavelet Transforms
%
% *Wavelets Workshop*
%
% June 4-7, 2008
% 
% University of St. Thomas
%
% *Catherine Beneteau*, *Caroline Haddad*, *David Ruch*, *Patrick Van Fleet*

%% Objective
% The purpose of this notebook is to introduce you to some basic
% wavelet-based denoising methods.

%% Conventions
% This m-file uses the DiscreteWavelets Toolbox.  Help is available on all
% functions in the toolbox by accessing the Help menu and clicking on the
% DiscreteWavelets Toolbox.  
%
% Demos are also available.  Click the Demos tab under Help to access all
% demos.

%% The Shrinkage Function
% The shrinkage function is defined in DiscreteWavelets and called
% ShrinkageFunction.  Here is the shrinkage function plotted for lambda=2.

x=-5:.001:5;
plot(x,ShrinkageFunction(x,2));

%%
% Let's try it on an example.

v = [-3.2 -2.5 1 1 1.1 5 3.2 1.9 -4 -6];

lambda = 2;
ShrinkageFunction(v,lambda)

%% A Sample Signal
% For our example, we will use a sample signal suggested by Donoho.  He
% calls this function the Heavisine function.
%
% Let's sample this function to create our true signal v and then ListPlot
% it.

n=2048;
x=0:1/n:1-1/n;
v=Heavisine(x);
close;
plot(x,v);

%% Creating Noise
% Now we create some Gaussian white noise.

noise = randn(1, n);
close;
plot(x,noise)

%%
% If you want, you can actually "play" the noise. Please be considerate of
% your neighbors  :-) 

sound(noise,11025);

%%
% We create the noisy vector by picking a noise level and simply adding v
% to a scaled version of the noise.

sigma = 0.5;
y = v + sigma*noise;
close;
plot(x,y,'Color',[.8 .6 .3]);

%% Estimating the Noise Level
% To estimate the noise level, we use the MAD function in the
% DiscreteWavelets package.  We have to first compute the wavelet 
% transform.  We'll use the Daubechies 4 filter for this example and
% compute 2 iterations.

h = Daub(4);
disp(sprintf('We will use the Daubechies 4-term filter.'));

i = 5;
wt = WT1D(y, h,i);
close;
WaveletVectorPlot(wt,i);

%%
% We need to extract the first highpass portion.  There is a command in
% DiscreteWavelets to do this:

wtlist=WaveletVectorToList(wt,i);
close;
plot(wtlist(1).hp);

%%
% Now we estimate the noise:

sigmahat = MAD(wtlist(1).hp)/.6745;
disp(sprintf('Our estimate of the noise is %f.',sigmahat));

%% VisuShrink
% To perform VISUShrink, we compute the universal tolerance and use it in
% the shrink functions on the highpass portions.  

hplen=n-length(wtlist(1).lp);

lambda = sigmahat*sqrt(2*log(hplen));
disp(sprintf('The value of the universal threshold is %f.', lambda));


for k=1:i
    wtlist(k).hp=ShrinkageFunction(wtlist(k).hp,lambda);
end

newwt=WaveletListToVector(wtlist,i);
close;
WaveletVectorPlot(newwt, i);

%%
% We apply the inverse wavelet transform to obtain our estimate of v:

vhat = IWT1D(newwt, h, i);
close;
plot(x,vhat,'Color',[0 1 0]);

%%
% Here is a visual comparison:

subplot(3,1,1);
plot(x,v);
title('True Signal');
subplot(3,1,2);
plot(x,y,'Color',[.8 .6 .3]);
title('Noisy Signal');
subplot(3,1,3);
plot(x,vhat,'Color',[0 1 0]);
title('Denoised Signal');

%% Extra Credit Problem
% Try denoising another of 4 basic functions.  Simply replace the Heavisine
% function above with the function of your choice.  All functions should be
% sampled on the interval [0,1].  Try changing the number of iteration of 
% the transformation.  You can also use Daub(6), Daub(8), Daub(10), 
% Coif(1), Coif(2) as possible filters.

% Doppler function
v=sqrt(x.*(1-x)).*sin(2*pi*1.05./(x+.05));
close;
plot(x,v,'Color',[1 0 0]);
title('Doppler function');

%%
% Squarewave function
v = sign(cos(8*pi*x));
close;
plot(x,v,'Color',[1 0 0]);
title('Squarewave function');

% Sawtooth function
xt=x:1:511;
v = [xt xt xt xt];
close;
plot(x,v,'Color',[.5 .5 1]);
title('Sawtooth function');

%% Audio Denoising

%% Import an Audio File
% Our first step is to load one of the audio clips included with the
% DiscreteWavelets Toolbox.  To do this, we use the AudioList function to
% list all the files and the AudioNames function to create absolute paths
% to each audio file.

AudioList();
d=AudioNames();

% Let's choose the second clip (bark.wav).

%%
% Since the second clip is an audio file in wav format, we use wavread to
% import it as a vector.  We also trim some elements off the end of the
% vector to ensure that the length of the resulting vector is divisible by
% 2^4.

[data,srate]=wavread(d{2});
data=ChopVector(data,4);
N=length(data);

disp(sprintf('The length of the audio clip is %i.',N));
disp(sprintf('The sampling rate for the clip is %i.',srate));

%% Play the Audio Clip
% In this cell, we play the audio clip and plot it.

plot(data);
title('Bark.wav');

wavplay(data,srate);

%% Add White Noise to the Clip
% We next add some white noise to the clip.  We can use the Matlab function
% randn to assist in this task.

noise=randn(N,1);

% Pick a noise level.
sigma = .1;

% Create the noisy audio clip.
v=data+sigma*noise;

% Plot the noisy vector.
clf;
plot(v);

wavplay(v,srate);

%% Compute the Discrete Wavelet Transformation
% We now compute four iterations of the discrete wavelet transformation of
% v.  We will use the D6 filter for this task.

its=4;
h=Daub(6);
wt=WT1D(v,h,its);

% Plot the transformation.
clf;
WaveletVectorPlot(wt,its);
title(sprintf('Wavelet Transformation - %i Iterations',its));

%% Estimate the Noise
% The next step in the denoising process is to estimate the noise.  Of
% course, we know the noise level, but it is good to see how well our
% estimator works.  We use the formula sigma ~ MAD(hp(1))/.6745 (see 
%Section 9.1 of the text).

% First we split the wavelet transform into various parts.
wtlist=WaveletVectorToList(wt,its);

% We compute the MAD of the first highpass portion and divide by .6745.
sigmahat=MAD(wtlist(1).hp)/.6745;

disp(sprintf('The noise is estimated to be %f.',sigmahat));

%% Use SureShrink to Denoise
% We will use SureShrink to denoise the audio clip.  
%
% The vector lambda is simply a vector of length its where lambda(1) is the
% tolerance for the first highpass portion, and so on. 
%
% Feel free to choose your own vector of lambdas!
%

for j=1:its
    % The highpass portions must have variance 1.
    wtlist(j).hp=wtlist(j).hp/sigmahat;
    % Find the SureShrink tolerance
    lambda=DonohoSure(wtlist(j).hp);
    % Perform the quantization
    wtlist(j).hp=ShrinkageFunction(wtlist(j).hp,lambda);
    % Rescale the highpass portions.
    wtlist(j).hp=wtlist(j).hp*sigmahat;
end

%%
% Plot the modified wavelet transformation

% First reconstruct the modified wavelet transform vector.
newwt=WaveletListToVector(wtlist,its);

% Plot the modified wavelet transformation.
clf;
WaveletVectorPlot(newwt,its);
title('Modified Wavelet Transformation');

%% Compute the Inverse Wavelet Transform 
% To obtain the denoised audio clip, we perform four iterations of the
% inverse wavelet transformation.

denoised=IWT1D(newwt,Daub(6),its);

% Plot the denoised audio clip.
plot(denoised);
title('Denoised Audio Clip');

% Play the denoised audio clip.
wavplay(denoised,srate);

%%
% Play the noisy clip.

wavplay(v,srate);

%%
% Play the original clip.

wavplay(data,srate);

%% Exercise - For Fun
% Open the Laurel and Hardy audio file above.  It is very large so you 
% might you might want to use only say one quarter  of the file.  Play the 
% audio file - it is noisy since it is an old recording.
%
% Your task:  By whatever means possible (but hopefully using wavelets!), 
% denoise the audio clip!

%%
close all;
displayEndOfDemoMessage(mfilename)