1
$\begingroup$

I am an experimental research engineer working on a Laser experiment. The details do not matter but due to practical limitations my data is influenced by the following transformation during data acquisition, and I was wondering if there is a mathematical method or an algorithm to "undo" its effect on the captured signal, that is to get the inverse. enter image description here

I want to get the y in terms of s. Thanks in advance

$\endgroup$
2
  • $\begingroup$ Are the $y_i$ equispaced samples, i. e. $y_i=y(iT)$? $\endgroup$ – Deve Dec 23 '13 at 9:06
  • $\begingroup$ Yes it is equispaced $\endgroup$ – Aladdin Dec 23 '13 at 18:42
1
$\begingroup$

Looks like a convolution. S is the convolution of y with another signal of w+1 consecutive ones. So, S is the result of a moving sum filter. But the moving sum filter completely removes certain frequencies in the signal. For example, if you use w=3 and have the following signal y:

 y: ... 0  1  0 -1  0  1  0 -1  0  1  0  -1 ...

Then S[i] for any i (the moving sum over 4 neighboring samples) will always be zero, no matter how you shift or scale y. And that makes the problem of recovering y from s something between impossible and hard.

In case your data is finite, it is possible to write your equations as a matrix-vector product and solve it with the usual algorithms. You can ask tools like Matlab or GNU Octave about how noisy the inverse would be. This is something you should check for since the problem is ill-conditioned. Your equations, assuming y and s are finite can be written like this:

$$ A \cdot y = s $$

Here, A would contain mostly zeros and w+1 consecutive ones per row, at least for first rows. At the end, you have less ones because y is supposed to be finite. Then, asking for the condition number is a matter of typing

cond(A)

in Matlab or Octave. Example for 6 samples and w=2

octave:1> A = toeplitz([1 0 0 0 0 0],[1 1 1 0 0 0])
A = 

   1   1   1   0   0   0
   0   1   1   1   0   0
   0   0   1   1   1   0
   0   0   0   1   1   1
   0   0   0   0   1   1
   0   0   0   0   0   1

octave:2> cond(A)
ans = 7.9053

Let's actually test it:

octave:3> y = [1 9 3 2 4 3]';
octave:4> s = A*y
s =

   13
   14
    9
    9
    7
    3

And now the inversion:

octave:5> y2 = A\s
y2 =

   1
   9
   3
   2
   4
   3

This inversion is actually easy to implement because the matrix exhibits a triangular structure. You can recover y sample by sample starting at the end and moving to the front.

This example was free of noise. Keep in mind that the condition number tells you how the relative error in s will be amplified (worst case) if you recover y from s. So, for example, if s contains a relative error of 3%, the recovered y may have a relative error as large as 3% * 7.9053 = 23.716%. Unfortunately, the condition number increases with the number of samples in this case. So, if you have to recover not six but 1000 samples and/or your w parameter is large, you're out of luck because the inversion will amplify the noise in s a lot, too.

$\endgroup$
4
  • $\begingroup$ Tried and the output was noise with peaks ! What does y2 = s'\A; mean? I tried it and it gave me something that vaguely looked like the correct answer (or may be not) $\endgroup$ – Aladdin Dec 23 '13 at 20:04
  • $\begingroup$ I wrote A\s, not s\A. You can think of a/b as a*inv(b). The blackslash is just a division in which the left hand side is inverted. So, A\s is something like inv(A)*s but it's actually preferable to inv(A)*s due to better numerical accuracy. $\endgroup$ – sellibitze Dec 24 '13 at 10:28
  • $\begingroup$ I know but A\s did not work for me s\A did :D I dont know why Also : Why w+1, and not just w ? Thank you again for your time and concern $\endgroup$ – Aladdin Dec 24 '13 at 20:27
  • $\begingroup$ @Aladdin: A\s should work if s is a column vector. As for w+1: A loop from i to i+w contains w+1 iterations. $\endgroup$ – sellibitze Dec 24 '13 at 20:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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