FFT-based fast polynomial rooting - Semantic Scholar
FFT-BASED FAST POLYNOMIAL ROOTING Leila Hoteit
Communication and Signal Processing section, Imperial College, London SW7 2BT, England, 1.hoteit 0 ic.ac.uk ABSTRACT A fast, accurate and robust approach is proposed for the computation of the roots of complex polynomials. The method is derived from the DFT-based differential cepstrum estimation for moving average signals. This minimudmaximum-phase polynomial factorisation is easily extended to a factorisation along an arbitrary circle. In an iterative fashion, we estimate the largest root modulus from the differential cepstrum, then factor out the associated root(s) from the polynomial. For band-limited signals with roots located along the unit circle, the polynomial origin is slightly perturbed prior to the application of the algorithm. On average, three Fast Fourier Transforms are required per polynomial root, offering a significant computational advantage in the root estimation of moderate to high order polynomials.
1. INTRODUCTION Numerous polynomial root finding algorithms can be found in the literature ( [ 6 ] ) :the eigenvalue-based method, the Jenkins-Traub approach, the Lehmer-Schur technique, etc. Although robust in practical terms, the QR method, used in the root finding technique for computing the eigenvalues of a companion matrix, is computationally demanding (order N 3 , where N is the polynomial degree) and is not guaranteed to converge to a solution. The JenkinsTraub approach yields very accurate approximations to the roots, but usually fails to find the zeros of a high degree polynomial. The Lehmer-Schur technique is a root detection algorithm, and can efficiently determine whether there are any polynomial roots within a circle of given centre and radius. Used as a root finding technique, it trades off high accuracy for the high degree of the polynomial considered. There exists a need for a root finding technique that is fast, accurate, and robust to high order, possibly ill-conditioned complex polynomials. The use of root finding techniques in the estimation of the cepstra of a short moving average sequence is suggested in [SI. Turning matters around, we propose the use of the FIT-based differential cepstrum estimation for the purpose of computing the roots of a moderate to high order complex polynomial. The differential cepstrum provides a minimum phase factorisation of the polynomial. We show how a polynomial can be factorised across an arbitrary circle through the simple transformations of translation and dilation. We then show how the largest root modulus can be accurately estimated from the differential cepstrum. The root finding method consists of an iterative estimation
of the largest magnitude, followed by an estimation of the associated root(s) through an FIT-based factorisation. The computational cost is approximately three Fast Fourier Transforms per root of a complex polynomial. The accuracy will depend on the chosen length of the FFT, and the potential presence of root clusters. The paper is organised as follows: in sections 2 and 3, we briefly review the complex and differential cepstra, and the FFTbased estimation of the cepstrum. In section 4,we show how the largest root modulus can be accurately estimated. This is followed by a description of the algorithm in section 5, and simulation results in section 6.
2. THE CEPSTRA AND THE ROOT SUMS 2.1. The Complex Cepstrum Consider a sequence 41 . whose z-transform expressed in polar form is: x ( Z ) = I X ( Z ) ~eJLX(z).The complex cepstrum a[.] of+] is defined as the inverse z-transform of:
The complex cepstrum exists if l o g [ X ( z ) ]has a convergent power series representation of the form: C:z--i[n]z-"; IzI = 1. The z-transform X ( z ) of a Moving Average (MA) sequence 41 . can be normalised to the form:
k= 1
(2)
k= I
where lakl and lbkl are both less than unity and correspond respectively to the roots inside the unit circle, and the reciprocal of the roots outside the unit circle. We assume for the moment that there are no roots on the unit circle.
2.2. The Differential Cepstrum The differential cepstrum was introduced by Polydoros and Fam in [ 5 ] , and is defined by:
{Z}
4 [ n ]4 z-1 {&(z)} 4 z-'
(3)
The differential cepstrum satisfies the difference equation:
+-(m
The author is now with Schlumberger Cambridge Research, CB3 OEL, England
0-7803-6293-4/00/$10.0002000 IEEE.
MO
M,
X(Z)=IAIn(I-akZ-'),(,-DkZ):
-
I ) . x [ m- 1 1 =
d,[k]x[m- k ] . k=--
3315
(4)
The relationship between the differentia! and the complex cepstrum, when the latter exists, is simple: X ' ( z ) = D,(z), yet the former has a well defined phase. Thus, by using the differential cepstrum, we avoid the phase ambiguity associated with the logarithm. Taking the inverse z-transform, when .?[n]exists, we get: d,[m] = - ( m - 1)i[m- 11.
Thus, the S-parameters defined in [7] are an alternative interpretation to the differential cepstrum. The interpretation is based upon the Newton Identities and their examination from a complex variable perspective.
(5)
Using relationship (9,we see that the differential cepstrum of a normalised MA sequence is given by:
3. POLYNOMIAL FACTORISATION USING THE FFT In [4], an implementation based on the Discrete Fourier Transform (DFT)of the differential cepstrum is proposed. This approach offers a minimudmaximum phase factorisation of an MA sequence without knowledge of its roots. For any integer m: If the z-transform of the Moving Average sequence is not normalised, and/or includes a phase shift, its differential cepstrum is still defined and can be calculated based on the following properties: e e
Scale invariance: y[n]= A . x [ n ] , A E C +4 . [ n ]= d , [ n ] .
Thus, when there are no roots on the unit circle, the contour integral for 5 r can be approximated by the summation:
Shift invariance and time delay measurement: y[n]= x [ n + r ] : r E Z : (7)
2.3. The S-parameters The S-parameter is a quantity defined by Stathaki [7] for polynomial sequences, and is equal to the mrh sum of the roots of the polynomial. In terms of ak and bk in equation 2, we have, for m E Z:
(14)
(8)
where the DFT and IDFT are respectively the Discrete Fourier Transform and its inverse. The Discrete Fourier Transforms have been zero-padded to N F , a power of 2, enabling the use of the computationally efficient FFT.
where S:'" and 5Frare the mrh sum of the zeros inside and outside the unit circle respectively. Thus the relationship to the differential cepstrum of the sequence is straightforward:
By taking the integral of P around different contours, F (4 we can perform various polynomial factorisation. For example, : of the roots
Communication and Signal Processing section, Imperial College, London SW7 2BT, England, 1.hoteit 0 ic.ac.uk ABSTRACT A fast, accurate and robust approach is proposed for the computation of the roots of complex polynomials. The method is derived from the DFT-based differential cepstrum estimation for moving average signals. This minimudmaximum-phase polynomial factorisation is easily extended to a factorisation along an arbitrary circle. In an iterative fashion, we estimate the largest root modulus from the differential cepstrum, then factor out the associated root(s) from the polynomial. For band-limited signals with roots located along the unit circle, the polynomial origin is slightly perturbed prior to the application of the algorithm. On average, three Fast Fourier Transforms are required per polynomial root, offering a significant computational advantage in the root estimation of moderate to high order polynomials.
1. INTRODUCTION Numerous polynomial root finding algorithms can be found in the literature ( [ 6 ] ) :the eigenvalue-based method, the Jenkins-Traub approach, the Lehmer-Schur technique, etc. Although robust in practical terms, the QR method, used in the root finding technique for computing the eigenvalues of a companion matrix, is computationally demanding (order N 3 , where N is the polynomial degree) and is not guaranteed to converge to a solution. The JenkinsTraub approach yields very accurate approximations to the roots, but usually fails to find the zeros of a high degree polynomial. The Lehmer-Schur technique is a root detection algorithm, and can efficiently determine whether there are any polynomial roots within a circle of given centre and radius. Used as a root finding technique, it trades off high accuracy for the high degree of the polynomial considered. There exists a need for a root finding technique that is fast, accurate, and robust to high order, possibly ill-conditioned complex polynomials. The use of root finding techniques in the estimation of the cepstra of a short moving average sequence is suggested in [SI. Turning matters around, we propose the use of the FIT-based differential cepstrum estimation for the purpose of computing the roots of a moderate to high order complex polynomial. The differential cepstrum provides a minimum phase factorisation of the polynomial. We show how a polynomial can be factorised across an arbitrary circle through the simple transformations of translation and dilation. We then show how the largest root modulus can be accurately estimated from the differential cepstrum. The root finding method consists of an iterative estimation
of the largest magnitude, followed by an estimation of the associated root(s) through an FIT-based factorisation. The computational cost is approximately three Fast Fourier Transforms per root of a complex polynomial. The accuracy will depend on the chosen length of the FFT, and the potential presence of root clusters. The paper is organised as follows: in sections 2 and 3, we briefly review the complex and differential cepstra, and the FFTbased estimation of the cepstrum. In section 4,we show how the largest root modulus can be accurately estimated. This is followed by a description of the algorithm in section 5, and simulation results in section 6.
2. THE CEPSTRA AND THE ROOT SUMS 2.1. The Complex Cepstrum Consider a sequence 41 . whose z-transform expressed in polar form is: x ( Z ) = I X ( Z ) ~eJLX(z).The complex cepstrum a[.] of+] is defined as the inverse z-transform of:
The complex cepstrum exists if l o g [ X ( z ) ]has a convergent power series representation of the form: C:z--i[n]z-"; IzI = 1. The z-transform X ( z ) of a Moving Average (MA) sequence 41 . can be normalised to the form:
k= 1
(2)
k= I
where lakl and lbkl are both less than unity and correspond respectively to the roots inside the unit circle, and the reciprocal of the roots outside the unit circle. We assume for the moment that there are no roots on the unit circle.
2.2. The Differential Cepstrum The differential cepstrum was introduced by Polydoros and Fam in [ 5 ] , and is defined by:
{Z}
4 [ n ]4 z-1 {&(z)} 4 z-'
(3)
The differential cepstrum satisfies the difference equation:
+-(m
The author is now with Schlumberger Cambridge Research, CB3 OEL, England
0-7803-6293-4/00/$10.0002000 IEEE.
MO
M,
X(Z)=IAIn(I-akZ-'),(,-DkZ):
-
I ) . x [ m- 1 1 =
d,[k]x[m- k ] . k=--
3315
(4)
The relationship between the differentia! and the complex cepstrum, when the latter exists, is simple: X ' ( z ) = D,(z), yet the former has a well defined phase. Thus, by using the differential cepstrum, we avoid the phase ambiguity associated with the logarithm. Taking the inverse z-transform, when .?[n]exists, we get: d,[m] = - ( m - 1)i[m- 11.
Thus, the S-parameters defined in [7] are an alternative interpretation to the differential cepstrum. The interpretation is based upon the Newton Identities and their examination from a complex variable perspective.
(5)
Using relationship (9,we see that the differential cepstrum of a normalised MA sequence is given by:
3. POLYNOMIAL FACTORISATION USING THE FFT In [4], an implementation based on the Discrete Fourier Transform (DFT)of the differential cepstrum is proposed. This approach offers a minimudmaximum phase factorisation of an MA sequence without knowledge of its roots. For any integer m: If the z-transform of the Moving Average sequence is not normalised, and/or includes a phase shift, its differential cepstrum is still defined and can be calculated based on the following properties: e e
Scale invariance: y[n]= A . x [ n ] , A E C +4 . [ n ]= d , [ n ] .
Thus, when there are no roots on the unit circle, the contour integral for 5 r can be approximated by the summation:
Shift invariance and time delay measurement: y[n]= x [ n + r ] : r E Z : (7)
2.3. The S-parameters The S-parameter is a quantity defined by Stathaki [7] for polynomial sequences, and is equal to the mrh sum of the roots of the polynomial. In terms of ak and bk in equation 2, we have, for m E Z:
(14)
(8)
where the DFT and IDFT are respectively the Discrete Fourier Transform and its inverse. The Discrete Fourier Transforms have been zero-padded to N F , a power of 2, enabling the use of the computationally efficient FFT.
where S:'" and 5Frare the mrh sum of the zeros inside and outside the unit circle respectively. Thus the relationship to the differential cepstrum of the sequence is straightforward:
By taking the integral of P around different contours, F (4 we can perform various polynomial factorisation. For example, : of the roots
