Two new algorithms are presented in the paper for the plus/minus factorization of a scalar discrete-time polynomial. The first method is based on the discrete Fourier transform theory (DFT) and its relationship to the Z-transform. Involving DFT computational techniques and the famous fast Fourier transform routine brings high computational efficiency and reliability. The method is applied in the case study of H 2optimal inverse dynamic filter to an audio equipment. The second numerical procedure originates in a symmetric spectral factorization routine, namely the Bauer's method of the 1950s. As a by-product, a recursive LU factorization procedure for Toeplitz matrices is devised that is of more general impact and can be of use in other areas of applied mathematics as well. Performance of the method is demonstrated by an l 1 optimal controller design example.where p þ ðzÞ has all roots inside and p À ðzÞ outside the unit disc. Clearly, the scalar plus/minus factorization is unique up to a scaling factor.Polynomial plus/minus factorization has many applications in control and signal processing problems. For instance, efficient algebraic design methods for time-optimal controllers [1], quadratically optimal filters for mobile phones [2,3], and l 1 optimal regulators [4], to name just a few, all recall the þ=À factorization as a crucial computational step.
EXISTING METHODSIn any case, the plus/minus factors for n55 cannot be achieved by a finite number of algebraic operations. This conclusion is due to the Galois's theorem stating that the roots of a polynomial of degree greater or equal to five cannot be expressed in a closed form. Therefore, all numerical algorithms for plus/minus factorization are iterative in nature and give just an approximation to the genuine factors. Some existing approaches to this problem are mentioned in this section.The most natural way is based on the computation of polynomial roots. Having determined the roots r 1 ; r 2 ; . . . ; r n of pðzÞ via any standard procedure for polynomial roots [5] and considering that pðzÞ=0 for all jzj ¼ 1 by assumption, one can divide the roots into two groups R þ ¼ fr i : m $ ðr i Þ ¼ 0; jr i j51g; R À ¼ fr i : m $ ðr i Þ ¼ 0; jr i j > 1g: Clearly, R þ and R À are the sets of roots of p þ ðzÞ and p À ðzÞ; respectively.Performance of this procedure heavily hinges on the accuracy of the computed polynomial roots. If these roots are distinct and separated enough, standard numerical routines [5] can determine them with good precision. However, it is well known that the relative accuracy of a computed root decreases as its multiplicity grows [5], and so does the accuracy of the spectral factor thus obtained.In addition, if the degree of the involved polynomial is high, say over 50, the very computation of the spectral factor coefficients is problematic due to rounding errors. It means that even if the desired roots of the spectral factor are evaluated with good accuracy, its particular coefficients, which are typically required in applications, are not accurate...