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)
Improve this question
$\endgroup$
10
  • 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$ Commented 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$ Commented 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$ Commented Jul 20, 2023 at 17:24
  • $\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$ Commented Jul 20, 2023 at 18:57
  • $\begingroup$ I am not sure but my best guess is that it could be a result of the periodic assumption inherent in the FFT. If your signal isn't periodic, the FFT implicitly "repeats" your signal to make it so, which can result in strange effects in the frequency domain. To mitigate these issues, you might want to carefully preprocess your signals (e.g., apply a window function, zero-pad to avoid circular convolution effects, etc.), and consider using a regularized deconvolution method to handle the division in the frequency domain. $\endgroup$ Commented Jul 20, 2023 at 20:14

1 Answer 1

Reset to default
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

Improve this answer
$\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$ Commented 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$ Commented 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$ Commented 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$ Commented 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$ Commented Jul 28, 2023 at 20:25

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.