2
$\begingroup$

Lets say I have 2 vectors (1D signals that are sigmoids): $s$ and $m$, both related through the relation: $m = s * r$, my goal here is to recover the vector $r$ (should look like a gaussian).

I tried to use Fourier transforms but it did not work very well (see here: https://dsp.stackexchange.com/questions/88801/convolution-and-fourier-transform-for-1d-signals), then tried to use deconvolve from scipy but once again I do not get the right gaussian.

Then someone told me on my last post to look into Toeplitz matrices so we have the following relationship (discrete convolution): $$ m[n] = \sum_{i=0}^{n-1} s[n-k] r[k]$$ we then have n equations with $k$ parameters $\rightarrow$ matrix problem to solve. You'll find below my code using this post: Relationship between discrete deconvolution and Toeplitz matrices. I have a problem with my python code as the output does not give a proper gaussian, and if I i fit the peak that looks like one I always get a 2/3 ratio of difference for the standard deviation (see the plots I give both $\sigma$), how could I explain this ?

import numpy as np
import scipy.signal
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def sigmoid(x, A, x0, k):
    return A / (1 + np.exp(k * (x - x0)))

def gaussian(x, mu, sigma):
    return np.exp(-((x - mu) ** 2) / (2 * sigma ** 2))

x = np.arange(0., 20.01, 0.01)

# Create a Gaussian signal
mu_gaussian = 10.0
sigma_gaussian = 1.0
y = gaussian(x, mu_gaussian, sigma_gaussian)

# Create a sigmoid-like filter
A_sigmoid = 1.0
x0_sigmoid = 10.0
k_sigmoid = -3.0
sigmo = sigmoid(x, A_sigmoid, x0_sigmoid, k_sigmoid)

sigmo2 = sigmoid(x, A_sigmoid, x0_sigmoid, -0.5)

plt.plot(sigmo, label = '1st sigmo')
plt.plot(sigmo2, label = 'convolved sigmo')
plt.legend()
plt.show()
# yc = scipy.signal.convolve(y, c, mode='full') / c.sum()
# ydc, remainder = scipy.signal.deconvolve(yc, c)
M = len(sigmo2)

# Create the first row of the Toeplitz matrix using the sigmoid values
first_row = np.concatenate((sigmo, np.zeros(M - 1)))
# Create the first column of the Toeplitz matrix with zeros and the first element of `sigmo`
first_col = np.concatenate(([sigmo[0]], np.zeros(M - 1)))

# Generate the Toeplitz matrix
A = scipy.linalg.toeplitz(first_col, first_row)

# # Calculate the product of transpose(A) and A
# product = np.dot(A.transpose(), A)
# # Calculate the inverse of the matrix product
# inverse_product = np.linalg.inv(product)
# # We get the moore penrose pseudo inverse
# Penrose = np.dot(inverse_product, A.transpose())

Penrose = np.linalg.pinv(A)
#Finally our retrieved gaussian
gauss_toep = np.dot(Penrose, sigmo_conv)

# Create new x values matching the length of gauss_toep
x_fit = np.linspace(0, 20, len(gauss_toep))

plt.plot(x_fit, gauss_toep)
plt.show()

popt_gauss, _ = curve_fit(gaussian, x_fit, gauss_toep, p0=[0.0001, 0.001, 15, 1])

print("$\sigma$ = ", popt_gauss[-1])

# Generate the fitted gaussian
gauss_fit = gaussian(x_fit, *popt_gauss)

plt.plot(x_fit, gauss_toep, label = r'Retrieved gaussian, $\sigma$ = {:.3f}'.format(sigma_gaussian))
plt.plot(x_fit, gauss_fit, label = r'Fitted retrieved, $\sigma$ = {:.3f}'.format(popt_gauss[-1]))
plt.legend()
plt.show()

Input and output sigmoids (the orange is the convolution of the blue and a gaussian)

In blue the output of the toeplitz and Least square method, orange the fit of the gaussian, we have a 2/3 difference between the original gaussian (used for the convolution)

$\endgroup$
5
  • $\begingroup$ Take a look at this $\endgroup$ Aug 4, 2023 at 12:46
  • $\begingroup$ I corrected thanks, now regarding the other post they are solving the linear equation without the toeplitz matrix but would it give the same result ? Cause the code I gave depends also on the window size we choose for the convolution , for a different window we would get different results for $r$ $\endgroup$
    – Michael
    Aug 4, 2023 at 12:51
  • $\begingroup$ Why don't you steal the MathJax code from the other post and express your question mathematically, rather than express it in Python? $\endgroup$ Aug 4, 2023 at 13:04
  • 2
    $\begingroup$ Does this answer your question? Compensating Loudspeaker frequency response in an audio signal This specifically shows and further derives how deconvolution is done using the Toeplitz matrix resulting in a least squares solution (along with an example). $\endgroup$ Aug 5, 2023 at 5:05
  • $\begingroup$ Well it surely helps for the process but I still find some problems when plotting my gaussian (I've edited the post to show my results) $\endgroup$
    – Michael
    Aug 7, 2023 at 16:32

1 Answer 1

1
$\begingroup$

You model is given by:

$$ \boldsymbol{m} = \boldsymbol{s} \ast \boldsymbol{r} $$

Where $ \ast $ is the convolution operator.
Now, since it is a linear operator, the above can be written in a matrix form:

$$ \boldsymbol{m} = \boldsymbol{S} \boldsymbol{r} $$

Where the matrix $ \boldsymbol{S} \in \mathbb{R}^{m \times n} $ is the model matrix.
Now, we need to take into account 2 properties of the matrix:

  1. Tall / Fat Matrix
    What are the relation between $ m $ and $ n $.
    We prefer $ m \gg n $. In other cases we'll must use regularization (Though it is recommended in any case if we have prior knowledge).

  2. Convolution Model
    Discrete convolutions have different way to handle boundaries and what samples are extracted from the result.
    The common forms of the convolution are: full, same and valid. Those are basically the different approaches to chose the output samples from the full convolution.
    Boundary are usually one of: constant, replicate and cyclic.
    Usually, when we apply explicit boundary conditions, we use the valid type of convolution on the padded signal.
    When no explicit choice of boundary conditions is made, we usually use the full convolution. You may look at Solving linear equation of discrete convolution kernels using black box model for the convolution, The Different Solutions for Filter Coefficients Estimation for Periodic Convolution and Full Convolution and How to Solve a Composition of Convolutions from Regularized Least Squares Model in Frequency Domain.

When we have the choice from above, it is easy to generate the matrix model and then basically solve a linear system.

The matrix, in some of the cases above (Like full or periodic) the model matrix is a Toeplitz matrix.
Being Toeplitz doesn't improve the results, it is only improvement to the calculations. Both applying the matrix or solving the system.
So if you have numerical issues with the explicit matrix, solving it as a Toeplitz system won't make a difference.

Solving specifically for a Gaussian Like kernel is done in 1D Deconvolution with Gaussian Kernel (MATLAB). Pay attention to the condition number which is a limiting factor of getting the solution.

$\endgroup$
5
  • $\begingroup$ Thank you for your help, I used your comments and links but I still face some weird problems when plotting my gaussian, do you have an idea of where the problem can come from looking at my code or plots maybe ? $\endgroup$
    – Michael
    Aug 7, 2023 at 16:31
  • $\begingroup$ @Michael, Did you check the condition number? $\endgroup$
    – Royi
    Aug 15, 2023 at 8:15
  • $\begingroup$ I tried to look at your post (1D Matlab) but, as shown in the python code, I did not checked it as I do not know how to calculate it and find the optimal one ? $\endgroup$
    – Michael
    Aug 15, 2023 at 17:34
  • $\begingroup$ @Michael, There is no optimal one. Yet if it is high, it means you won't be able to recover the data accurately. $\endgroup$
    – Royi
    Aug 15, 2023 at 18:46
  • $\begingroup$ I calculated it with the formula: $\sqrt{\mathbf{S}^T \cdot \mathbf{S}}$ where $\mathbf{S}$ but lets say I get 40 is there a way to tell that this value is acceptable and thus my problem is not ill posed ? $\endgroup$
    – Michael
    Aug 15, 2023 at 20:03

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