# wmulden

Wavelet multivariate denoising

## Syntax

``````[x_den,npc,nestcov,dec_den,pca_params,den_params] = wmulden(x,level,wname,npc_app,npc_fin,tptr,sorh)``````
``````[___] = wmulden(x,level,wname,"mode",extmode,npc_app,npc_fin,tptr,sorh)``````
``````[___] = wmulden(dec,npc_app)``````
``````[dec,pca_params] = wmulden("estimate",dec,npc_app,npc_fin)``````

## 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`).

`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``` 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_app ≤ P```, 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_fin ≤ P```, 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)}.$

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

 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.

 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.