# Deconvolving a 1d Signal Using a Lookup Table of Kernels

assuming I measure a signal that has different PSFs per position in time.

for example:

t = linspace(0,20);

% "ground truth" signal will be like:

x = @(t0) exp(-(t-t0).^2/0.1)  ;

% some made up impulse response (psf) that depends on t0 will be:

h = @(t0) diff(diff(  exp(-(fftshift(t)).^2./(0.1*t0) )));

% the convovled signal:
y = @(t0) conv(x(t0),h(t0) ,'same');

% now if I have signal from two positions, I get:
plot(t,y(3)+y(15))


Note that the two peaks now are distorted differently as function of their position.

What methods can I use here, given that I have such a PSF lookup table, such as that h = @(t0)... above, to deconvolve my 1D signal, even though it will behave differently in different positions as seen in the plot? just a standard deconvolution wont work here.

EDIT : trying to clarify the question further. I am looking after a way to "deconvolve" the signal that is distorted by such position dependent PSF. So instead of these two features I will be able to trace back the original signal (this case just two peaks). Using standard de-convolution schemes will not work well because they assume an effective single PSF, and here we have a "family" of PSFs. Is there a way to solve it? I was hoping, for example, that extending the dimension of the PSF will allow to accommodate such effect, or maybe using other tools to "train" a system to understand it.

EDIT 2: Here's a file that show an example of x - the ground truth signal, y - the signal convoluted by the positions dependent psfs (or kernels) and, psfs - an array of kernels per position.

• Could you please make the question clearer. I am not sure what are you after. – Royi Jul 20 '20 at 7:01
• I just added an edit, I hope this clarify further my question. – bla Jul 20 '20 at 7:11
• I still don't understand it. Could you explain the math? At the end, if the output is given by a convolution what basically are you after, the parameter t0? – Royi Jul 21 '20 at 9:05
• Instead of having one PSF in the convolution of the entire signal y(t)=conv(x(t),psf) we have the set of PSFs (one for each location in t). So, you can think of this as breaking the convolution over the entire signal grid t to many contributions, for example, each k-th element in the signal is dictated by a different PSF : y(t_k ) = conv (y( t_k) , PSF_k ) , etc. I am after a deconvolution \ deblur of the signal, to remove the effect of the different PSFs. In the example above, such de-blur will result in just two peaks, of the same widths and amplitude in the positions 3 and 15. – bla Jul 21 '20 at 21:34
• I answered many Deconvolution / Deblurring questions here and I still don't get your question. If your kernel is time variant the operation isn't convolution hence the problem isn't deconvolution. In images usually what we have is the assumption the kernel is spatially variant yet it changes on areas much bigger than its support. Then basically there is an ensemble of sub problems each is deconvolution. If in your case each sample of time has different kernel than this is a completely different operation from deconvolution. – Royi Jul 22 '20 at 6:26

The way I understand the problem is each sample of the output is a linear combination of the samples of the input.

Hence it is modeled by:

$$\boldsymbol{y} = H \boldsymbol{x}$$

Where the $$i$$ -th row of $$H$$ is basically the instantaneous kernel of the $$i$$ -th sample of $$\boldsymbol{y}$$.

The problem above is highly ill poised.

In the classic convolution case we know the operator matrix, $$H$$, has a special form (Excluding the borders) - Circulant Matrix. With some other assumptions (Priors) one could solve this ill poised problem to some degree.

Even in the case of Spatially Variant Kernels in Image Processing, usually, some form is assumed (Usually being block circulant matrix, and the number of samples of each kernel is larger than the number of samples in the support of the kernel).

Unless you add some assumptions and knowledge into your model the solution will be Garbage In & Garbage Out:

numInputSamples = 12;
numOutputSamples = 10;

mH = rand(numOutputSamples, numInputSamples);
mH = mH ./ sum(mH, 2); %<! Assuming LPF with no DC change

vX = randn(numInputSamples, 1);

vY = mH * vX;

mHEst = vY / vX;


See the above code. You will always have a perfect solution yet it will have nothing to do with mH.

Now, if I get it right, you say I don't know $$H$$ perfectly, but what I have is a pre defined options.

So let's say we have a matrix $$P \in \mathbb{R}^{k \times n}$$ which in each row has a pre defined combination:

$$H = R P$$

Where $$R$$ is basically a row selector matrix, namely it has a single element with value $$1$$ in each row and the rest is zero.

Something like:

mP = [1, 2, 3; 4, 5, 6; 7, 8, 9];
mH = [1, 2, 3; 7, 8, 9; 7, 8, 9; 4, 5, 6; 4, 5, 6];

% mH = mR * mP;

mR = mH / mP;


So our model is:

\begin{aligned} \arg \min_{R, \boldsymbol{x}} \quad & \frac{1}{2} {\left\| R P \boldsymbol{x} - \boldsymbol{y} \right\|}_{2}^{2} \\ \text{subject to} \quad & R \boldsymbol{1} = \boldsymbol{1} \\ & {R}_{i, j} \geq 0 \quad \forall i, j \\ \end{aligned}

It is still exceptionally hard (Non convex) problem but with some more knowledge it can be solved by utilizing alternating methods where we break the optimization problem as:

1. Set $$\hat{\boldsymbol{x}}^{0}$$.
2. Solve $$\hat{R}^{k + 1} = \arg \min_{R} \frac{1}{2} {\left\| R P \hat{\boldsymbol{x}}^{k} - \boldsymbol{y} \right\|}_{2}^{2}$$ subject to $$R \boldsymbol{1} = \boldsymbol{1}, \; {R}_{i, j} \geq 0 \; \forall i, j$$.
3. Solve $$\hat{\boldsymbol{x}}^{k + 1} = \arg \min_{\boldsymbol{x}} \frac{1}{2} {\left\| \hat{R}^{k + 1} P x - \boldsymbol{y} \right\|}_{2}^{2}$$.
4. Check for convergence, if not go to (2).

Now each sub problem is convex and easy to solve.

Yet I still recommend you add better assumptions / priors.

Such as the minimal number of contiguous samples that must have the same PSF (Similar to 2D in images where we say each smooth area is smoothed by a single PSF).

Remark

We didn't employ the fact each element in $$R$$ is either 0 or 1 as the straight forward use of that will create a Non Convex sub problem.
In case the number of PSF's is small we can use MIP solvers. But the model above assumed each row is a PSF so for large number of samples even in case we have small number of PSF'w the matrix is actually built by shifting each PSF as well. So we'll have large number in any case.

Another trick might be something like Solving Unconstrained 0-1 Polynomial Programs Through Quadratic Convex Reformulation.
Yet the simplest method would be "projecting" $$R$$ into the space (Which is not convex, hence projection is not well defined). One method might be setting the largest value to 1 and zero the rest.

## Update

In comments you made it clear you know the kernel per output sample.
Hence the model is simpler:

$$\boldsymbol{y} = A \boldsymbol{x} + \boldsymbol{n}$$

The least squares solution is simply $$\boldsymbol{x} = {H}^{-1} \boldsymbol{y}$$.
For better conditioning and noise regularization (Actually prior about the data, but that's for another day) you can solve:

$$\hat{\boldsymbol{x}} = {\left( {A}^{T} A + \lambda I \right)}^{-1} {A}^{T} \boldsymbol{y}$$

This is a MATLAB code for proof of concept:

load('psfs.mat');

mA = psfs;
vY = y;
vX = x;

vParamLambda    = [1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 1];
numParams       = length(vParamLambda);

vValMse = zeros(numParams, 1);

mAA = mA.' * mA;
vAy = mA.' * vY;
mI  = eye(size(mA));

for ii = 1:numParams
paramLambda = vParamLambda(ii);
vEstX = (mAA + paramLambda * mI) \ vAy;

vValMse(ii) = mean((vEstX(:) - vX(:)) .^ 2);
end

figure();
hL = plot(vParamLambda, 10 * log10(vValMse));
xlabel('Value of \lambda');
ylabel('MSE [dB]');


This is the result:

• Thank you for the answer . I actually do know the instantaneous kernel (see the line "given that I have such a PSF lookup table" in the question). How will that change your answer? Not only I know the instantaneous kernel, it only changes at some local part of the signal. – bla Jul 22 '20 at 16:56
• Have you looked on the 2nd part of the answer? Do you mean you know for each samples of the output the PSF used? – Royi Jul 22 '20 at 17:02
• yes, you are addressing my point in the 2nd part. I am trying to understand how to translate the approach you write into code. some notations I'm not sure I understand, what x hat 0? the first element of x/norm(x)? what algorithms would you use to solve step 2 and 3? will adding the assumption that only between positions 0 and 8 I have a smoothly changing psf, and afterwards it can be approximated as a single psf help in that regard? is that considered a small # for the MIP solver you mentioned? – bla Jul 22 '20 at 18:12
• @bla, First, could you answer if you know the exact PSF per sample of the output ($\boldsymbol{y}$) or just it is one of the list you have? The $k$ above means iteration. We iterate doing step 2 and 3. I will point that. – Royi Jul 23 '20 at 4:20
• first I want to thank you again for the patience and effort. I deeply appreciate it. second, yes, if I understand you correctly, I have a mapping of what PSF would work at which position (or array element) of y. So, if y is a vector composed of delta functions at different position, I'll get what I showed in the question. Naturally, if these deltas are closer than the envelope width of the PSFs they generate, it will hard to point what were their original positions. This is what I'm trying to solve. – bla Jul 24 '20 at 2:28
1. If the signal is oversampled and the PSF variation corresponds (approximately) to a smooth local compression/expansion, perhaps you can resample y so as to make the PSF approximately LTI, then apply conventional methods (somewhat akin to homomorphic processing)

2. If the input signal is convolved with a small discrete set of PSFs, perhaps you can devonvolve the entire signal with all of them, then pick the output that corresponds to that region?

3. As a MATLAB guy I found this snippet interesting: http://eeweb.poly.edu/iselesni/lecture_notes/least_squares/LeastSquares_SPdemos/deconvolution/html/deconv_demo.html perhaps you can get by with something ala (depending on the numerical properties of your convolution matrix and your complexity requirements):

x = randn(3,1);

h = randn(3,3);

y = h*x;

x_hat = h \ (y+eps);

• tried #2 with limited success, and as for #1, unfortunately, the change is not as smooth but I actually can well characterize it by measuring it independently, cant this help? can you demonstrate how homomorphic processing will work here? – bla Jun 24 '20 at 6:53
• I am not sure that homomorphic is the right term here. That usually refers to making a nonlinear system solvable by linear means by applying a nonlinear pre/post-processor. What I suggested in 1 would be more like applying space/time-variant processing so as to make a s/t-variant system appear to be LTI. en.wikipedia.org/wiki/Homomorphic_filtering – Knut Inge Jun 24 '20 at 7:06
• Ivan Selesnick is an exellent DSP guy, before being a Matlab guy – Laurent Duval Jul 22 '20 at 12:20

As an answer probably require to have more details on the look-up table (smoothed and regularity of the kernels), here is a couple of recent papers, including a review:

• thank you for the relevant and insightful papers. some are under a pay wall unfortunately. Are the 2D methods presented there compatible with m 1D signal question as far as you can tell? – bla Jul 22 '20 at 6:37
• I mentioned the DOI links for stable referencing. Most have a preprint somewhere. 2D methods can transfer to 1D, it is more about kernel shape than dimensions. In 1D, the references I know of are more in geophysics, with papers more difficult to get. In both meanings – Laurent Duval Jul 22 '20 at 12:18