5
$\begingroup$

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:

enter image description here

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 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 0's 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)? eg. F(x)[0]/F(h)[0] .... and so on (tried doing this but got weird results, maybe due to improper padding of h)

Also is there any other method than F(x)/F(h) like finding g where i can directly get the result by convulating x and g?

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

$\endgroup$
5
  • $\begingroup$ Will you accept MATLAB code? Could you post the data (As csv file) of both the temperature and the kernel? $\endgroup$
    – Royi
    Jul 17 at 12:54
  • $\begingroup$ @Royi added data.csv and the impuslse response Help in matlab is fine too thanks $\endgroup$ Jul 17 at 13:01
  • $\begingroup$ 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. $\endgroup$
    – Hilmar
    Jul 17 at 19:05
  • $\begingroup$ @MitulAgrawal, Are you sure the data is correct? Do you have the code that generated it? $\endgroup$
    – Royi
    Jul 18 at 9:45
  • $\begingroup$ @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 :) $\endgroup$ Jul 18 at 11:23
4
$\begingroup$

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:

enter image description here

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:

enter image description here

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.

$\endgroup$
9
  • $\begingroup$ Nice fan club, started mine $\endgroup$ Jul 20 at 7:16
  • $\begingroup$ @OverLordGoldDragon, What do you mean? $\endgroup$
    – Royi
    Jul 20 at 10:59
  • $\begingroup$ This - neat, really, a win-win. But it's also why I avoid +1 your posts. If you're unlucky, might chime in. $\endgroup$ Jul 20 at 11:30
  • $\begingroup$ 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? $\endgroup$
    – Royi
    Jul 20 at 17:51
  • $\begingroup$ I'm saying you've secured a baseline vote with Watchers. In an ideal world, I'd +1 if I simply find an answer legitimate, regardless of its existing votes - but in this network I can shit bricks and still get nill. Worse, it's rife with absurd vote-to-effort ratio posts. This makes me withhold my vote if anyone already has 2, unless it's exceptional. Unlike with others, I do find your posts to be high quality and fair effort - but put simply, you're already 'compensated' for it, and in a somewhat unfair manner (though I don't mind). $\endgroup$ Jul 21 at 12:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.