Main Content

wmulden

Wavelet multivariate denoising

    Description

    [x_den,npc,nestcov,dec_den,pca_params,den_params] = wmulden(x,level,wname,npc_app,npc_fin,tptr,sorh) returns a denoised version x_den of the input matrix x. The 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.

    example

    [___] = wmulden(x,level,wname,"mode",extmode,npc_app,npc_fin,tptr,sorh) uses the extension mode extmode for the discrete wavelet transform (DWT).

    [___] = wmulden(dec,npc_app) uses the wavelet decomposition structure dec.

    [dec,pca_params] = wmulden("estimate",dec,npc_app,npc_fin) returns the wavelet decomposition dec and the principal components estimates pca_params.

    Examples

    collapse all

    Load a multivariate signal x together with the original signals (x_orig) and true covariance matrix (covar).

    load ex4mwden

    Set the denoising method parameters.

    level = 5;
    wname = "sym4";
    tptr = "sqtwolog";
    sorh = "s";

    Set the PCA parameters to select the number of retained principal components automatically by Kaiser's rule.

    npc_app = "kais";
    npc_fin = "kais";

    Perform multivariate denoising.

    [x_den,npc,nestcov] = wmulden(x,level,wname, ...
        npc_app,npc_fin,tptr,sorh);

    Display the original, observed, and denoised signals. The first function, which is irregular, is correctly recovered while the second function, more regular, is well denoised.

    kp = 0;
    for i = 1:4 
        subplot(4,3,kp+1)
        plot(x_orig(:,i))
        ylim([-9 12])
        title(["Original Signal ",num2str(i)])
        subplot(4,3,kp+2)
        plot(x(:,i))
        ylim([-9 12])
        title(["Observed Signal ",num2str(i)])
        subplot(4,3,kp+3)
        plot(x_den(:,i)) 
        ylim([-9 12])
        title(["Denoised Signal ",num2str(i)])
        kp = kp+3;
    end

    Figure contains 12 axes objects. Axes object 1 with title Original Signal 1 contains an object of type line. Axes object 2 with title Observed Signal 1 contains an object of type line. Axes object 3 with title Denoised Signal 1 contains an object of type line. Axes object 4 with title Original Signal 2 contains an object of type line. Axes object 5 with title Observed Signal 2 contains an object of type line. Axes object 6 with title Denoised Signal 2 contains an object of type line. Axes object 7 with title Original Signal 3 contains an object of type line. Axes object 8 with title Observed Signal 3 contains an object of type line. Axes object 9 with title Denoised Signal 3 contains an object of type line. Axes object 10 with title Original Signal 4 contains an object of type line. Axes object 11 with title Observed Signal 4 contains an object of type line. Axes object 12 with title Denoised Signal 4 contains an object of type line.

    The second output argument gives the number of retained principal components for PCA for approximations and for final PCA.

    npc
    npc = 1×2
    
         2     2
    
    

    The third output argument contains the estimated noise covariance matrix using the MCD bases on the matrix of finest details.

    nestcov
    nestcov = 4×4
    
        1.0784    0.8333    0.6878    0.8141
        0.8333    1.0025    0.5275    0.6814
        0.6878    0.5275    1.0501    0.7734
        0.8141    0.6814    0.7734    1.0967
    
    

    Compare the estimated noise covariance with the true values. The estimation is satisfactory since the values are close to the true values given by covar.

    covar
    covar = 4×4
    
        1.0000    0.8000    0.6000    0.7000
        0.8000    1.0000    0.5000    0.6000
        0.6000    0.5000    1.0000    0.7000
        0.7000    0.6000    0.7000    1.0000
    
    

    Input Arguments

    collapse all

    Input data, specified as a matrix. The input matrix x contains P signals of length N stored column-wise, where N > P.

    Wavelet Decomposition Parameters

    Level of wavelet decomposition, specified as a positive integer.

    Wavelet, specified as a character vector or string scalar. wname can specify an orthogonal or biorthogonal wavelet. For a list of supported wavelets, see wfilters.

    Data Types: char | string

    Wavelet decomposition structure of signals to be denoised, specified as a structure. dec is assumed to be the output of mdwtdec. dec can be replaced with x, wname, and level.

    Example: dec = mdwtdec("c",x,level,wname)

    Extension mode used when performing the DWT, specified as one of the following:

    mode

    DWT Extension Mode

    'zpd"

    Zero extension

    "sp0"

    Smooth extension of order 0

    "spd" (or "sp1")

    Smooth extension of order 1

    "sym" or "symh"

    Symmetric extension (half point): boundary value symmetric replication

    "symw"

    Symmetric extension (whole point): boundary value symmetric replication

    "asym" or "asymh"

    Antisymmetric extension (half point): boundary value antisymmetric replication

    "asymw"

    Antisymmetric extension (whole point): boundary value antisymmetric replication

    "ppd"

    Periodized extension (1)

    "per"

    Periodized extension (2)

    If the signal length is odd, wextend adds to the right an extra sample that is equal to the last value, and performs the extension using the "ppd" mode. Otherwise, "per" reduces to "ppd".

    The global variable managed by dwtmode specifies the default extension mode.

    Principal Components Parameters

    Principal components selection method for approximations at level level, specified as one of these.

    • If npc_app is an integer, npc_app sets the number of retained principal components for approximations at level level in the wavelet domain. npc_app must satisfy 0 ≤ npc_appP, where P is the number of columns in x.

    • If npc_app is "kais" or "heur", the wmulden function selects the number of retained principal components using Kaiser's rule or the heuristic rule automatically.

      • Kaiser's rule keeps the components associated with eigenvalues greater than the mean of all eigenvalues.

      • The heuristic rule keeps the components associated with eigenvalues greater than 0.05 times the sum of all eigenvalues.

    • Setting npc_app is "none" is equivalent to setting npc_app equal to P, where P is the number of columns in x.

    Final PCA selection method after wavelet reconstruction, specified as one of these.

    • If npc_fin is an integer, npc_fin sets the number of retained principal components for final PCA after wavelet reconstruction. npc_fin must satisfy 0 ≤ npc_finP, where P is the number of columns in x.

    • If npc_fin is "kais" or "heur", the wmulden function selects the number for final PCA using Kaiser's rule or the heuristic rule automatically.

      • Kaiser's rule keeps the components associated with eigenvalues greater than the mean of all eigenvalues.

      • The heuristic rule keeps the components associated with eigenvalues greater than 0.05 times the sum of all eigenvalues.

    • Setting npc_fin is "none" is equivalent to setting npc_fin equal to P, where P is the number of columns in x.

    Denoising Parameters

    Threshold selection rule to apply to the wavelet decomposition of x.

    • "rigsure" — Use the principle of Stein's Unbiased Risk.

    • "heursure" — Use a heuristic variant of Stein's Unbiased Risk.

    • "sqtwolog" — Use the universal threshold 2ln(length(x)).

    • "minimaxi" — Use minimax thresholding. (See thselect for more information.)

    • "penalhi", "penalme", "penallo" — Use Birgé-Massart strategy. For more information, see wthrmngr.

    Type of thresholding to perform:

    • "s" — Soft thresholding

    • "h" — Hard thresholding

    Output Arguments

    collapse all

    Denoised data, returned as a matrix.

    Selected numbers of retained principal components, returned as a vector.

    Estimated noise covariance matrix obtained using the minimum covariance determinant (MCD) estimator.

    Wavelet decomposition of the denoised data, returned as a structure with the following fields:

    • dirDec — Direction indicator: 'r' (row) or 'c' (column)

    • level — Level of wavelet decomposition

    • wname — Wavelet name

    • dwtFilters — Structure with four fields: LoD, HiD, LoR, and HiR

    • dwtEXTM — DWT extension mode

    • dwtShift — DWT shift parameter (0 or 1)

    • dataSize — Size of x

    • ca — Approximation coefficients at level level

    • cd — Cell array of detail coefficients, from level 1 to level level

    The coefficients ca and cd{k}, for k from 1 to level, are matrices and are stored in rows if dirdec = 'r' or in columns if dirdec = 'c'.

    Principal components estimates, returned as a structure such that:

    pca_params.NEST = {pc_NEST,var_NEST,NESTCOV}
    pca_params.APP  = {pc_APP,var_APP,npc_APP}
    pca_params.FIN  = {pc_FIN,var_FIN,npc_FIN}
    
    where

    • pc_XXX is a P-by-P matrix of principal components.

      The columns are stored in descending order of the variances.

    • var_XXX is the principal component variances vector.

    • NESTCOV is the covariance matrix estimate for detail at level 1.

    Denoising parameters, returned as a structure.

    • den_params.thrVAL is a vector of length level which contains the threshold values for each level.

    • den_params.thrMETH is a character vector containing the name of the denoising method (tptr).

    • den_params.thrTYPE is a character variable containing the type of the thresholding (sorh).

    Algorithms

    The multivariate denoising procedure is a generalization of the one-dimensional strategy. It combines univariate wavelet denoising in the basis where the estimated noise covariance matrix is diagonal and non-centered Principal Component Analysis (PCA) on approximations in the wavelet domain or with final PCA.

    The robust estimate of the noise covariance matrix given by the minimum covariance determinant estimator based on the matrix of finest details.

    References

    [1] Aminghafari, Mina, Nathalie Cheze, and Jean-Michel Poggi. “Multivariate Denoising Using Wavelets and Principal Component Analysis.” Computational Statistics & Data Analysis 50, no. 9 (May 2006): 2381–98. https://doi.org/10.1016/j.csda.2004.12.010.

    [2] Rousseeuw, Peter J., and Katrien Van Driessen. “A Fast Algorithm for the Minimum Covariance Determinant Estimator.” Technometrics 41, no. 3 (August 1999): 212–23. https://doi.org/10.1080/00401706.1999.10485670.

    Version History

    Introduced in R2006b