General Solution
We have a Deconvolution problem with known operator.
One way to define the objective function is:
$$ \arg \min_{x} \frac{1}{2} {\left\| h \ast x - y \right\|}_{2}^{2} = \arg \min_{x} f \left( x \right) $$
The are 3 method to solve this:
- Use Gradient Descent by deriving the gradient of $ f \left( \cdot \right) $.
- Write the problem in Matirx Form and either use direct solver (One could use Toeplitz System Solver) of this Least Squares problem or again use iterative solver in Matrix Form.
- Solve this in Fourier Frequency Domain.
I previously solved similar problems using the Matrix Form so this time we'll solve it in the convolution form.
The Gradient
Clearly the derivative of Inner Product is related to the adjoint of the operator.
In the case of Convolution the Adjopint is the Correlation Operator.
Let's define the Convolution Operation $ \ast $ with MATLAB Code.
So $ y = h \ast x $ will be in MATLAB:
vY = conv(vX, vH, 'valid');
Defining $ f \left( \cdot \right) $ in MATLAB would be:
hObjFun = @(vX) 0.5 * sum((conv(vX, vH, 'valid') - vY) .^ 2);
With the Derivative being:
vG = conv2((conv2(vX, vH, 'valid') - vY), vH(end:-1:1), 'full');
Where when we convolve with the flipped version of a vector it means we're basically doing correlation.
Solution in MATLAB
So, one must pay attention that since we use valid
convolution the size of $ y $ and $ x $ are different (Depends on the size of $ h $).
Usually we chose $ h $ to have odd number of elements to have its radius to be defined.
Then the MATLAB solution would be:
vX = zeros(size(vY, 1) + (2 * kernelRadius), 1);
vObjVal(1) = hObjFun(vX);
for ii = 1:numIterations
vG = conv2((conv2(vX, vH, 'valid') - vY), vH(end:-1:1), 'full');
vX = vX - (stepSize * vG);
vObjVal(ii + 1) = hObjFun(vX);
end
Specific Solution
In your case you asked for $ h $ to be Box Blur.
So I did the above with this model. Let's go through results.
First, let's see if the solution we got to is a real solution.
In order to examine that, we will convolve the result of the optimization with $ h $ and compare result to the input data. We expect it to be similar to the input.
As can be seen above it is, indeed, a perfect reconstruction of the input data.
We can see this by the objective value as well:
Yet the result is not good:
So, why is the result so bad?
There are 2 options:
- The SNR is not good enough to solve the inverse problem as is. It requires additional regularization (Wiener / Tikhonov (They are the same)).
- The model of Box Blur doesn't fit.
So in my code I implemented Wiener / Tikhonov Regularization (See paramLambda
in my code).
The full code is available on my StackExchange Signal Processing Q51460 GitHub Repository.