# wmulden

Wavelet multivariate denoising

## Description

example

[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.

[___] = 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).

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

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 $\sqrt{2\mathrm{ln}\left(\text{length}\left(x\right)\right)}.$

• "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