# Deblurring 1D data using direct inverse filtering

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] 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)/F(h) 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.

• Will you accept MATLAB code? Could you post the data (As csv file) of both the temperature and the kernel?
– Royi
Jul 17, 2021 at 12:54
• @Royi added data.csv and the impuslse response Help in matlab is fine too thanks Jul 17, 2021 at 13:01
• Dividing by the transfer function is a bad idea since you will get "divide by zero" or "divide by too small a number" problems. $h$ is a lowpass filter and hence it's not invertible. You can't recover the original signal. Jul 17, 2021 at 19:05
• @MitulAgrawal, Are you sure the data is correct? Do you have the code that generated it?
– Royi
Jul 18, 2021 at 9:45
• @Royi yes the data is correct (x[n] : true values and y[n] : recieved values), Unfortunately I don't have the code that generated it. Thanks for the Answer :) Jul 18, 2021 at 11:23

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).

• Nice fan club, started mine Jul 20, 2021 at 7:16
• @OverLordGoldDragon, What do you mean?
– Royi
Jul 20, 2021 at 10:59
• This - neat, really, a win-win. But it's also why I avoid +1 your posts. If you're unlucky, might chime in. Jul 20, 2021 at 11:30
• I am not sure I get it. Why do the codes make you avoid +1? If I see your answer and it is valuable (Like all posts of your I saw so far) I up vote. What do you mean by Win Win?
– Royi
Jul 20, 2021 at 17:51
• I am not sure what you mean with unfair manner? You think people coming from GitHub to vote? I think your voting calculations mean the others, by your perception, are doing right by being very hard on the +1. Be generous and others will be with you.
– Royi
Jul 21, 2021 at 14:41