Deconvolution of a 1D Signal with Known Kernel (Square Wave)

Question

I have a signal measured from a radiation detector in a narrow beam of radiation. The peaks I get are quasi-gaussian in shape, see attached picture.

The signal is not a function of time, rather a function of distance. The x-axis is is mm and the y-axis is in arbitrary detector response.

The detector used to measure this signal had a finite width, which is contributing to a broadening of the peaks. What I want to do is deconvolve a square wave of width equal to that of the detector from the signal, hence removing some of the broadening effect.

I am hoping to do this in matlab, however I am having trouble using the deconv function due to each data set being two vectors, x and y, and the fact that each data set is a function of linear distance not time.

Any ideas on how to go about this?

Answer 1

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:

1. Use Gradient Descent by deriving the gradient of $ f \left( \cdot \right) $.
2. 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.
3. Solve this in Fourier Frequency Domain.

I previously solved similar problems using the Matrix Form (See my answer to Deconvolution of 1D Signals Blurred by Gaussian Kernel) 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:

1. The SNR is not good enough to solve the inverse problem as is. It requires additional regularization like Wiener / Tikhonov (They are the same). I implemented it in the code so you have Wiener Filter to play with and understand its working.
2. 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 (Look at the SignalProcessing\Q51460 folder).