3
$\begingroup$

My problem is the following, I have 3 curves/signals (1D) , the measure, the signal and the resolution of my detector: $\mathcal{M},\mathcal{S},\mathcal{R}$, knowing that : \begin{equation}\label{eq:four1} \mathcal{M} = \mathcal{S} * \mathcal{R} \end{equation} with $*$ being the convolution product, and my goal is to find out what $\mathcal{R}$ is, I could 'deconvoluate' but don't know how to so I though about the properties of the Fourier transform ($\mathfrak{F}$) in the frequency domain: \begin{align} \mathfrak{F}\left(\mathcal{S} * \mathcal{R}\right) &= \mathfrak{F}\left(\mathcal{S}\right) \cdot \mathfrak{F}\left(\mathcal{R}\right) \\ &= \mathfrak{F}\left(\mathcal{M}\right) \end{align} And so we get our resolution (a gaussian): \begin{equation}\label{eq:resol_eq} \mathcal{R} = \mathfrak{F}^{-1} \left( \frac{\mathfrak{F}\left(\mathcal{M}\right)}{\mathfrak{F}\left(\mathcal{S}\right)}\right) \end{equation} ($\cdot$ is the 'classical product) Is this something correct to do? In my case I am working with 1 D curves, both $\mathbb{M}$ and $\mathbb{S}$ are sigmoids and to get their fft I use numpy (python library): np.fft.fft(sigmoid_curve) and then I use ifft to get the inverse Fourier transform and finally get my $\mathbb{R}$.

When I test it I get indeed a gaussian but with the wrong $\sigma$ (I only care about this parameter) and I get some sort of repetitive pattern as you can see in the Figure below (light blue curve). Maybe it comes from the fact that I use fft's and not a proper Fourier transform ? Thanks

Looking at the plot, my goal is to find the best gaussian (purple) so that the green curve (S) fits the red one (M). Plot python

Thanks to Gideon's answer I could come up with this:

#time scales for the sigmoids
x_s = np.arange(-100, 100, 0.01)
x_r = np.arange(-950, 105, 0.01)

# source and measure
s = 1/(1 + np.exp(-x_s)) + 1
m = 1/(1 + np.exp(-.5*x_r)) + 1.1 #is the convolution of m and s 

#direct solution as in the presented example
S = linalg.toeplitz(s, np.concatenate([[s[0]], np.zeros("r???".size-1)]))
assert np.allclose(m, np.matmul(S, "r ???"))
r_estimated = linalg.lstsq(S, m)[0]
assert np.allclose(r_estimated, "r???")

#fast and efficeint solution (Levinson-Durbin recursion)
r_estimated_levin= linalg.solve_toeplitz((s, np.zeros(s.size)), m)[:"r???".size]

since I do not have access to r (I am actually looking for it) how could I do this ? like so ?

S = linalg.toeplitz(np.concatenate([[s[0]], np.zeros(len(m) - 1)]), s)
$\endgroup$
11
  • 1
    $\begingroup$ One issue that immediately pops out is in this expression: \begin{equation}\label{eq:resol_eq} \mathbb{R} = \mathfrak{F}^{-1} \left( \frac{\mathfrak{F}\left(\mathbb{M}\right)}{\mathfrak{F}\left(\mathbb{S}\right)}\right) \end{equation} If $\mathfrak{F}\left(\mathbb{S}\right)$ takes small values in certain ranges then you should try and exclude those ranges from your division. In essence, dividing by a very small number ($\lim_{x\rightarrow 0} x$) may lead to extremely large values or undefined results. $\endgroup$ Jul 20, 2023 at 17:16
  • $\begingroup$ Oh yes I forgot about that what a shame great idea thank you, apart from that the 'formula' is correct right ? $\endgroup$
    – Michael
    Jul 20, 2023 at 17:21
  • $\begingroup$ Yes but in the line before that you have a typo I believe. Convolution in one domain is multiplication in the other domain. $\endgroup$ Jul 20, 2023 at 17:24
  • $\begingroup$ Indeed thanks again I just edited it $\endgroup$
    – Michael
    Jul 20, 2023 at 17:26
  • $\begingroup$ So I checked and now I never divide by 0 but still got the same figure, moreover the gaussian I am expecting has to be centered in the inflexion point of my sigmoid which is the case so I am now wondering why there are many gaussian with the wrong $\sigma$ $\endgroup$
    – Michael
    Jul 20, 2023 at 18:57

1 Answer 1

3
$\begingroup$

$$\mathbb{M} = \mathbb{S} * \mathbb{R}$$

might be explicitly written as

$$\mathbb{M}[n]=\sum_{m=0}^{N-1}\mathbb{S}[n-m]\mathbb{R}[m]$$

So, assuming $n>m$, you have $n$ equations with $m$ parameters to find. The solution might be found using the least squares method.

To solve this system, it is common to use the matrix form, with the Toeplitz matrix.

Edges could be handled with respect to the edge conditions or neglected.

edit

Following the request, I attached a code.

import numpy as np
from scipy import linalg, __version__
import matplotlib.pyplot as plt

#semitimescales
x = np.arange(-100, 100, 0.01)
x2 = np.arange(-15, 15, 0.01)

# source (I have added a bias to allow further inversion of the Toeplitz)
s = 1/(1 + np.exp(-x)) + 1
# semi-goussian
r = np.exp(-x2**2/30)
r /= r.sum()
#convolution
m = np.convolve(s, r)[:x.size]

#plotting
plt.figure()
plt.plot(x, s)
plt.plot(x2, r/r.max())
plt.plot(x, m)
plt.grid()
plt.legend(['S', 'Normalized R', 'M'])
plt.show()


#direct solution as in the presented example
S = linalg.toeplitz(s, np.concatenate([[s[0]], np.zeros(r.size-1)]))
assert np.allclose(m, np.matmul(S, r))
r_estimated = linalg.lstsq(S, m)[0]
assert np.allclose(r_estimated, r)

#fast and efficeint solution (Levinson-Durbin recursion)
r_estimated = linalg.solve_toeplitz((s, np.zeros(s.size)), m)[:r.size]
assert np.allclose(r_estimated, r)

print(np.__version__)
print(__version__)

enter image description here

$\endgroup$
8
  • $\begingroup$ Okay but then using the Toeplitz matrix I can perform a 'deconvolution' to obtain the $\mathbb{S}$ ? Could you give me some insights on how to use the least square methods to obtain my signal ? $\endgroup$
    – Michael
    Jul 25, 2023 at 16:46
  • $\begingroup$ @Michael, Have a look at dsp.stackexchange.com/questions/23350, dsp.stackexchange.com/questions/74419 and dsp.stackexchange.com/questions/55284. You have a code attached there. $\endgroup$
    – Royi
    Jul 25, 2023 at 19:03
  • $\begingroup$ Yes I've seen those posts but I don't know how to implement those matrices calculations for 1D signals like you did but in python, in my case i need to find my matrix (NxN where N is the length of the region I am interested in: the inflexion point of my sigmoid right ? )but then I how can I get the $\mathbb{R}$ matrix using the equation above ? $\endgroup$
    – Michael
    Jul 26, 2023 at 21:26
  • $\begingroup$ @Michael, I have attached a code. Please tell me if there is anything preventing you from accepting the answer... $\endgroup$ Jul 27, 2023 at 18:24
  • $\begingroup$ @GideonGenadiKogan I did not had enough credit to validate.. now I do , thank you very much it was very basic indeed using spicy functions... really impressive I tried with convolution processes but always had a problem with the boundaries and all but now it works, I compared r and r_estimated and they are indeed the same (with different amplitudes). $\endgroup$
    – Michael
    Jul 28, 2023 at 20:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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