1

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?

Justin G
  • 31
  • 1
  • Have you tried the deconv function: http://www.mathworks.com/help/matlab/ref/deconv.html ? – Amal Jun 22 '16 at 15:46

1 Answers1

1

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.

enter image description here

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) ')
Peter K.
  • 25,714
  • 9
  • 46
  • 91
  • 1
    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
  • To get it to work for multiple signals, you'll have to do something like apply the above algorithm along the time series. Not knowing the precise location might be more problematic. Can you update your question with more details as to how the extra distorted signals enter and/or how long your time series is? That might allow me to modify the algorithm to work better for your more general case. – Peter K. Jun 23 '16 at 14:33