Re: Roots() and factorizing large polynomials

[Dupuis]

Borge wrote:
> 
> roots_x = multroot(x); % x is impulse response of digital FIR filter
> manipulate(roots_x); % Invert roots with ||>1 to obtain minimum-phase
> roots
> 
> % Reassemble new impulse response from roots:
> asm_x = [1];
> for n=1:length(roots_x)
>         for m=1:roots_x(n,2)
>             asm_x = conv(asm_x, [1 -roots_x(n,1)]);
>         end;
> end
> 
> This reassembly method is very crude, since it requires as many
> convolutions as there are roots. The output grows to a very large
> number, and it is not real. So numerical precision must play a large
> role here. I can reduce the number by using iterative, branched
> convolutions. Or would you suggest an alternative method?
> 

Your code seems right. You have indeed to convolve all the roots together.
Maybe you could do
1) extract the roots using multroot()
2) reflect the ones outside the unit circle inside it (i.e. stabilise them)
3) verify afterwards that the new set of roots extracted from a new call to
multroot() match the roots you put at stage 2

Now, if you have complex-conjugate pairs, you have to be very carefull with
them. The convolution should proceed as
- if the x(i) is real, convolves m times with [1 -x(i)]
- if the x(i) is complex, do m times: asm_x = conv(asm_x, [1 -2*real(xi) 
(real(x(i))^2+imag(x(i))^2])
This way, at each stage you use only pure real polynomials. If x(i) is
complex, do not apply the previous formula a second time on x(i) complex
conjugate.

This approach avoids the builtup of small errrors, which will causes a
imaginary part to appear.

Regards

Pascal
-- 
View this message in context: 
http://old.nabble.com/Roots%28%29-and-factorizing-large-polynomials-tp27502322p27591166.html
Sent from the Octave - General mailing list archive at Nabble.com.