%DCMDDENO Demonstrates denoising tools in the Wavelet Toolbox. 

%       M. Misiti, Y. Misiti, G. Oppenheim, J.M. Poggi 12-Mar-96.
%	Last Revision: 01-Oct-97.
%       Copyright (c) 1995-98 by The MathWorks, Inc.
% $Revision: 1.6 $

% Demo initialization ====================
if ~exist('DRVShowGUIFlag'), figNumber = 0; end;

echo off
clc;
titre   = '1-D De-noising using wavelet';

if wshowdrv('init',figNumber)
        if figNumber==0
		active_fig = wdfigutl('create',titre);
        else
                set(figNumber,'Name',titre);
                wshowdrv('setactfig',figNumber,figNumber);
                democmdm('startshow',figNumber,'deno');
                return
        end
end
if exist('wait_mode')~=1 , wait_mode = 1; end
if figNumber , active_fig = wshowdrv('getactfig',figNumber); end

comment = str2mat( ...
                '% Set signal to noise ratio and set rand seed.', ...
                '        snr = 3; init = 2055615866;', ...
                ' ', ...
                '% Generate original signal and a noisy version adding',...
                '% a standard Gaussian white noise.',...
                '        [xref,x] = wnoise(3,11,snr,init);'...
          );
wshowdrv('disp',figNumber,comment); 

% Set signal to noise ratio and set rand seed.
        snr = 3; init = 2055615866;
  
% Generate original signal and a noisy version adding
% a standard Gaussian white noise.
        [xref,x] = wnoise(3,11,snr,init);

% Plot signals.
figure(active_fig)
h1 = subplot(311);
plot(xref,'r'), axis([1 2048 -10 10]);
title('Original signal');
figure(active_fig)
h2 = subplot(312);
plot(x,'g'), axis([1 2048 -10 10]);
title(['Noisy signal - Signal to noise ratio = ',num2str(fix(snr))]);
set(h1,'Xtick',[],'Xticklabel',[]);
if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
                '% Denoise noisy signal using soft heuristic SURE thresholding', ...
                '% and scaled noise option, on detail coefficients obtained', ...
                '% from the decomposition of x, at level 5 by sym8 wavelet.', ...
                '% Generate original signal and a noisy version adding', ...
                '% a standard Gaussian white noise.',...
                '        lev = 5;', ...
                '        xd  = wden(x,''heursure'',''s'',''one'',lev,''sym8'');'...
                );
wshowdrv('disp',figNumber,comment); 

% Denoise noisy signal using soft heuristic SURE thresholding
% and scaled noise option, on detail coefficients obtained
% from the decomposition of x, at level 5 by sym8 wavelet.
        lev = 5;
        xd  = wden(x,'heursure','s','one',lev,'sym8');

% Plot signal.
figure(active_fig)
set(h2,'Xtick',[],'Xticklabel',[]);
subplot(313); plot(xd); axis([1 2048 -10 10])
title('Denoised signal - heuristic SURE');
if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
                '% Denoise noisy signal using soft SURE thresholding.', ...
                ' ', ...
                '        xd = wden(x,''rigrsure'',''s'',''one'',lev,''sym8'');'...
                );
wshowdrv('disp',figNumber,comment);

% Denoise noisy signal using soft SURE thresholding.
        xd = wden(x,'rigrsure','s','one',lev,'sym8');

% Plot signal.
figure(active_fig)
subplot(313); plot(xd); axis([1 2048 -10 10]);
title('Denoised signal - SURE');
if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
                '% Denoise noisy signal using fixed form threshold with', ...
                '% a single level estimation of noise standard deviation.', ...
                ' ', ...
                '        xd = wden(x,''sqtwolog'',''s'',''sln'',lev,''sym8'');'...
                );
wshowdrv('disp',figNumber,comment);

% Denoise noisy signal using fixed form threshold with
% a single level estimation of noise standard deviation.
        xd = wden(x,'sqtwolog','s','sln',lev,'sym8');

% Plot signal.
figure(active_fig)
subplot(313); plot(xd); axis([1 2048 -10 10])
title('Denoised signal - Fixed form threshold');

if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
                '% Denoise noisy signal using fixed minimax threshold with', ...
                '% a multiple level estimation of noise standard deviation.', ...
                ' ', ...
                '        xd = wden(x,''minimaxi'',''s'',''sln'',lev,''sym8'');'...
                );
wshowdrv('disp',figNumber,comment);

% Denoise noisy signal using minimax threshold with
% a multiple level estimation of noise standard deviation.
        xd = wden(x,'minimaxi','s','sln',lev,'sym8');

% Plot signal.
figure(active_fig)
subplot(313); plot(xd); axis([1 2048 -10 10])
title('Denoised signal - Minimax');
comment = str2mat( ...
                '% If many trials are necessary, it is better to perform', ...
                '% decomposition one time and threshold it many times :', ...
                ' ', ...
                '        % decomposition.', ...
                '        [c,l] = wavedec(x,lev,''sym8'');', ...
                ' ', ...
                '        % threshold the decomposition structure [c,l].', ...
                 '        xd = wden(c,l,''minimaxi'',''s'',''sln'',lev,''sym8'');'...
                );
wshowdrv('disp',figNumber,comment);

% If many trials are necessary, it is better to perform 
% decomposition one time and threshold it many times:
   % decomposition.
   [c,l] = wavedec(x,lev,'sym8');
   % threshold the decomposition structure [c,l].
   xd = wden(c,l,'minimaxi','s','sln',lev,'sym8');

if wshowdrv('pause',figNumber,wait_mode), return; end;

% Reset axes.
%------------
delete(findobj(active_fig,'type','axes'));

comment = str2mat( ...
                '% Load electrical signal and select a part.', ...
                ' ', ...
                '        load leleccum; indx = 2600:3100;', ...
                '        x = leleccum(indx);' ...
                );
wshowdrv('disp',figNumber,comment);

% Load electrical signal and select a part.
        load leleccum; indx = 2600:3100; 
        x = leleccum(indx);

figure(active_fig)
subplot(211);
plot(indx,x,'r'), title('Original signal')

if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
                '% Use wdencmp for signal de-noising.', ...
                ' ', ...
                '       % find default values (see ddencmp).', ...
                '        [thr,sorh,keepapp] = ddencmp(''den'',''wv'',x);', ...
                ' ', ...
                '        % denoise signal using global thresholding option.', ...
                '        xd = wdencmp(''gbl'',x,''db3'',2,thr,sorh,keepapp);' ...
                );
wshowdrv('disp',figNumber,comment);

% Use wdencmp for signal de-noising.
        % find default values (see ddencmp).
        [thr,sorh,keepapp] = ddencmp('den','wv',x);

        % denoise signal using global thresholding option.
        xd = wdencmp('gbl',x,'db3',2,thr,sorh,keepapp);

figure(active_fig)
subplot(212);
plot(indx,xd), title('Denoised signal')
if wshowdrv('pause',figNumber,wait_mode), return; end;

% Reset axes.
%------------
delete(findobj(active_fig,'type','axes'));

comment = str2mat( ...
                '% Some trial examples without commands counterpart.' ...
                );
wshowdrv('disp',figNumber,comment);

init = 2055615866;
snr  = 5;  % signal to noise ratio.

[xref,x] = wnoise(1,11,snr,init);
indx = linspace(0,1,length(x));
figure(active_fig)
subplot(211);
plot(indx,x,indx,xref)
title('Noisy and Original signals')
xlabel(['Signal to noise ratio = ',num2str(fix(snr))])
xd = wden(x,'heursure','s','one',5,'sym8');
figure(active_fig)
subplot(212);
plot(indx,xd,indx,xref), xlabel('Denoised and Original signals')
if wshowdrv('pause',figNumber,wait_mode), return; end;

% Reset axes.
%------------
delete(findobj(active_fig,'type','axes'));

snr = 4;  % signal to noise ratio.
[xref,x] = wnoise(2,11,snr,init);
indx = linspace(0,1,length(x));
figure(active_fig)
subplot(211);
plot(indx,x,indx,xref)
title('Noisy and Original signals')
xlabel(['Signal to noise ratio = ',num2str(fix(snr))])
xd = wden(x,'rigrsure','s','one',5,'sym4');
figure(active_fig)
subplot(212);
plot(indx,xd,indx,xref), xlabel('Denoised and Original signals')
if wshowdrv('pause',figNumber,wait_mode), return; end;

% Reset axes.
%------------
delete(findobj(active_fig,'type','axes'));

snr = 3;  % signal to noise ratio.
[xref,x] = wnoise(3,11,snr,init);
indx = linspace(0,1,length(x));
figure(active_fig)
subplot(211);
plot(indx,x,indx,xref)
title('Noisy and Original signals')
xlabel(['Signal to noise ratio = ',num2str(fix(snr))])
xd = wden(x,'sqtwolog','s','one',5,'sym8');
figure(active_fig)
subplot(212);
plot(indx,xd,indx,xref), xlabel('Denoised and Original signals')
if wshowdrv('pause',figNumber,wait_mode), return; end;

% Reset axes.
%------------
delete(findobj(active_fig,'type','axes'));

snr = 2;  % signal to noise ratio.
[xref,x] = wnoise(3,11,snr,init);
indx = linspace(0,1,length(x));
figure(active_fig)
subplot(211);
plot(indx,x,indx,xref)
title('Noisy and Original signals')
xlabel(['Signal to noise ratio = ',num2str(fix(snr))])
xd = wden(x,'sqtwolog','s','one',5,'sym8');
figure(active_fig)
subplot(212);
plot(indx,xd,indx,xref), xlabel('Denoised and Original signals')
if wshowdrv('pause',figNumber,wait_mode), return; end;

% Reset axes.
%------------
delete(findobj(active_fig,'type','axes'));

snr = 4;  % signal to noise ratio.
[xref,x] = wnoise(4,11,snr,init);
indx = linspace(0,1,length(x));
figure(active_fig)
subplot(211);
plot(indx,x,indx,xref)
title('Noisy and Original signals')
xlabel(['Signal to noise ratio = ',num2str(fix(snr))])
xd = wden(x,'minimaxi','s','one',5,'sym4');
figure(active_fig)
subplot(212);
plot(indx,xd,indx,xref), xlabel('Denoised and Original signals')
if wshowdrv('pause',figNumber,wait_mode), return; end;

% Reset axes.
%------------
hdl_a = wshowdrv('set_axes',figNumber,[2 2]);

if figNumber , fig = figNumber; else , fig = active_fig; end
set(fig,'Name','2-D De-noising using wavelet');

comment = str2mat( ...
                '% Load original image.', ...
                ' ', ...
                '        load sinsin', ...
                '        % X contains the original image.'...
                );
wshowdrv('disp',figNumber,comment);

% Load original image.
        load  sinsin
% Dim axes
row = 2; col = 2;

% Plot Original Image.
figure(active_fig)
colormap(pink)
sm = size(map,1);
h1 = subplot(row,col,1);
image(wcodemat(X,sm)), title('Original Image')
axis('square')
set(h1,'Xtick',[],'Xticklabel',[],'Ytick',[],'Yticklabel',[]);
if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
                '% Generate noisy image.', ...
                ' ', ...
                '        init = 2055615866; randn(''seed'',init);', ...
                '        x = X + 18*randn(size(X));' ...
                );
wshowdrv('disp',figNumber,comment);

% Generate noisy image. 
        init = 2055615866; randn('seed',init); 
        x = X + 18*randn(size(X));

% Plot Noisy Image.
figure(active_fig)
h2 = subplot(row,col,2);
image(wcodemat(x,sm)), title('Noisy Image')
axis('square')
set(h2,'Xtick',[],'Xticklabel',[],'Ytick',[],'Yticklabel',[]);
if wshowdrv('pause',figNumber,wait_mode), return; end;

comment = str2mat( ...
                '% Use wdencmp for image de-noising.', ...
                ' ', ...
                '        % find default values (see ddencmp).', ...
                '        [thr,sorh,keepapp] = ddencmp(''den'',''wv'',x);', ...
                ' ', ...
                '        % denoise image using global thresholding option.', ...
                '        xd = wdencmp(''gbl'',x,''sym4'',2,thr,sorh,keepapp);' ...
                );
wshowdrv('disp',figNumber,comment);

% Use wdencmp for image de-noising.
        % find default values (see ddencmp).
        [thr,sorh,keepapp] = ddencmp('den','wv',x);

        % denoise image using global thresholding option.
        xd = wdencmp('gbl',x,'sym4',2,thr,sorh,keepapp);

% Plot De-noised Image.
figure(active_fig)
h3 = subplot(row,col,3);
image(wcodemat(xd,sm)), title('De-noised Image')
axis('square')
set(h3,'Xtick',[],'Xticklabel',[],'Ytick',[],'Yticklabel',[]);

if figNumber==0
        wshowdrv('pause',figNumber,wait_mode);
        delete(active_fig);
end
