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);
addpath(genpath('RawData'));
addpath(genpath('AuxiliaryFunctions'));
%% Loading Data
load('input.mat');
load('handel.mat');
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);
mGradientKernel = convmtx(vGradientKernel, numSamples);
% Enforcing size to match Convolution using 'same' property
mBlurKernel = mBlurKernel(1:numSamples, :);
mGradientKernel = mGradientKernel(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);
mGradientKernel(1, end:end) = vGradientKernel(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).
