room impulse response analysis and polynomial factoring

i am doing some research on acoustics and particularly the nature of
room reverberation. this led me to some harsh number crunching using
matlab.
i am working with polynomials having hundreds of terms and very high
powers in the range of 800 to 2000.
i am particularly interested in the position of the roots relative to
the unit circle. Thus i am using roots() a matlab funtion to compute
the roots of the polynomial. This seems to be working fine (altough i
presume i have rounding errors).
Now my intention is to position all roots inside the unti circle by
computing the reciprocal of the roots outside and then recompute the
coefficients of the new polynomial using the poly() function
This however is leading to a false result. The functions roots() and
poly(), seem to work well when i tested them with small polynomials as
shown below:

polycoeff =

    1.0000    0.6000    0.4000    0.2000    0.1000

>> polyroots=roots(polycoeff)

polyroots =

   0.1340 + 0.5544i
   0.1340 - 0.5544i
  -0.4340 + 0.3450i
  -0.4340 - 0.3450i

>> originalpoly=poly(polyroots)

originalpoly =

    1.0000    0.6000    0.4000    0.2000    0.1000

>>

Is this a problem of rounding? I even tried using some other code i
found on matlab's file exchange poly1() and roots1() but still same
false result occur. The coefficients of the resulting polynomial range
from very small 10^-40 to 10^128 which is not the case with my
polynomials which have the largest coeff at 1. Any help greatly
appreciated.

jeremy

jeremyscerri@gmail.com wrote:
> i am doing some research on acoustics and particularly the nature of
> room reverberation. this led me to some harsh number crunching using
> matlab.
> i am working with polynomials having hundreds of terms and very high
> powers in the range of 800 to 2000.

Ouch. You are pushing the limits here, no doubt about that.
You would need to look at very specialized computation routines.
I have seen mentioned a paper on an algorithm for rooting a high
order (> 100) polynomial, but I can't remember where to find it.

> i am particularly interested in the position of the roots relative to
> the unit circle. Thus i am using roots() a matlab funtion to compute
> the roots of the polynomial. This seems to be working fine (altough i
> presume i have rounding errors).
> Now my intention is to position all roots inside the unti circle by
> computing the reciprocal of the roots outside and then recompute the
> coefficients of the new polynomial using the poly() function
> This however is leading to a false result. The functions roots() and
> poly(), seem to work well when i tested them with small polynomials as
> shown below:
>
> polycoeff =
>
>     1.0000    0.6000    0.4000    0.2000    0.1000
>
> >> polyroots=roots(polycoeff)
>
> polyroots =
>
>    0.1340 + 0.5544i
>    0.1340 - 0.5544i
>   -0.4340 + 0.3450i
>   -0.4340 - 0.3450i
>
> >> originalpoly=poly(polyroots)
>
> originalpoly =
>
>     1.0000    0.6000    0.4000    0.2000    0.1000
>
> >>
>
> Is this a problem of rounding? I even tried using some other code i
> found on matlab's file exchange poly1() and roots1() but still same
> false result occur. The coefficients of the resulting polynomial range
> from very small 10^-40 to 10^128 which is not the case with my
> polynomials which have the largest coeff at 1. Any help greatly
> appreciated.

Rounding errors are a big issue. Matlab uses an eigenvalue
decomposition for finding the roots, which is a O(n^3) algorithm.
So basically, in a polynomial of degree 2000 you run through on
the order of 2000^3 -- almost ten billion -- FLOPS. If you assume
you accumulate 1 eps numerical error per FLOP, there is
almost no accuracy left when the routine finishes.

The bottom line is that your application is really pushing the limits.
It is, at best, very difficult to achieve what you want.

Rune

Rune Allnor wrote:
> jeremyscerri@gmail.com wrote:

>> i am doing some research on acoustics and particularly the
>> nature of room reverberation. this led me to some harsh number
>> crunching using matlab.
>> i am working with polynomials having hundreds of terms and very
>> high powers in the range of 800 to 2000.
> 
> Ouch. You are pushing the limits here, no doubt about that.
> You would need to look at very specialized computation routines.
> I have seen mentioned a paper on an algorithm for rooting a high
> order (> 100) polynomial, but I can't remember where to find it.

Here's a paper that talks about degrees in the hundred thousands:
http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=1253552


Martin

-- 
Quidquid latine scriptum sit, altum viditur.

On 5 Mar 2006 12:52:29 -0800, "jeremyscerri@gmail.com"
<jeremyscerri@gmail.com> wrote:

>i am doing some research on acoustics and particularly the nature of
>room reverberation. this led me to some harsh number crunching using
>matlab.
>i am working with polynomials having hundreds of terms and very high
>powers in the range of 800 to 2000.
  (snipped)
>
>Is this a problem of rounding? I even tried using some other code i
>found on matlab's file exchange poly1() and roots1() but still same
>false result occur. The coefficients of the resulting polynomial range
>from very small 10^-40 to 10^128 which is not the case with my
>polynomials which have the largest coeff at 1. Any help greatly
>appreciated.
>
>jeremy

Hi Jeremy,

  Ya' might check out the November 2003
issue of the IEEE Signal Processing magazine.
There's an article there discussing "factoring 
polynomials" whose orders are in the 
hundreds of thousands!

Good Luck,
[-Rick-]
 

thank you very much for your input, especially the one related to
matlab eigenvalue decomposition routine, will check out the IEEE paper.
jeremy

jeremyscerri@gmail.com wrote:
> thank you very much for your input, especially the one related to
> matlab eigenvalue decomposition routine, will check out the IEEE paper.
> jeremy
> 

Also:

   http://www.dsp.rice.edu/software/fvhdp.shtml


Bob
-- 

"Things should be described as simply as possible, but no simpler."

                                              A. Einstein

thanks bob, for your very valid input. will try out the code. I have
also been directed to this other site
http://ccrma.stanford.edu/~jos/filters/Creating_Minimum_Phase_Filters.html

apparently this complex cepstrum thing avoids polynomial factoring
altogether and still serves my purpose of moving all outside zeros
within the unit circle.
 i have already tried it and it does what it says, the magnitude
spectrum is unchanged to the expense of changing the phase. I have also
convolved this new minimum phase impulse response with some anechoic
sources and I couldn't perceive anything different when compared to the
convolution of the original mixed phase IR with same sources. Does
anyone know any references to the human perception and phase changes, i
once read something about it but lost the links.

thanks
jeremy

jeremyscerri@gmail.com wrote:
> thanks bob, for your very valid input. will try out the code. I have
> also been directed to this other site
> http://ccrma.stanford.edu/~jos/filters/Creating_Minimum_Phase_Filters.html
> 
> apparently this complex cepstrum thing avoids polynomial factoring
> altogether and still serves my purpose of moving all outside zeros
> within the unit circle.

Right!  If that is your sole purpose for root finding, then by all means 
avoid it with this superb method.  A Matlab function I wrote several 
years ago to implement this follows my sig.

>  i have already tried it and it does what it says, the magnitude
> spectrum is unchanged to the expense of changing the phase. I have also
> convolved this new minimum phase impulse response with some anechoic
> sources and I couldn't perceive anything different when compared to the
> convolution of the original mixed phase IR with same sources. Does
> anyone know any references to the human perception and phase changes, i
> once read something about it but lost the links.

It is very murky and undecided.  If you were dealing with HRTF's you 
would lose all time of arrival information and be left with only 
relative amplitude.  It would be an interesting experiment to see if 
such HRTF's lost the ability to simulate spatial positioning.  If said 
HRTF's also included room acoustics then I would expect considerable 
loss because of the IR length compression.

If you are simply seeking magnitude response contouring then I think 
minimum phase is the best way to go.  I've done a fair bit of comparing 
minimum with linear phase equalizing and while the difference is subtle 
I definitely prefer minimum phase.


Bob

*-----------------------------------------------
function y = minphase(x)
%minphase - minimum phase from mixed phase impulse responses
%
%   y = minphase(x) returns the minimum phase signals derived from x.
%
%   x can be a column array of impulse responses in which y will be
%   a column array of the minimum phase equivalents.
%
%   References:
%		http://www.nauticom.net/www/jdtaft/minphase.htm
%
% Author Bob Cain arcane@arcanemethods.com

%Double the impulse response lengths with zeros
	x = [x; zeros(size(x))];
%Number of impulse responses
	m = size(x,2);
%Size of each
	n = size(x,1);
%The weighting function
	w = [ones(n/2,m); -ones(n/2,m)];
%The log magnitude
	lm=log(abs(fft(x)));
%The work
	y = ifft(exp(lm+i*imag(fft(w.*ifft(lm)))));
%Return half the result
	y = y(1:n/2,:);
*-----------------------------------------------
-- 

"Things should be described as simply as possible, but no simpler."

                                              A. Einstein