# How to deconvolve matrix when a model of noise exists?

I have a matrix of N rows of time-series data. There is a specific noise contaminating measurement of the data that I have some information about.

The noise in the data can be modeled as a poisson distribution that blurs signal from a given column in the matrix to adjacent columns. For example, if the original data should be a single peak surrounded by no signal:

0    0    0    1    0    0    0


The measured signal distributed slightly asymmetrically resulting in something like this:

0.001    0.005    0.1    0.5    0.2    0.001    0


If I have a good model of how the noise is distributing the data between the columns, how can I use this information to deconvolve the matrix into an approximation of the original signal?

This isn't (as far as I am aware) a standard approach to deconvolution, but it seems to address your problem:

• Assume that your original signal $s$ is distorted by a circulant matrix $\mathbf{C}$ made up of the (shifted) components of your measured signal $m$: $$m = \mathbf{C} s$$
• Then just find the generalized inverse of $\mathbf{C}$, denoted $\mathbf{C}^\dagger$ and form your estimate $\hat{s}$ as: $$\hat{s} = \mathbf{C}^\dagger m$$

The image below shows:

• $\color{black}{\tt black}$ : the original signal in your question.
• $\color{red}{\tt red}$ : the distorted signal in your question.
• $\color{green}{\tt green}$ : 100 realizations of inverting a noisy version of your distorted signal.

As you can see, the green estimates bounce around the original noiseless signal reasonably well. R Code Below

# 31682

# http://stackoverflow.com/a/15796694/12570
circ<-function(x) {
n<-length(x)
matrix(x[matrix(1:n,n+1,n+1,byrow=T)[c(1,n:2),1:n]],n,n)
}

original <- c(0,0,0,1,0,0,0)
distorted <- c(0.001,0.005,0.1,0.5,0.2,0.001,0)

inverse_matrix <- ginv( circ(distorted[c(4,3,2,1,7,6,5)]) )

Nruns <- 100

output <- inverse_matrix %*% distorted

plot(original, type="l",
ylim=c(min(c(original,distorted,noisy_distorted, output)),max(c(original,distorted,noisy_distorted, output))),
xlim=c(1,length(output)), lwd=5)
lines(noisy_distorted, col="red", lwd=5)

for (runNo in 1:Nruns)
{
noisy_distorted <- distorted + rnorm(length(distorted), 0, 0.01)
output <- inverse_matrix %*% noisy_distorted
lines(output, col="green")
}

lines(original, lwd=5)

title('Original, distorted, and estimated originals (100 realizations) ')

• Thanks, I will play with this and get back to you. I appreciate your help and time. – Justin G Jun 22 '16 at 21:37
• @JustinGardin let me know if there are other issues. – Peter K. Jun 22 '16 at 22:10
• I don't understand the reasoning behind the line: ginv( circ(distorted[c(4,3,2,1,7,6,5)]) ) I understand that you are permuting the distorted array, turning this array into a circulant matrix, then finding the inverse of that matrix. Why do you have the arrangement in this order 4,3,2,1,7,6,5? Does it start with 4 because 4 is the location of the original signal? Sorry, I don't have great mathematics skills and the wikipedia page only makes superficial sense to me. – Justin G Jun 23 '16 at 10:12
• Yes, that's all: the first row of the matrix has to be rotated and reversed to get the right output for your input. The matrix multiplication is just convolution. – Peter K. Jun 23 '16 at 10:22
• Lets imagine I don't know precisely which point is the original signal, or that there are multiple signals in one timeseries vector. I tried multiple test distorted vectors, and this seems to work for one signal. Is this a property of this solution? For example distorted <- c(0.001,0.005,0.1,0.5,0.2,0.001,0,0.001,0.005,0.1,0.5) inverse_matrix <- ginv( circ(distorted[c(4,3,2,1,12,11,10,9,8,7,6,5)]) ) output <- inverse_matrix %*% distorted only returns one peak. Moreover, if I change the values to lower the intensity of the signal, the identified peak still returns a 1. – Justin G Jun 23 '16 at 10:53