Multidimensional signal restoration

In multidimensional signal processing, Multidimensional signal restoration refers to the problem of estimating the original input signal from observations of the distorted or noise contaminated version of the original signal using some prior information about the input signal and /or the distortion process. Multidimensional signal processing systems such as audio, image and video processing systems often receive as input, signals that undergo distortions like blurring, band-limiting etc. during signal acquisition or transmission and it may be vital to recover the original signal for further filtering. Multidimensional signal restoration is an inverse problem, where only the distorted signal is observed and some information about the distortion process and/or input signal properties is known. A general class of iterative methods have been developed for the multidimensional restoration problem with successful applications to multidimensional deconvolution, signal extrapolation  and denoising.

Definition
In general, the multidimensional signal restoration problem can be represented by an equation of the form,


 * $$ y(n_1, n_2,....,n_m) = D[ x(n_1, n_2, .....,n_m)]$$

where $$y(n_1, n_2,....,n_m)$$ represents the observed m-dimensional distorted output signal, $$x(n_1, n_2,....,n_m)$$ represents the m-dimensional undistorted input signal and $$D[\cdot]$$ represents the distortion operator acting upon the input signal. $$D[\cdot]$$ can be used to model a wide range of transformations such as blurring, additive noise, time limiting, band limiting etc. of multidimensional signals.

A simple straightforward solution to above equation is of the form,


 * $$ x(n_1, n_2,....,n_m) = D^{-1}[ y(n_1, n_2, .....,n_m)]$$

where $$D^{-1}[\cdot]$$ is the inverse distortion operator.

However, in most cases of practical use, it may be extremely difficult to implement the inverse distortion operator $$D^{-1}[\cdot]$$ or the such an inverse distortion operator may not exist and even in situations where the distortion operator $$D[\cdot]$$ is known and its inverse can be approximately implemented, the resultant reconstructed signal $$\tilde x(n_1, n_2,....,n_m)$$ can have very large reconstruction errors due to the inaccuracies present in the estimation of the inverse operator $$D^{-1}[\cdot]$$. A general class of iterative methods based on the idea of successive approximation is used to estimate the unknown input signal $$x(n_1, n_2,....,n_m)$$.

Generalized constrained iIterative signal restoration
Since a simple approach of recovering the input signal $$x(n_1, n_2,....,n_m)$$ by implementing the inverse distortion operator $$D^{-1}[\cdot]$$ on the observed signal $$y(n_1, n_2,....,n_m)$$ is not practical, specific iterative restoration algorithms were developed for certain types of distortions like blurring, finite signal domain support, finite frequency domain support of signals etc. making certain assumptions about the properties of the input signal such as finite time/spatial-domain support, non-negativity etc. A generalized iterative method that can model the above-mentioned distortions and signal domain specific constraints was later developed.

The general iterative solution based on successive approximation can have the following form,


 * $$ \tilde x_{k+1}(n_1, n_2,....,n_m) = F[ \tilde x_k(n_1, n_2, .....,n_m)]$$

where $$ \tilde x_{k+1}(n_1, n_2,....,n_m)$$ is the estimate of the input signal at iteration $$ k+1$$, $$ \tilde x_{k}(n_1, n_2,....,n_m)$$ is the estimate at iteration $$ k$$ and $$ F[\cdot]$$ represents the iteration operator that relates signal estimate at iteration $$k+1$$ to the signal estimate at iteration $$k$$.

In many cases, certain signal domain properties of the input signal to be reconstructed are known and can usually be modelled as a constraint. The constraint can be defined by the constraint operator $$C$$, such that


 * $$ x(n_1, n_2,....,n_m) = Cx(n_1, n_2,....,n_m) $$

only if $$ x(n_1, n_2,....,n_m) $$ satisfies the constraint $$C$$. It has been shown that such a constraint operator formulation can be used to model signal domain properties like non-negativity, finite frequency domain support, finite spatial domain support. The observed signal $$y(n_1, n_2,....,n_m)$$ can thus be represented in terms of the distortion operator $$D$$ and signal domain constraint $$C$$ as


 * $$ y(n_1, n_2,....,n_m) = DC x(n_1, n_2,....,n_m)$$

where the concatenation $$DC$$ represents the sequence of enforcing a signal domain constraint followed by a distortion operation on the input signal $$ x(n_1, n_2,....,n_m) $$. Under the assumption that the conditions for uniqueness and convergence of the iterative solution are met, the generalized constrained iterative restoration solution is given as


 * $$ \tilde x_{k+1}(n_1, n_2,....,n_m) = \tilde x_0(n_1,n_2,...,n_m) + (I-\lambda D)Cx_{k}(n_1,n_2,....,n_m) $$


 * $$ \tilde x_0(n_1,n_2,....,n_m) = \lambda y $$

where $$\lambda$$ is a constant to control convergence rate, $$I$$ is the identity matrix and $$\tilde x_0(n_1,n_2,...,n_m)$$ is the initial estimate of $$ x(n_1, n_2,....,n_m) $$.

Constrained iterative deconvolution
In cases where the distortion operator is both linear and shift invariant, the distortion of the input signal can be easily modelled as a convolution


 * $$ y(n_1, n_2,....,n_m) = x(n_1, n_2, .....,n_m) \ast \ast h(n_1,n_2,..,n_m) $$

where $$h(n_1,n_2,..,n_m) $$ represents the impulse response of the linear shift-invariant distortion filter. Under the assumption of linear shift-invariance, the general signal restoration problem can be transformed into a deconvolution problem with the following easily implementable iterative solution,


 * $$ \tilde x_{k+1}(n_1,n_2,...,n_m) = \lambda y(n_1,n_2,....,n_m) + g(n_1,n_2,....,n_m)\ast \ast C[\tilde x_k(n_1,n_2,....,n_m)] $$


 * $$ g(n_1,n_2,....,n_m) = \delta(n_1,n_2,....,n_m) - \lambda h(n_1,n_2,....,n_m) $$

where $$ \delta(n_1,n_2,...,n_m)$$ is an m-dimensional impulse and $$ \lambda$$ is a constant for controlling the rate of convergence. Although this solution can be easily implemented by convolution, the iterations converge to a solution only when $$H(\omega_1,\omega_2,....,\omega_m) \ne 0 $$, where $$ H(\omega_1,\omega_2,....,\omega_m)$$ represents the frequency response of the distortion filter $$h(n_1,n_2,...,n_m)$$.

By imposing a signal domain constraint of finite extent support and positivity over the finite region of support, the constrained iterative deconvolution solution can be guaranteed to converge. Such a signal domain constraint can be realistically imposed for many cases of practical use. For example, in the case of image deblurring, the blur kernel can be assumed to have a positive impulse response over a finite region of support.

Signal restoration from phase
In certain multidimensional signal processing applications, the phase of the frequency domain response of the input signal may be preserved even after undergoing distortion. For phase preserving distortions, it is possible to uniquely restore a multidimensional signal entirely from the phase of its Fourier transform as long as conditions of uniqueness and convergence are met.

The idea of recovering a signal from the phase of the frequency domain response of the input signal is a particularly useful result for images( 2-D signals). Assuming a phase preserving distortion and the existence of a unique solution for recovering a signal from its phase, the phase-based signal restoration algorithm takes the form of an iterative transformation between frequency domain and signal domain, where a frequency domain constraint (phase preservation) is first enforced on the Fourier Transform of the current estimate of the signal, followed by a spatial domain constraint (finite region of support) that is enforced in the signal domain on the current estimate of the signal.

Signal restoration from magnitude
Similar to the phase-based restoration of an unknown input signal, it is also possible to restore a signal from the magnitude of the frequency domain response of the observed signal. In certain optical systems, it is much more easier to measure the magnitude of the signal or the magnitude of its Fourier transform, but it is very difficult to precisely measure the phase of either the signal or its Fourier transform. Such cases represent a magnitude preserving distortion acting on the input signal.

Assuming a magnitude preserving distortion and the existence of a unique solution for recovering a signal from its magnitude, the magnitude-based signal restoration algorithm takes the form of an iterative transformation between frequency domain and signal domain, where a frequency domain constraint (magnitude preservation) is first enforced on the Fourier Transform of the current estimate of the signal, followed by a spatial domain constraint (finite region of support) that is enforced in the signal domain on the current estimate of the signal.

Although the magnitude-based signal restoration approach is very similar to the phase-based recovery approach, the convergence of the magnitude-based recovery approach to an acceptable result is much more difficult to achieve. In general, starting with a zero phase estimate or a random phase estimate for the magnitude-based signal recovery approach may not result in convergence, where as in the case of the phase-based signal recovery approach, even starting with a constant unit magnitude for the Fourier transform of the estimated signal, results in convergence. Random or zero phase initialization for magnitude-based signal recovery of images, usually does not result in acceptable reconstruction results even after a large number of iterations. On the other hand, starting with an initial phase information that is a noisy or heavily quantized version( but not random or uniform phase) of the original phase information, results in a very quick convergence of the magnitude-based signal recovery approach. It has been shown that an image can be perfectly recovered from the magnitude of its Fourier transform and 1-bit quantization of the original signal phase information (i.e. the initial estimate for phase of the Fourier transform can have only two values, either $$0$$ or $$\pi$$).