www.gusucode.com > wavelet 源码程序 matlab案例代码 > wavelet/denoisingsignalsdemo.m

    %% Denoising Signals and Images
% This example discusses the problem of signal recovery from noisy data. The
% general denoising procedure involves three steps. The basic version of
% the procedure follows the steps described below:
%
% * Decompose: Choose a wavelet, choose a level N. Compute the wavelet
% decomposition of the signal at level N. 
%
% * Threshold detail coefficients: For each level from 1 to N, select a
% threshold and apply soft thresholding to the detail coefficients. 
%
% * Reconstruct: Compute wavelet reconstruction using the original
% approximation coefficients of level N and the modified detail
% coefficients of levels from 1 to N. 
%
% Two points must be addressed in particular: 
%
% * how to choose the threshold, 
%
% * and how to perform the thresholding. 

% Copyright 2006-2012 The MathWorks, Inc.

%% Soft or Hard Thresholding?
% Thresholding can be done using the function WTHRESH which returns soft or
% hard thresholding of the input signal. Hard thresholding is the simplest
% method but soft thresholding has nice mathematical properties. Let thr
% denote the threshold.
thr = 0.4; 
%%
% Hard thresholding can be described as the usual process of setting to
% zero the elements whose absolute values are lower than the threshold. The
% hard threshold signal is x if |x|>thr, and is 0 if |x|<=thr. 
y = linspace(-1,1,100); 
ythard = wthresh(y,'h',thr); 
%% 
% Soft thresholding is an extension of hard thresholding, first setting to
% zero the elements whose absolute values are lower than the threshold, and
% then shrinking the nonzero coefficients towards 0. The soft threshold
% signal is sign(x)(|x|-thr) if |x|>thr and is 0 if |x|<=thr. 
ytsoft = wthresh(y,'s',thr);
subplot(1,3,1)
plot(y)
title('Original')
subplot(1,3,2)
plot(ythard)
title('Hard Thresholding')
subplot(1,3,3)
plot(ytsoft)
title('Soft Thresholding')
%%
% As can be seen in the figure above, the hard procedure creates
% discontinuities at x = 
% $\pm$
% t, while the soft procedure does not. 

%% Threshold Selection Rules
% Recalling step 2 of the de-noise procedure, the function THSELECT
% performs a threshold selection, and then each level is thresholded. This
% second step can be done using WTHCOEF, directly handling the wavelet
% decomposition structure of the original signal. Four threshold selection
% rules are implemented in the function THSELECT. Typically it is
% interesting to show them in action when the input signal is a Gaussian
% white noise. 
rng default;
y = randn(1,1000);
%%
% Rule 1: Selection using principle of Stein's Unbiased Risk Estimate (SURE)
thr = thselect(y,'rigrsure')
%%
% Rule 2: Fixed form threshold equal to sqrt(2*log(length(y)))
thr = thselect(y,'sqtwolog')
%%
% Rule 3: Selection using a mixture of the first two options
thr = thselect(y,'heursure')
%%
% Rule 4: Selection using minimax principle
thr = thselect(y,'minimaxi')
%%
% Minimax and SURE threshold selection rules are more conservative and
% would be more convenient when small details of the signal lie near the
% noise range. The two other rules remove the noise more efficiently. 
%% 
% Let us use the "blocks" test data credited to Donoho and Johnstone as
% first example. Generate original signal xref and a noisy version x adding 
% a standard Gaussian white noise. 
sqrt_snr = 4;      % Set signal to noise ratio
init = 2055615866; % Set rand seed
[xref,x] = wnoise(1,11,sqrt_snr,init);
%%
% De-noise noisy signal using soft heuristic SURE thresholding on detail
% coefficients obtained from the decomposition of x, at level 3 by sym8
% wavelet. 
scal = 'one'; % Use model assuming standard Gaussian white noise.
xd = wden(x,'heursure','s',scal,3,'sym8');
Nx = length(x);
subplot(3,1,1),plot(xref), xlim([1 Nx])
title('Original Signal')
subplot(3,1,2),plot(x), xlim([1 Nx])
title('Noisy Signal')
subplot(3,1,3),plot(xd), xlim([1 Nx])
title('De-noised Signal - Signal to noise ratio = 4')
%%
% Since only a small number of large coefficients characterize the original
% signal, the method performs very well.

%% Dealing with Non-White Noise
% When you suspect a nonwhite noise, thresholds must be rescaled by a
% level-dependent estimation of the level noise. As a second example, let
% us try the method on the highly perturbed part of an electrical signal.
% Let us use db3 wavelet and decompose at level 3. To deal with the
% composite noise nature, let us try a level-dependent noise size
% estimation. 
load leleccum; indx = 2000:3450; 
x = leleccum(indx); % Load electrical signal and select part of it. 
%%
% Find first value in order to avoid edge effects. 
deb = x(1);
%%
% De-noise signal using soft fixed form thresholding 
% and unknown noise option. 
scal = 'mln'; % Use a level-dependent estimation of the level noise
xd = wden(x-deb,'sqtwolog','s',scal,3,'db3')+deb;
Nx = length(x);
subplot(2,1,1),plot(x), xlim([1 Nx])
title('Original Signal')
subplot(2,1,2),plot(xd), xlim([1 Nx])
title('De-noised Signal')
%%
% The result is quite good in spite of the time heterogeneity of the nature
% of the noise after and before the beginning of the sensor failure around
% time 2450.

%% Image Denoising
% The denoising method described for the one-dimensional case applies also
% to images and applies well to geometrical images. The two-dimensional
% denoising procedure has the same three steps and uses two-dimensional
% wavelet tools instead of one-dimensional ones. For the threshold
% selection, prod(size(y)) is used instead of length(y) if the fixed form
% threshold is used. 
%%
% Generate a noisy image.
load  woman 
init = 2055615866; rng(init); 
x = X + 15*randn(size(X));
%%
% In this case fixed form threshold is used with estimation of level noise,
% thresholding mode is soft and the approximation coefficients are kept.
[thr,sorh,keepapp] = ddencmp('den','wv',x);
thr 
%%
% De-noise image using global thresholding option.
xd = wdencmp('gbl',x,'sym4',2,thr,sorh,keepapp);
figure('Color','white')
colormap(pink(255)), sm = size(map,1);
image(wcodemat(X,sm)), title('Original Image')
figure('Color','white')
colormap(pink(255))
image(wcodemat(x,sm)), title('Noisy Image')
%%
image(wcodemat(xd,sm)), title('De-Noised Image')

%% Summary
% denoising methods based on wavelet decomposition is one of the most
% significant applications of wavelets. More details can be found in the
% section of the documentation entitled "More About the Thresholding
% Strategies."