Main Content

wpbmpen

Penalized threshold for wavelet packet denoising

collapse all in page

    Syntax

    thr = wpbmpen(t,sigma,alpha)
    thr = wpbmpen(t,sigma,alpha,ARG)

    Description

    thr = wpbmpen(t,sigma,alpha) returns a global threshold for denoising. t is a wavelet packet tree corresponding to the wavelet packet decomposition of the signal or image you want to denoise. sigma and alpha are parameters in a penalization method used to obtain the threshold. For more information, see Penalized Criterion.

    example

    thr = wpbmpen(t,sigma,alpha,ARG), where ARG is arbitrary, also plots three curves (see Penalized Criterion):

    • 2×sigma^2×t×(alpha + log(n/t))

    • sum(c(k)^2,k≤t)

    • crit(t)

    Examples

    collapse all

    Open Live Script

    Denoise Signal

    Load a noisy chirp signal.

    load noischir
x = noischir;

    Obtain the wavelet packet decomposition of the signal at level 5 using the sym6 wavelet.

    wname = "sym6";
lev = 5;
t = wpdec(x,lev,wname);

    Estimate the noise standard deviation from the detail coefficients at level 1, corresponding to node index 2.

    det1 = wpcoef(t,2);
sigma = median(abs(det1))/0.675;

    Use wpbmpen to select the global thresholding for denoising, using the recommended value of alpha.

    alpha = 2;
thr = wpbmpen(t,sigma,alpha)
    thr = 
4.5740

    Use wpdencmp to denoise the signal using the threshold. Use soft thresholding and keep the approximation.

    keepapp = 1;
xd = wpdencmp(t,"s","nobest",thr,keepapp);

    Plot the original and denoised signals.

    tiledlayout(2,1)
nexttile
plot(x)
title("Original Signal")
nexttile
plot(xd)
title("Denoised Signal")

    Figure contains 2 axes objects. Axes object 1 with title Original Signal contains an object of type line. Axes object 2 with title Denoised Signal contains an object of type line.

    Denoise Image

    Load an image.

    load flower
img = flower.Noisy;

    Obtain the wavelet packet decomposition of the image at level 3 using the coif2 wavelet.

    wname = "coif2";
lev = 3;
t = wpdec2(img,lev,wname);

    Estimate the noise standard deviation from the detail coefficients at level 1.

    det1 = [wpcoef(t,2) wpcoef(t,3) wpcoef(t,4)];
sigma = median(abs(det1(:)))/0.6745;

    Use wpbmpen to select the global thresholding for denoising.

    alpha = 1.1;
thr = wpbmpen(t,sigma,alpha)
    thr = 
0.2050

    Use wpdencmp to denoise the image using the threshold. Use soft thresholding and keep the approximation.

    keepapp = 1;
xd2 = wpdencmp(t,"s","nobest",thr,keepapp);

    Plot the original and denoised images.

    figure
colormap(gray)
tiledlayout(1,2)
nexttile
imagesc(img)
title("Original Image")
nexttile
imagesc(xd2)
title("Denoised Image")

    Figure contains 2 axes objects. Axes object 1 with title Original Image contains an object of type image. Axes object 2 with title Denoised Image contains an object of type image.

    Input Arguments

    collapse all

    Wavelet packet tree corresponding to the wavelet packet decomposition of the signal or image to be denoised, specified as a wptree object.

    Example: t = wpdec(x,5,"sym4")

    Standard deviation of the zero mean Gaussian white noise in the denoising model, specified as a scalar. For more information, see wnoisest and Penalized Criterion.

    Data Types: double

    Tuning parameter for the penalty term, specified as a scalar greater than 1. The sparsity of the wavelet representation of the denoised signal or image grows with alpha. Typically, alpha equals 2. For more information, see Penalized Criterion.

    Data Types: double

    Output Arguments

    collapse all

    Global threshold for denoising, returned as a scalar. thr is obtained by a wavelet coefficients selection rule using a penalization method provided by Birgé-Massart. For more information, see Penalized Criterion.

    More About

    collapse all

    Penalized Criterion

    The global threshold thr is obtained by a wavelet coefficients selection rule using a penalization method provided by Birgé-Massart.

    The global threshold minimizes the penalized criterion given by the following:

    Let t* be the minimizer of

    crit(t) = -sum(c(k)^2,k≤t) + 2×sigma^2×t×(alpha + log(n/t)),
    where c(k) are the wavelet coefficients sorted in decreasing order of their absolute value, and n is the number of coefficients; then thr = |c(t*)|.

    References

    [1] Birgé, Lucien, and Pascal Massart. “From Model Selection to Adaptive Estimation.” In Festschrift for Lucien Le Cam, edited by David Pollard, Erik Torgersen, and Grace L. Yang, 55–87. New York, NY: Springer New York, 1997. https://doi.org/10.1007/978-1-4612-1880-7_4.

    Version History

    Introduced before R2006a

    See Also

    | | |