Deblurring 1D data using direct inverse filtering

Question

In my assignment I have been given recorded temperature over a period of time (193 values) and the impulse response (5 values with n=0 corresponding to the middle value.)

Data: data.csv

h = [1/16 4/16 6/16 4/16 1/16]

given question:

This is what I need to do for deblurring -> take Fourier transformation of x and h > F(x)/F(h) > inverse Fourier transformation

This is my code for the Fourier transformation :

def dtft (x) : 
    N = len(x)
    m = [[np.exp(-2j * np.pi * i * j / N) for i in range(N)] for j in range(N)]  
    return [sum(m[i][j]*x[j] for j in range(N)) for i in range(N)]

I am facing difficulty in the F(x)/F(h) part.

Firstly do I add zeros in h to make its length equal to x? How do I make use of the statement - "n=0 corresponds to the middle value in h?"

Secondly for F(x)/F(h) should I divide the corresponding frequencies of F(x) by F(h)? For example: F(x)[0]/F(h)[0] and so on. (I tried doing this but got weird results, maybe due to improper padding of h.)

Is there any other method than F(x)/F(h) like finding g where I can directly get the result by convolving x and g?

Any help would be appreciated as I am new to the topic.

Answer 1

Doing direct division in Frequency Domain means you are assuming cyclic / periodic boundary conditions for your data which I don't think is the proper assumption for your data.

The model I'd pursue would be as I described in my answer to Deconvolution of Synthetic 1D Signals - How To?
As one can see in that question the marked answer by @Hilmar is wrong. Another good reference is Deconvolution of a 1D Time Domain Wave Signal Convolved with Series of Rect Signals.

The Model

$$ f \left( x \right) = \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2} $$

Where $ h $ is the filter signal, $ x $ is the signal we're after and $ y $ is the measured signal.

The optimization problem is given by:

$$ \arg \min_{x} f \left( x \right) = \arg \min_{x} \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2} $$

Assuming White Noise finding the signal $ x $ which minimizes the function above will yield the optimal solution.

In order to solve it we'll convert the problem into a matrix form:

$$ \arg \min_{x} f \left( x \right) = \arg \min_{x} \frac{1}{2} \left\| H x - y \right\|_{2}^{2} $$

We can build the matrix $ H $ using the function I built: CreateConvMtx1D(). I used convolution shape of same as it seems to be implied by the question.

Sensitivity

Since the solution is given by $ {\left( {H}^{T} H \right)}^{-1} {H}^{T} y $ the condition number of $ {H}^{T} H $ is practically the sensitivity of the model.

The Solution

The Condition Number of $ H $ is relatively high (~4.6209e+07). Hence we have 2 options:

1. Understand the problem is highly sensitive, use Float64 for calculations and hope the noise level is low. Otherwise the solution will have high noise since the condition number implies about the noise simplification (Also about how close the filter is to zero in frequency domain).
2. Use regularization as a combined step of denoising and deconvolution (Denoising and Sharpening in question lingo).

For the regularized case (Equivalent of Wiener Deconvolution, See How Is the Formula for the Wiener Deconvolution Derived?) the optimization problem becomes:

$$ \arg \min_{x} f \left( x \right) = \arg \min_{x} \frac{1}{2} {\left\| H x - y \right\|}_{2}^{2} + \frac{ {\lambda}^{2} }{2} {\left\| x \right\|}_{2}^{2} $$

The solution is given by $ {\left( {H}^{T} H + {\lambda}^{2} I \right)}^{-1} {H}^{T} y $.

The Question

All above is correct yet each case requires its own treatment.

Looking at the signals:

One can see the filter didn't do much (Even look at the sharp rise). Yet the noise level is very high (The STD is ~2).

Looking at at the DFT of the filter above:

The above figure (DC Component was removed) only confirms, the LPF has almost no effect in the frequencies the data really exist. So deconvolution will only "fight" noise which will yield very bad results (As it will be amplified).

For this kind of problem the best thing (Which is the essence of question) to do is reduce the noise. You may do that with a filter. Yet a sharp filter for this problem will require many coefficients and won't be proper for this limited number of samples.

This is a noise governed problem and without a good prior for the data, deconvolution won't work.

The code is available at my StackExchange Signal Processing Q76333 GitHub Repository (Look at the SignalProcessing\Q76333 folder).