division of power series and (direct) deconvolution

[Music Engineer]

yesterday night, when reading in the book "Numerical Methods for Scientists and Engineers" by R.W.Hamming, the author majorly perplexed me on page 48. the goal was to evaluate ln(1-x)/ln(1+x) by means of separate power series for numerator and denominator. the expansion is given by:

Code: Select all

 ln(1-x)     -(x + x^2/2 + x^3/3 + x^4/4 + ...)
--------- = ------------------------------------
 ln(1+x)       x - x^2/2 + x^3/3 - x^4/4 + ...

now he says: "Divide out the two series formally to get:"

Code: Select all

 ln(1-x)     
--------- = -(1 + x + x^2/2 + 5x^3/12 + ...)
 ln(1+x) 

and i was just like: WTF?! divide out formally? is that supposed to be obvious? at least for me, it wasn't. now i think, i have figured out what has to be done. i think, it involves a deconvolution of two sequences of numbers. i wrote a little latex document about that, in order to not have to go through it again:

www.rs-met.com/documents/dsp/PowerSeries.pdf

maybe it can be done in much simpler terms and i'm just not seeing the wood for the trees. if so, i'd be happy if someone could point me in the right direction.

anyway - i was rewarded with a nice little direct (time-domain) deconvolution algorithm. i guess, i re-invented the wheel (yet again) - but i'm not aware of any material on direct (i.e. not FFT based) deconvolution.

[Ivan_C]

I think you might find information on Google with "Taylor series division"

I have found something here : http://en.wikibooks.org/wiki/Calculus/Taylor_series

[Ivan_C]

I'm doing the same thing than in the link I have given. So we are looking for these terms :

Code: Select all

 ln(1-x)     
--------- =  c0 + c1 x + c2 x^2 + c3 x^3 + ...
 ln(1+x) 

I can write that also this way :

ln(1-x) = (c0 + c1 x + c2 x^2 + c3 x^3 + ...) ln(1+x)

Or :

ln(1-x) = (c0 + c1 x + c2 x^2 + c3 x^3 + ...) (x - x^2/2 + x^3/3 - x^4/4 + ...)

And :

ln(1-x) = c0 x + (c1 - c0/2) x^2 + (c2 - c1/2 + c0/3) x^3
+ (c3 - c2/2 + c1/3 - c0/4) x^4 + ...

If I identify the terms of this relationship to the ones from the initial Taylor serie of ln(1-x), I can find the solution :

c0 = -1
c1 - c0/2 = -1/2 => c1 = -1
c2 - c1/2 + c0/3 = -1/3 => c2 = -1/2
c3 - c2/2 + c1/3 - c0/4 = -1/4 => c3 = -5/12

Cool I have the same results (I was afraid a little)

Anyway, I have read your PDF. b0 is null in this example right ?

[Ivan_C]

I'm going to try to generalize a little. So I have two polynomials I want to divide :

P1 = b0 + b1 x + b2 x^2 + b3 x^3 + ...
P2 = a0 + a1 x + a2 x^2 + a3 x^3 + ...

I want to calculate P1 / P2 = c0 + c1 x + c2 x^2 + c3 x^3 + ... = P3

Using the previous method :

P1 = c0 a0 + (c1 a0 + c0 a1) x + (c2 a0 + c1 a1 + c0 a2) x^2
+ (c3 a0 + c2 a1 + c1 a2 + c0 a3) x^3 + ...

So, we get the following identities :

b0 = c0 a0
b1 = (c1 a0 + c0 a1)
b2 = (c2 a0 + c1 a1 + c0 a2)
b3 = (c3 a0 + c2 a1 + c1 a2 + c0 a3)
...

At the end, we get a linear equation system to solve. It can be done using matrix representations.

Code: Select all

    | b0 |   | a0 0  0  0  |   | c0 |
B = | b1 | = | a1 a0 0  0  | * | c1 | = M * C
    | b2 |   | a2 a1 a0 0  |   | c2 |
    | b3 |   | a3 a2 a1 a0 |   | c3 |

And finally :

C = M ^ -1 * B

... to be continued (things are not that easy when M is not invertible)

[Music Engineer]

C = M ^ -1 * B

yep - and taking advantage of the triangular structure results in this "deconvolution" recursion

[Music Engineer]

Anyway, I have read your PDF. b0 is null in this example right ?

i actually assume it to be nonzero because i divide by it. if one or more of the low-index coefficients are zero, the whole equation system changes a little - some must be skipped and division occurs by the first nonzero coefficient etc.

in code, it currently looks like that, over here:

Code: Select all

/** Deconvolves the impulse response h out of the signal y 
resulting in the signal x which has a length of 
yLength-hLength+1. */

template <class T>
void rsDeConvolve(T *y, int yLength, T *h, int hLength, T *x)
{
  int m = rsFirstIndexWithNonZeroElement(h, hLength);
  if( m == -1 )
  {
    // h is all zeros - return an all-zero x signal:
    rsFillWithZeros(x, yLength-hLength+1);
    return;
  }
  T scaler = T(1) / h[m];
  x[0]     = scaler * y[m];
  for(int n = 1; n < yLength-hLength+1; n++)
  {
    x[n] = y[n+m];
    for(int k = m+1; k <= rsMin(hLength-1, n+m); k++)
      x[n] -= h[k] * x[n-k+m];
    x[n] *= scaler;
  }
}

i did some unit tests - it seems to work.

[Ivan_C]

division occurs by the first nonzero coefficient etc.

Makes sense

[Music Engineer]

deleted - wanted to edit a post but quoted myself instead. (happens every once in a while)

[xoxos]

and i was just like: WTF?! divide out formally? is that supposed to be obvious? at least for me, it wasn't.

[kryptonaut]

You can also do it by laying it out in a similar way to a traditional long-division sum, only using the coefficients of powers of x instead of digits (which are coefficients of powers of 10).

[Music Engineer]

You can also do it by laying it out in a similar way to a traditional long-division sum, only using the coefficients of powers of x instead of digits (which are coefficients of powers of 10).

yes - apparently, if polynomial multiplication is done by convolving the coefficient arrays, polynomial long division (which generalizes numeric long division) ought to be carried out by deconvolution. it's just that, in deconvolution (in the way given in the code above), we ignore the possibility for a remainder.