Least-squares Signal Declipping

Goal: implement the least-squares signal declipping method.

Introduction

Clipping occurs when a detector (or an analog-to-digital converter, or an amplifier) attempts to process a signal that exceeds its dynamic range. This results in the signal being "flattened" or "clipped" at the maximum and/or minimum representable values of the detector. This distortion degrades the signal quality.

Declipping is the process of trying to reconstruct the original, unclipped signal from its clipped version.

Least-squares Declipping

Let us say that we are given a clipped signal, the size-`N` vector `y`, where `n` elements (with indices `m_1, m_2, ..., m_n`) are clipped (where `n < N`), that is, replaced with `y_max` if they exceeded `y_max` (or with `y_min` if they fell below `y_min`).

We need to find the size-`n` vector `z` of the reconstructed values that were clipped such that the total reconstructed vector `x` is given as:

• `x_i = y_i`, if `i != m_k` for the unclipped elements, and
• `x_(m_k) = z_k`, for `k=1, ..., n` for the reconstructed elements.

That is, the reconstructed vector `x` is searched in the form

`x = tilde(y) + Mz`, where `tilde(y)` is the clipped signal where the clipped elements (that is, the elements equal to `y_max` or `y_min`) are replaced with zeros. The matrix `M` inserts the elements of the vector `z` into the clipped places. Specifically, matrix `M` is an `N times n` matrix with all elements equal zero except for `M_(m_k,k) = 1`. The vector `z` of the reconstructed elements is then obtained by minimizing the square form `L(z) = ||D(tilde(y) + Mz)||^2`, where `D` is the third-order finite difference matrix (see the lecture note). This form minimizes the third derivative, forcing the declipped signal `z` to be in the form of a parabola.

You should solve this least-squares minimization problem using your own implementation of the QR-factorization method as described in the lecture note.

Task

1. Implement a function with the signature

   vector declip(vector y, double y_min, double y_max)

   that takes a clipped signal y and the clipping limits y_min and y_max returns the reconstructed signal x.
2. Test your implementation on a synthetic clipped signal (e.g., a sine wave clipped at `y_max=0.8` and `y_min=-0.8`); plot the original, clipped, and declipped signals for comparison.
3. Apply your implementation on a more complicated clipped signal (at your imagination).