Exercise 2

Generalized Fourier Improvisations on Signal De-Noising

Question 1 (theoretical warming-up).

Assume the standard nonparametric regression model $$ y_i=g(x_i)+\epsilon_i,\;i=1,\ldots n $$

Following the generalized Fourier series approach, we choose a certain orthonormal basis $\{\phi_j(x)\}$ and try to estimate the generalized Fourier coefficients $c_{gj}=\int g(x)\phi_j(x)dx$ of the unknown function $g(x)$ from the empirical generalized Fourier coefficients $c_{yj}, j=1,\ldots,n$ of the data ${\bf y}=(y_1,\ldots,y_n)^T$. Consider the following penalized estimator $\hat{\bf c}_g$ of ${\bf c}_g$: $$ \hat{\bf c}_g=\arg \min_{\bf c} \left\{\sum_{j=1}^n (c_{yj}-c_j)^2+\lambda Pen({\bf c}) \right\}, $$ where $Pen({\bf c})$ is the penalty function and $\lambda$ is the tuning parameter controlling the trade-off between goodness-of-fit and penalty loss.

Show that

  1. $l_2$-penalty $Pen({\bf c})=||{\bf c}||^2_2=\sum_{j=1}^n c_j^2$ leads to linear shrinkage estimator
  2. $l_1$-penalty $Pen({\bf c})=||{\bf c}||_1=\sum_{j=1}^n |c_j|$ yields soft thresholding estimator
  3. $l_0$-penalty $Pen({\bf c})=||{\bf c}||_0=\#\{c_j \neq 0\}$ corresponds to hard thresholding estimator

Question 2.

Generate a noisy sample of size 1024 from the Bumps function available in DJ.EX function from the package wavethresh (see Computational Notes below) with the root signal to noise ratio 5 at uniformly designed points $x_i=\frac{i-1}{n}, i=1,...,1024$ (see help(DJ.EX()) for details). Plot the original (noiseless) Bumps signal and the generated noisy one. Now we (actually you) will try to de-noise the noisy signal using wavelets.
  1. Perform the Discrete Wavelet Transform (DWT) of the noisy data and compare it with the DWT of the noiseless Bumps signal. To "feel" wavelets better try several mother wavelets from Daubechies family by taking different filter.number (see Computational Comments). Note that filter.number=1 corresponds to the Haar basis. Comment the plots of the resulting DWT coefficients.
  2. Threshold noisy wavelet coefficients using the universal threshold of Donoho and Johnstone and one of the adaptive thresholding rules (e.g., SURE, FDR or CV). Plot thresholded wavelet coefficients and compare them with those from 2.1.
  3. Reconstruct de-noised estimates by inverse DWT. Compare the two estimates and explain the differences.
  4. Compare the wavelet estimators with a smoothing spline or kernel estimator.
Try several noisy samples from the Bumps function to strengthen your conclusions.

Question 3.

The file NMR.dat contains a nuclear magnetic resonance signal of length 1024.
  1. Plot the data and make some initial conclusions.
  2. Let's start with the Fourier cosine series. Perform the Discrete Cosine Transform (DCT) of the data. Unfortunately, there is no R function for DCT but you can either look for various R packages in CRAN or simply write it by yourself (try to avoid explicit loops as much as you can otherwise your program might run slowly for such a large data set). The DCT of vector ${\bf y}$ is given by $$ c_0=\bar{\bf y},\;c_j=\frac{1}{n}\sum_{i=1}^n \sqrt{2}\cos\left(\frac{\pi j (i-1/2)}{n}\right)y_i,\;j=1,\ldots,n-1, $$ while the inverse DCT is $$ y_i=c_0+\sum_{j=1}^{n-1}\sqrt{2} \cos\left(\frac{\pi j (i-1/2)}{n}\right) c_j, \;i=1,\ldots,n $$

    Look at the DCT of the original noisy signal and try to de-noise it - truncate series, select largest coefficients (by fixing the number of non-zero coefficients or by thresholding), etc. - improvise as much as you want. Give the resulting DCT estimator and comment the results.

  3. Now let's see what wavelets can do with this data. Perform the DWT of the signal, threshold the DWT coefficients by any method you have chosen (you can try several mother wavelets, various thresholds, etc.) and compare them with those of the original signal.
  4. Compare the resulting wavelet estimator with the DCT estimator and smoothing spline. Comment the results.
Computational Note for R users.

To perform wavelet analysis you can use the R wavethresh package written by Guy Nason. In fact, there are 3 basic functions you will need in wavethresh:

  • >y.wd=wd(y,filter.number= ,...)
    performs DWT of ${\bf y}$ and creates a .wd object $y.wd$ you'll need further. The parameter filter.number selects the mother wavelets you use in wavelet decomposition. For example, filter.number=1 corresponds to the Haar mother wavelet. The more is the filter number, the more regular is the resulting mother wavelet. To plot the resulting wavelet coefficients simply do
    >plot(y.wd)
  • >ythresh.wd=threshold.wd(y.wd,policy= ,type= ,...)
    thresholds noisy wavelet cofficients according to the type ("hard" or "soft") and the policy ("universal", "cv", "fdr", "sure", etc.) you ask. You can plot the thresholded coefficients by
    >plot(ythresh.wd)
  • >y.wr=wr(ythresh.wd)
    performs the inverse DWT and derives the resulting thresholded wavelet estimate.
See help files for details of their use. More details can be found in the book Wavelet Methods in Statistics with R by Guy Nason available also in the library.