deconvolution in time?

Hello!
I have a output signal y(k) and a plant P(z). The signal y(k) contains 
some noise and I know the PSD of the noise. But lets assume that for now 
that the noise is qhite and the system looks like this:

Y(z)=U(z)P(z)+E(z)

u(k)                      y(k)
---------->[P(z)]-->[+]--->
                      |
                e(k)  |
               -------|

I want to restore the signal u(k) from y(k). I can do this off-line in 
freq domain, no problem, but I need this on-line. I think I need 
something related to Wiener filters.

The deconvolution P^(-1)(z) may be unstable, which is the biggest 
problem for me. I am not worried that it may be no casual.

I read that there are Kalman filter based approaches to this, but I 
can't find some practical explanation I can understand in reasonable 
amount of time.

Can somebody give me some reading directions - books, papers. Better 
papers, because I can download most of them in PDF.

Thanks a lot!


--JS

Dear Julian,

some years before I thought about a similar problem and found a
solution:

The convolution can be considered as a multiplication of a vector y(k)
and a triangular (Toeplitz) matrix that is built up from the plant
P(k). For a better understanding I try to develop this in formulae.

y(k) = x(k) * P(k)          // I use * as symbol for convolution

y(k) = SUM(x(k) P(k-n))     // This is the classic discrete time
convolution formula: n is the 'running' variable, the sum is built up
for all n from 0 to k-1.

the formula above is the said vector/matrix multiplication: y(k) =
x(k)P(k) - You did correct by using a capital P for the plant to symbol
a matrix :-)

y(0)      |P(0)   0    0  ...      0|      x(0)
y(1)      |P(1)  P(0)  0 ..        0|      x(1)
y(2)      |P(2)  P(1)  P(0)..      0|      x(2)
.      =  |P(3)  P(2)  P(1) P(0).. 0|  x     .
.         |P(4)     ...            0|        .
.         | .........              0|        .
y(k-1)    |P(k-1) P(k-2) ...    P(0)|      x(k-1)

The Matrix (hope it got not scrambled too much ;-)) is very
straigtforward: the diagonal is occuied with all P(0)es.

Now you want to extract the vector x(k). The first 'idea' is to
'devide' the vector y(k) by the matrix P(k). O.K. this is done by
matrix inversion of P(k):

x(k) = y(k) P(k)^-1

The matrix inversion is done by multiplying the inverse determinant
with the matrix itself. Before you scratch your head rethinking this
sentence here's the formula:

P(k)^-1 = 1/det(P(k))  P(k)

To have a matrix invertable there have to be some conditions to be
fulfilled: The determinant of the matrix P(k) must not be zero
(otherwise the inverse matrix would have all infinities :-)). How to
determine the determinant? Before you read the entire math library of
your university: All upward and downward diagonal products that have to
be summed up and substracted in the crazy manner you have learned at
High School will result to ZERO since you have the upper right triangle
of your Toeplitz Matrix occupied with zeroes! The only exception is the
diagonal occupied with all P(0)es. Thus the determinant equals the
power of the P(0)es!

det(P(k)) = P(0)^k                  //crazy

Now you've got everything to do the deconvolution-in-time. If you have
some math tools available at hand like Maple, Matlab or Mathematica the
software can do the work for you. .. I have a working Maple worksheet
available you can better make use of this. If you like to brew your own
software look at numerical recipes
http://www.library.cornell.edu/nr/cbookcpdf.html

I hope I could help you... Please ask me if you did not understand
this!

To all others: If I did something wrong in my 'assumptions' above -
please correct me...

All the best

Christian






Julian Stoev wrote:

> Hello!
> I have a output signal y(k) and a plant P(z). The signal y(k) contains
> some noise and I know the PSD of the noise. But lets assume that for now
> that the noise is qhite and the system looks like this:
>
> Y(z)=U(z)P(z)+E(z)
>
> u(k)                      y(k)
> ---------->[P(z)]-->[+]--->
>                       |
>                 e(k)  |
>                -------|
>
> I want to restore the signal u(k) from y(k). I can do this off-line in
> freq domain, no problem, but I need this on-line. I think I need
> something related to Wiener filters.
>
> The deconvolution P^(-1)(z) may be unstable, which is the biggest
> problem for me. I am not worried that it may be no casual.
>
> I read that there are Kalman filter based approaches to this, but I
> can't find some practical explanation I can understand in reasonable
> amount of time.
>
> Can somebody give me some reading directions - books, papers. Better
> papers, because I can download most of them in PDF.
> 
> Thanks a lot!
> 
> 
> --JS

"Julian Stoev" <stoev@deadspam.com> wrote in message
news:436afce1$0$41145$14726298@news.sunsite.dk...
> Hello!
> I have a output signal y(k) and a plant P(z). The signal y(k) contains
> some noise and I know the PSD of the noise. But lets assume that for now
> that the noise is qhite and the system looks like this:
>

Should that read 'the noise is shite'?

McC

"Julian Stoev" <stoev@deadspam.com> wrote in message
news:436afce1$0$41145$14726298@news.sunsite.dk...
> Hello!
> I have a output signal y(k) and a plant P(z). The signal y(k) contains
> some noise and I know the PSD of the noise. But lets assume that for now
> that the noise is qhite and the system looks like this:
>
> Y(z)=U(z)P(z)+E(z)
>
> u(k)                      y(k)
> ---------->[P(z)]-->[+]--->
>                       |
>                 e(k)  |
>                -------|
>
> I want to restore the signal u(k) from y(k). I can do this off-line in
> freq domain, no problem, but I need this on-line. I think I need
> something related to Wiener filters.
>
> The deconvolution P^(-1)(z) may be unstable, which is the biggest
> problem for me. I am not worried that it may be no casual.
>
> I read that there are Kalman filter based approaches to this, but I
> can't find some practical explanation I can understand in reasonable
> amount of time.
>
> Can somebody give me some reading directions - books, papers. Better
> papers, because I can download most of them in PDF.
>
> Thanks a lot!
>
>
> --JS
Do you know u(k) for some short period (as in equalisation) or is this blind
deconvolution?

McC

Christian Langen wrote:
...
> The matrix inversion is done by multiplying the inverse determinant
> with the matrix itself. Before you scratch your head rethinking this
> sentence here's the formula:
> 
> P(k)^-1 = 1/det(P(k))  P(k)

Are you really sure about this part?

Real_McCoy wrote:

> Do you know u(k) for some short period (as in equalisation) or is this blind
> deconvolution?

No, I don't know u(k). And IMHO, this is not a blind deconvolution. In 
blind deconvolution you don't know P(z) and u(k). In my case I have a 
pretty good idea about P(z), assume that I know P(z) 100%. So this is 
just deconvolution.

BTW, I start now to read the answer by Christian Langen. It looks what I 
may need, but I just start reading it. Thanks a lot!


--JS

Hello Christian,

Your approach is based on the assumption, that P(k) is FIR filter? I did 
not mention it, but I used P(z). ;-)

My case is related to IIR description of the system. For example

      5.641e-006 z^3 + 0.1712 z^2 + 0.8146 z - 0.08121
P(z)=------------------------------------------------
            z^3 - 2.095 z^2 + 1.231 z - 0.1322

It does not look that I can apply your inversion to the IIR case?

Thanks a lot!

--JS



Christian Langen wrote:
> Dear Julian,
> 
> some years before I thought about a similar problem and found a
> solution:
> 
> The convolution can be considered as a multiplication of a vector y(k)
> and a triangular (Toeplitz) matrix that is built up from the plant
> P(k). For a better understanding I try to develop this in formulae.
> 
> y(k) = x(k) * P(k)          // I use * as symbol for convolution
> 
> y(k) = SUM(x(k) P(k-n))     // This is the classic discrete time
> convolution formula: n is the 'running' variable, the sum is built up
> for all n from 0 to k-1.
> 
> the formula above is the said vector/matrix multiplication: y(k) =
> x(k)P(k) - You did correct by using a capital P for the plant to symbol
> a matrix :-)
> 
> y(0)      |P(0)   0    0  ...      0|      x(0)
> y(1)      |P(1)  P(0)  0 ..        0|      x(1)
> y(2)      |P(2)  P(1)  P(0)..      0|      x(2)
> .      =  |P(3)  P(2)  P(1) P(0).. 0|  x     .
> .         |P(4)     ...            0|        .
> .         | .........              0|        .
> y(k-1)    |P(k-1) P(k-2) ...    P(0)|      x(k-1)
> 
> The Matrix (hope it got not scrambled too much ;-)) is very
> straigtforward: the diagonal is occuied with all P(0)es.
> 
> Now you want to extract the vector x(k). The first 'idea' is to
> 'devide' the vector y(k) by the matrix P(k). O.K. this is done by
> matrix inversion of P(k):
> 
> x(k) = y(k) P(k)^-1
> 
> The matrix inversion is done by multiplying the inverse determinant
> with the matrix itself. Before you scratch your head rethinking this
> sentence here's the formula:
> 
> P(k)^-1 = 1/det(P(k))  P(k)
> 
> To have a matrix invertable there have to be some conditions to be
> fulfilled: The determinant of the matrix P(k) must not be zero
> (otherwise the inverse matrix would have all infinities :-)). How to
> determine the determinant? Before you read the entire math library of
> your university: All upward and downward diagonal products that have to
> be summed up and substracted in the crazy manner you have learned at
> High School will result to ZERO since you have the upper right triangle
> of your Toeplitz Matrix occupied with zeroes! The only exception is the
> diagonal occupied with all P(0)es. Thus the determinant equals the
> power of the P(0)es!
> 
> det(P(k)) = P(0)^k                  //crazy
> 
> Now you've got everything to do the deconvolution-in-time. If you have
> some math tools available at hand like Maple, Matlab or Mathematica the
> software can do the work for you. .. I have a working Maple worksheet
> available you can better make use of this. If you like to brew your own
> software look at numerical recipes
> http://www.library.cornell.edu/nr/cbookcpdf.html
> 
> I hope I could help you... Please ask me if you did not understand
> this!
> 
> To all others: If I did something wrong in my 'assumptions' above -
> please correct me...
> 
> All the best
> 
> Christian

"Julian Stoev" <stoev@deadspam.com> wrote in message
news:436e9b60$0$41144$14726298@news.sunsite.dk...
> Real_McCoy wrote:
>
> > Do you know u(k) for some short period (as in equalisation) or is this
blind
> > deconvolution?
>
> No, I don't know u(k). And IMHO, this is not a blind deconvolution. In
> blind deconvolution you don't know P(z) and u(k). In my case I have a
> pretty good idea about P(z), assume that I know P(z) 100%. So this is
> just deconvolution.
>
> BTW, I start now to read the answer by Christian Langen. It looks what I
> may need, but I just start reading it. Thanks a lot!
>
>
> --JS

The answer is a straight forward Weiner filter for Deconvolution.You can
improve on this by doing Smoothing.

McC

You got me: The matrix P(k) has to be mirrowed, i.e. the rows and
columns have to be changed

|P(0)  0   0 ...  0|
|P(1) P(0)   ...  0|
|  ....            |
|P(k-1) ....   P(0)|

has to be substituted by

|P(0) P(1).... P(k-1)|
|0    P(0) ....P(k-2)|
|0     ....          |
|0     ....    P(0)  |

An internal student in our company did the same mistake like me some
years before :-).
Nevertheless the factor 1/det(P(k)) has to be determined as
described...

Regards

Christian

Dear Julian,

sorry my home-brewn deconvolution-in-time algorithm (updated some
minutes before) doesn't work as recursive filter approach. Use a Wiener
or Kalman filter instead if you want this....

All the best

Christian