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

    %% Multivariate Wavelet Denoising
% The purpose of this example is to show the features of multivariate
% denoising provided in Wavelet Toolbox(TM).
% 
% Multivariate wavelet denoising problems deal with models of the form
%
% $$X(t) = F(t) + e(t)$$
%
% where the observation _*X*_ is _*p*_-dimensional, _*F*_ is
% the deterministic signal to be recovered, and _*e*_ is a
% spatially-correlated noise signal. This example uses a number of noise
% signals and performs the following steps to denoise the deterministic
% signal.

% Copyright 2007-2012 The MathWorks, Inc.

%% Loading a Multivariate Signal
% To load the multivariate signal, type the following code at the MATLAB(R)
% prompt:
 load ex4mwden
 whos
%%
% Usually, only the matrix of data |x| is available. Here, we also have the
% true noise covariance matrix |covar| and the original signals |x_orig|. 
% These signals are noisy versions of simple combinations of the two
% original signals. The first signal is "Blocks" which is irregular, and
% the second one is "HeavySine" which is regular, except around time 750.
% The other two signals are the sum and the difference of the two original
% signals, respectively. Multivariate Gaussian white noise exhibiting
% strong spatial correlation is added to the resulting four signals, which
% produces the observed data stored in |x|.

%% Displaying the Original and Observed Signals
% To display the original and observed signals, type:
kp = 0;
for i = 1:4 
    subplot(4,2,kp+1), plot(x_orig(:,i)); axis tight;
    title(['Original signal ',num2str(i)])
    subplot(4,2,kp+2), plot(x(:,i)); axis tight;
    title(['Observed signal ',num2str(i)])
    kp = kp + 2;
end
%%
% The true noise covariance matrix is given by: 
covar

%% Removing Noise by Simple Multivariate Thresholding
% The denoising strategy combines univariate wavelet denoising in the
% basis, where the estimated noise covariance matrix is diagonal with
% noncentered Principal Component Analysis (PCA) on approximations in the
% wavelet domain or with final PCA. 
% 
% First, perform univariate denoising by typing the following lines to set 
% the denoising parameters:
level = 5;
wname = 'sym4';
tptr  = 'sqtwolog';
sorh  = 's';
%%
% Then, set the PCA parameters by retaining all the principal components:
npc_app = 4;
npc_fin = 4;
%%
% Finally, perform multivariate denoising by typing: 
x_den = wmulden(x, level, wname, npc_app, npc_fin, tptr, sorh);

%% Displaying the Original and Denoised Signals
% To display the original and denoised signals type the following:
clf
kp = 0;
for i = 1:4 
    subplot(4,3,kp+1), plot(x_orig(:,i)); axis tight; 
    title(['Original signal ',num2str(i)])
    subplot(4,3,kp+2), plot(x(:,i)); axis tight;
    title(['Observed signal ',num2str(i)])
    subplot(4,3,kp+3), plot(x_den(:,i)); axis tight;
    title(['Denoised signal ',num2str(i)])
    kp = kp + 3;
end

%% Improving the First Result by Retaining Fewer Principal Components
% We can see that, overall, the results are satisfactory. Focusing on the
% two first signals, note that they are correctly recovered, but we can
% improve the result by taking advantage of the relationships between the
% signals, leading to an additional denoising effect. 
% 
% To automatically select the numbers of retained principal components
% using Kaiser's rule, which retains components associated with eigenvalues 
% exceeding the mean of all eigenvalues, type: 
npc_app = 'kais';
npc_fin = 'kais';
%%
% Perform multivariate denoising again by typing:
[x_den, npc, nestco] = wmulden(x, level, wname, npc_app, ...
     npc_fin, tptr, sorh);

%% Displaying the Number of Retained Principal Components
% The second output argument |npc| is the number of retained principal
% components for PCA for approximations and for final PCA.
npc
%%
% As expected, because the signals are combinations of two original
% signals, Kaiser's rule automatically detects that only two principal
% components are of interest. 
% 

%% Displaying the Estimated Noise Covariance Matrix
% The third output argument |nestco| contains the estimated noise
% covariance matrix: 
nestco
%%
% As it can be seen by comparing it with the true matrix covar given
% previously, the estimation is satisfactory. 

%% Displaying the Original and Final Denoised Signals
% To display the original and final denoised signals type:
kp = 0;
for i = 1:4 
    subplot(4,3,kp+1), plot(x_orig(:,i)); axis tight;
    title(['Original signal ',num2str(i)])
    subplot(4,3,kp+2), plot(x(:,i)); axis tight;
    title(['Observed signal ',num2str(i)])
    subplot(4,3,kp+3), plot(x_den(:,i)); axis tight;
    title(['Denoised signal ',num2str(i)])
    kp = kp + 3;
end
%%
% These results are better than those previously obtained. The first
% signal, which is irregular, is still correctly recovered, while the
% second signal, which is more regular, is better denoised after this
% second stage of PCA.

%% Learning More About Multivariate Denoising
% You can find more information about multivariate denoising, including
% some theory, simulations, and real examples, in the following reference:
% 
% M. Aminghafari, N. Cheze and J-M. Poggi (2006), "Multivariate denoising
% using wavelets and principal component analysis," Computational
% Statistics & Data Analysis, 50, pp. 2381-2398.