# Use MATLAB to Restore a Signal from a Given Degraded Signal Using Tikhonov Regularization

Anyone could share how to develop an application in MATLAB to restore the signal from a given degraded signal using Tikhonov regularization i.e restoring the signal $f$ via solving

$$\min || g - f \circ h ||^2 + \lambda ||\sigma \circ f ||^2$$

where $\sigma = [1,-1]$. and $\circ$ denotes the circulant convolution operator. The setting of the value epsilon is dependent on the signal to noise ratio.

Given an audio file with noisy convoluted measurement,

The following routine are expected Tikdeconv.m

Tikdeconv De-convolue a signal using Tikhonov regularisation J=deblurw(I,PSF,LAMBDA) de-convolute input signal $I$ by PSF using Tikhonov regularisation with regularisation parameter LAMBDA returning de-convoluted signal J. The assumption is that the signal $I$ was created by convolving a signal $J$ with kernel PSF and possible by adding /guassian white noise. You may not call or use any code from the following built-in function in MATLAB:deconvwr,deconv,deconv2

Any feedback on how to solve this would be appreciated. Thank you

• the information of the low pass convolution kernel is also given Oct 15, 2015 at 2:34
• MA4272 :) Good luck with your assignment... Oct 15, 2015 at 4:28
• Hi, I don't understand the term $\sigma = \left[ 1, -1 \right]$, do you mean the derivative of the signal? We're talking about 1D signal, right?
– Royi
Oct 17, 2015 at 15:45
• @Drazick σ=[1,−1] is the gradient and yes its 1D signal Oct 18, 2015 at 15:03
• I see. I think I know how to solve it. Hopefully tomorrow I will have time to write it down.
– Royi
Oct 18, 2015 at 18:04

The idea is to represent all operation sing Matrices.
Once it is done, it is easy to solve the problems as a Least Squares problems.

The way to represent Convolution Operation using a Matrix is by Toeplitz Matrix.
For 1D it is pretty straight forward to do (Just pay attention to boundary).

So let's take the simple model in the comment:

$$g = f \circ h + n \rightarrow g = H f + n$$

Now, the LS solution is easy:

$$\arg \min_{f} {\left \| H f - g \right \|}^{2}_{2} = { \left( {H}^{T} H \right) }^{-1} {H}^{T} g$$

The matrix form for the regularized model is:

$$\arg \min_{f} {\left \| H f - g \right \|}^{2}_{2} + \lambda {\left \| \Sigma f \right \|}^{2}_{2}$$

Where $$H$$ is the matrix Toeplitz Matrix generated from the Blur Kernel and $$\Sigma$$ is the Toeplitz matrix generated form of the derivative operator (In this case, keeps the signal smooth to regulate the High Pass property of the first term).

This is actually Tikhonov Regularization and it is easily solved.

$$\arg \min_{f} {\left \| H f - g \right \|}^{2}_{2} + \lambda {\left \| \Sigma f \right \|}^{2}_{2} = { \left( {H}^{T} H + \lambda {\Sigma}^{T} \Sigma \right) }^{-1} {H}^{T} g$$

All the above is pretty straight forward to calculate in MATLAB once the the $$H$$ and $$\Sigma$$ matrices are formed (Given the Convolution Kernel, the function convmtx can assist you generate the matrix).

The way to create those matrices is pretty simple. Given a signal of the length $$N$$ (Samples) implies those matrices are $$N \times N$$ (Assuming the output has the same number of samples).
Where the $$i$$ row of the matrix applies the linear combination for the $$i$$ output sample.

For instance, given the Gradient Kernel and a signal of length 5 the matrix becomes:

$$\begin{bmatrix} 1 & -1 & 0 & 0 & 0 \\ 0 & 1 & -1 & 0 & 0 \\ 0 & 0 & 1 & -1 & 0 \\ 0 & 0 & 0 & 1 & -1 \\ 0 & 0 & 0 & 0 & 1 \end{bmatrix}$$

One must pay attention to the rows on the beginning and the end to apply the required boundary conditions.
In the case above I chose zero boundary conditions (Last row).

One should note, that in practice, the convolution matrix is not feasible due to memory limitations.
Hence it must be represented by a Sparse Matrix (Since it is mostly zeros) or the problem should be solved using a solver which doesn't require the matrix explicitly (Like Conjugate Gradient).

Another note would be that if Circular Convolution is after (Imitating DFT Based Convolution) the matrices should not only be Toeplitz but Circulant Matrices which means they enforce periodic boundary conditions (See the circular property in MATLAB padarray function).

Here is a MATLAB Code for the exercise.
I only used 800 samples due to memory constraints.

% Audio Deconvolution
% See http://dsp.stackexchange.com/questions/26433

%% General Parameters and Initialization

% The follwoing code must be executed before any of the following cells.
% The cells must be executed one after one in the order they are written.
clear();
close('all');

set(0, 'DefaultFigureWindowStyle', 'docked');
defaultLooseInset = get(0, 'DefaultAxesLooseInset');
set(0, 'DefaultAxesLooseInset', [0.05, 0.05, 0.05, 0.05]);

titleFontSize   = 14;
axisFotnSize    = 12;
stringFontSize  = 12;

thinLineWidth   = 2;
normalLineWidth = 3;
thickLineWidth  = 4;

smallSizeData   = 36;
mediumSizeData  = 48;
bigSizeData     = 60;

randomNumberStream = RandStream('mlfg6331_64', 'NormalTransform', 'Ziggurat');
subStreamNumber = 57162;
set(randomNumberStream, 'Substream', subStreamNumber);
RandStream.setGlobalStream(randomNumberStream);

vInputSignal    = input;
vRefSignal      = handel;

vBlurKernel     = [0.0545; 0.2442; 0.4026; 0.2442; 0.0545];
vGradientKernel = [-1; 1]; % Forward Finite Differences
% vGradientKernel = [1; -1]; % Backward Finite Differences

%% Parameters
% Limiting the Number of Samples due to Memory Constraints
numSamplesToProcess = 800;

numSamples = size(vInputSignal, 1);

numSamples = min([numSamplesToProcess, numSamples]);

vInputSignal    = vInputSignal(1:numSamples);
vRefSignal      = vRefSignal(1:numSamples);

% Creating the Operators in Matrix Form
mBlurKernel     = convmtx(vBlurKernel, numSamples);

% Enforcing size to match Convolution using 'same' property
mBlurKernel     = mBlurKernel(1:numSamples, :);

% Enforcing Circulant Matrix (Like DFT Based Convolution)
mBlurKernel(1, (end - 3):end) = vBlurKernel(1:4);
mBlurKernel(2, (end - 2):end) = vBlurKernel(1:3);
mBlurKernel(3, (end - 1):end) = vBlurKernel(1:2);
mBlurKernel(4, (end - 0):end) = vBlurKernel(1:1);

regFactor = 0.02; % Regularization Factor - Lambda

%% Least Squares Solution

vLsRestoredSignal = (mBlurKernel.' * mBlurKernel) \ (mBlurKernel.' * vInputSignal);

%% Regularized (Tikhonov) Least Squares Solution

vRegLsRestoredSignal = ((mBlurKernel.' * mBlurKernel) + (regFactor * (mGradientKernel.' * mGradientKernel))) \ (mBlurKernel.' * vInputSignal);

%% Display Results

hFigure = figure();
hAxes   = axes();
hLineSeries = plot(1:numSamples, [vRefSignal, vInputSignal]);
set(hLineSeries, 'LineWidth', normalLineWidth);
set(get(hAxes, 'Title'), 'String', ['Reference and Input Signal'], ...
'FontSize', titleFontSize);
set(get(hAxes, 'XLabel'), 'String', 'Samples Index', ...
'FontSize', axisFotnSize);
set(get(hAxes, 'YLabel'), 'String', 'Samples Value', ...
'FontSize', axisFotnSize);
hLegend = legend({['Reference Signal'], ['Input Signal']});

hFigure = figure();
hAxes   = axes();
hLineSeries = plot(1:numSamples, [vRefSignal, vLsRestoredSignal, vRegLsRestoredSignal]);
set(hLineSeries, 'LineWidth', normalLineWidth);
set(get(hAxes, 'Title'), 'String', ['Restored Signals and the Reference Signal'], ...
'FontSize', titleFontSize);
set(get(hAxes, 'XLabel'), 'String', 'Samples Index', ...
'FontSize', axisFotnSize);
set(get(hAxes, 'YLabel'), 'String', 'Samples Value', ...
'FontSize', axisFotnSize);
hLegend = legend({['Reference Signal'], ['Least Squares Restored Signal'], ['Regularized Least Squares Restored Signal']});

hFigure = figure();
hAxes   = axes();
hLineSeries = plot(1:numSamples, [vInputSignal, vLsRestoredSignal, vRegLsRestoredSignal]);
set(hLineSeries, 'LineWidth', normalLineWidth);
set(get(hAxes, 'Title'), 'String', ['Restored Signals and the Input Signal'], ...
'FontSize', titleFontSize);
set(get(hAxes, 'XLabel'), 'String', 'Samples Index', ...
'FontSize', axisFotnSize);
set(get(hAxes, 'YLabel'), 'String', 'Samples Value', ...
'FontSize', axisFotnSize);
hLegend = legend({['Input Signal'], ['Least Squares Restored Signal'], ['Regularized Least Squares Restored Signal']});

%% Restore Defaults
set(0, 'DefaultFigureWindowStyle', 'normal');
set(0, 'DefaultAxesLooseInset', defaultLooseInset);


The code is available at GitHub - GitHub StackExchangeCodes (Look at the SignalProcessing\26433 folder).

• thank you very much! but could you explain further in particular on how to form H and Σ matrices? any references or examples would really help! Oct 19, 2015 at 8:43
• Hi Richard, Please mark the question as answered. I will add about the construction in few hours. You can see here for instance - ecsutton.ece.ufl.edu/dip/handouts/chapt05p2.pdf (For images) or in the Wikipedia link I added. Anyhow, I will give an example later.
– Royi
Oct 19, 2015 at 8:49
• This could also be helpful - dsp.stackexchange.com/questions/15096/….
– Royi
Oct 19, 2015 at 8:54
• The problem stipulates circular convolution, so the matrix H should be circulant (not Toeplitz). Oct 23, 2015 at 16:34
• @LukeP, You are right, didn't notice. I updated the code accordingly.
– Royi
Oct 24, 2015 at 8:12