VIEWS: 10 PAGES: 6 POSTED ON: 11/23/2011
FFT BASED ALGORITHM FOR POLYNOMIAL PLUS-MINUS FACTORIZATION ˇ Martin Hromˇ´k, Michael Sebek cı Centre for Applied Cybernetics Czech Technical University, Prague, Czech Republic e-mail: m.hromcik@c-a-k.cz, m.sebek@c-a-k.cz Keywords: Polynomial design methods, numerical algo- One of them is based on direct computation of roots. Using rithms. standard methods for polynomial roots evaluation, see [5, 18] for instance, one can separate the stable and unstable roots of Abstract Ô´×µ directly and construct the plus and minus parts from re- lated ﬁrst order factors or, alternatively, emply a more efﬁcient In this report a new algorithm is presented for the plus/minus recursive procedure based on the matrix eigenvalue theory [18]. factorization of a scalar discrete-time polynomial. The method Alternative algorithm relies on polynomial spectral factoriza- is based on the discrete Fourier transform theory (DFT) and tion and gratest polynomial divisor computation. If Õ ´Þ µ is its relationship to the Z-transform. Involving DFT computa- the spectral factor of the symmetric product Ô´Þ µÔ´Þ ½ µ then tional techniques, namely the famous fast Fourier transform the greatest common divisor of Ô´Þ µ and Õ ´Þ µ is obviously the routine (FFT), brings high computational efﬁciency and reli- plus factor of Ô´Þ µ. The minus factor can be derived similarly ability. The effectiveness of the proposed algorithm is demon- from Ô´Þ ½ µ and Õ ´Þ ½ µ. As opposed to the previous approach strated by a particular practical application. Namely the prob- based on direct roots computation which typically makes prob- lem of computing an À ¾ -optimal inverse dynamic ﬁlter to an lems for higher degrees and/or roots multiplicities, this proce- audio equipment is considered as it was proposed by M. Ster- dure relies on numerically reliable algorithms for polynomial nad and colleagues in [16] to improve behavior of moderate spectral factorization [11, 2]. Unfortunately, the polynomial quality loudspeakers. Involved spectral factorization can be greatest common divisor computation is much more sensitive. converted into plus-minus factorization in a special case which As a result, both these techniques do not work properly for high in turn is resolved by our new method. degrees (say over 50). In this report we will introduce a completely new approach to 1 Introduction the problem, inspired by our work on efﬁcient algorithms for This paper describes a new method for the plus-minus factor- polynomial spectral factorization, see [2]. It is based on the ization of a discrete-time polynomial. Given a polynomial in DFT theory and provides both a fruitful view on the relation the Þ variable, between DFT and the -transform theory, and a powerful com- putational tool in the form of the fast Fourier transform algo- Ô´Þ µ Ô¼ · Ô½ Þ · Ô¾ Þ ¾ · ¡ ¡ ¡ · ÔÒ Þ Ò rithm. without any roots on the unit circle, its plus/minus factorization 3 Discrete Fourier Transform is deﬁned as Ô´Þ µ Ô· ´Þ µÔ ´Þ µ If Ô Ô¼ Ô½ ÔÆ is a vector of complex numbers, then its Ý (1) direct DFT is given by the vector Ý ¼ Ý½ ÝÆ , where where Ô· ´Þ µ has all roots inside and Ô ´Þ µ outside the unit Æ ¾ Ô Æ ·½ disc. Clearly, the scalar plus/minus factorization is unique up to a scaling factor. Ý (2) ¼ Polynomial plus/minus factorization has many applications in The vector Ý is called the image of vector Ô. Conversely, if control and signal processing problems. For instance, efﬁ- Ý Ý¼ Ý½ Ý is given, then its inverse DFT recovers the original vector Ô cient algebraic design methods for time-optimal controllers [7], Æ quadratically optimal ﬁlters for mobile phones [14, 15], and Ð ½ Ô Ô ¼ ½Ô , where Æ optimal regulators [3], to name just a few, all recall the +/- fac- Æ ½ ¾ Æ ·½ torization as a crucial computational step. Ô Ý (3) Æ ·½ ¼ 2 Existing Methods DFT is of great interest in various engineering ﬁelds. For its From the computational point of view, nevertheless, the task is relationship to Fourier series of sampled signals, DFT is fre- not well treated. There are two quite natural methods. quently used in signal processing. One of the experimental identiﬁcation methods employs DFT as well [9]. The close re- analytic for ½ Þ and ½ · Þ respectively. lationship of DFT to interpolation is also well known and was At this time the necessity of the degree shift yielding the used recently to solve some tasks of the polynomial control two-sided polynomial Ô can be explained. According to the theory [12] and to treat robustness analysis problems of certain Cauchy’s theorem of argument [6], the curve Ô´Þ µ for Þ ½ kind [13]. encircles the origin in the complex plane as many times as is the For numerical computation of DFT, the efﬁcient recursive FFT number of roots of Ô´Þ µ lying in the complex unit disc. Hence algorithm was developed by Cooley and Tukey in 1965 [4]. If the logarithms cannot be applied directly as its imaginary part, the length of the input is a power of two, a faster version of reading the phase of Ô´Þ µ, would not be continuous. An easy FFT (sometimes called radix-2 FFT) can be employed [4]. In solution to avoid this situation is to move the desired number of general, the FFT routine features a highly beneﬁcial compu- roots of Ô´Þ µ from inﬁnity to zero by performing proper degree tational complexity and involves Ç´Æ ÐÓ ´Æ µµ multiplications shift. and additions for a vector of length Æ . Once Ü· ´Þ µ and Ü ´Þ µ are computed, the plus/minus factors Thanks to the importance of DFT mentioned above, the FFT Ô· Ô are recovered as Ü· ´Þ µ Ü ´Þ µ algorithms are naturally available as built-in functions of many computing packages (M ATLAB ÌÅ , M ATHEMATICA ÌÅ etc.). Ô· Ô· ·Ô· Þ ½ ·¡ ¡ ¡ Ô ¼ ½ Ô ·Ô Þ ·¡ ¡ ¡ ¼ ½ This is another good reason for employing the procedure pro- posed in this paper. Since Ü· ´Þ µ is analytic in ½ Þ , so is Ô· ´Þ µ and hence it can be expanded according to (3). Moreover, as a result of 4 Plus-minus Factorization and DFT exponential function, Ô · ´Þ µ is nonzero in ½ Þ . In other words, it has all its zeros inside the unit disc and is therefore 4.1 Theory Schur stable. Note also that Ô · ´Þ µ has to be a (ﬁnite) polyno- mial of degree (due to the uniqueness of the solution to the Given a polynomial problem which is known to be a polynomial) though Ò´Þ µ is an inﬁnite power series. Similar reasoning proves the Ô factor Ô´Þ µ Ô¼ · Ô½ Þ · ¡ ¡ ¡ · Ô Þ desired properties. nonzero for Þ ½, we ﬁrst apply a direct degree shift to arrive at a two-sided polynomial 4.2 Numerical Algorithm Ô´Þ µ Ô¼ Þ Æ · ¡ ¡ ¡ · Ô Þ Æ Numerical implementation follows the ideas considered above. A polynomial Ô´Þ µ is represented by its coefﬁcients Ô where Æ is the number of roots of Ô´Þ µ lying inside the unit ¼ Ö or, equivalently, by function values È in the Fourier in- circle. Now, instead of solving equation (1), we look for terpolating points ¾ Ê ¼ Ê, where Ê Ô· ´Þ µ Ô· · Ô· Þ ½ · ¡ ¡ ¡ · Ô· Þ Æ and Ô ´Þ µ Ô · ¾Ê·½ . Accordingly, a power series can be approximated by a Þ · ¡ ¡ ¡ · Ô Þ Æ such that Æ Ô½ ¼ ½ ¼ Æ ﬁnite set of its coefﬁcients or by its values in a ﬁnite number of interpolation points on the unit circle. Some operations of the Ô´Þ µ Ô· ´Þ µÔ ´Þ µ (4) procedure, namely the decomposition of Ò´Þ µ into Ü · ´Þ µ and Ü ´Þ µ, are performed in the time domain (operations on coef- Relation between the pairs Ô · Ô and Ô· Ô are obvious. ﬁcients), while the others (evaluation of logarithmic and expo- In order to solve the equation (4), logarithm is applied. As nential functions) are executed in the frequency domain (opera- Ô´Þ µ Ô· ´Þ µ and Ô ´Þ µ are all analytic and nonzero in ½ tions with values over Þ ½). Mutual conversion between the Þ ½ · the logarithms exist. Let us denote them as two domains is mediated by the shifted discrete Fourier trans- ÐÒ Ô´Þ µ Ò´Þ µ ÐÒ Ô· ´Þ µ Ü· ´Þ µ ÐÒ Ô ´Þ µ Ü ´Þ µ Here form operator deﬁned as Ò´Þ µ, obtained from the given Ô´Þ µ, is a Laurent inﬁnite power Ê Ê Ü ½ series Ü ¾Ê · ½ Ê Ê Ò´Þ µ ¡¡¡·Ò Þ ·Ò½ ¼ ½ · ¡ ¡ ¡ · Ò ½Þ which approximates the Z-transform by dealing with Ê It can be directly decomposed, ·Ê instead of inﬁnite ½ ·½, and with Þ Ò´Þ µ Ü· ´Þ µ · Ü ´Þ ½ µ Ê ·Ê instead of continuum Þ · . with power series The accuracy of results depends on the number of interpolation Ò¼ points ¾Ê · ½ involved in the computation. This number can Ü· ´Þ µ Ü· · Ü· Þ ½ · ¡ ¡ ¡ ¼ ½ ½ · ¡ ¡ ¡ · Ò ½Þ be considered as a simple tuning knob of the computational ¾ process. Ò¼ Ü ´Þ µ Ü · Ü Þ · ¡ ¡ ¡ ¼ ½ · Ò½ Þ · ¡¡¡ (5) Resulting numerical routine looks then as follows: ¾ Algorithm 1: Scalar discrete-time plus-minus Step 8 - Exponential function: factorization. To get the plus/minus factors, the exponential functions Ô· ´Þ µ Ü· ´Þ µ and Ô ´Þ µ Ü ´Þ µ remain to be evalu- Input: Scalar polynomial ated. First we compute the values of Ô · ´Þ µ and Ô ´Þ µ at Ô´Þ µ Ô¼ · Ô½ Þ · ¡ ¡ ¡ · Ô Þ , nonzero for Þ ½. the Fourier points: È · · ¼ · Ê . Similar steps Polynomials Ô ´Þ µ and Ô ´Þ µ, · apply for the minus part. Output: the plus and minus factors of Ô´Þ µ. Step 9 - Inverse FFT (II): Finally, the coefﬁcients Ô· Ô· ¼ Ô· of Ô· ´Þ µ are Ê Step 1 - Choice of the number of interpolation points. recovered by inverse FFT performed on the vector È . · Decide about the number Ê. Ê approximately ½¼ to ¼ The resulting approximation to the plus factor Ô · ´Þ µ then times larger than is recommended up to our practical equals Ô· ´Þ µ Ô· · Ô·½ Þ ½ · ¡ ¡ ¡ · Ô·Æ Þ Æ Proceed ¼ experience. with the minus part accordingly. Step 2 - Degree shift. Step 10 - Finalization: Find out the number Æ of zeros of Ô´Þ µ inside the unit disc. Convert the plus-minus factors Ô · ´Þ µ and Ô ´Þ µ of Ô´Þ µ A modiﬁcation of well known Schur stability criterion can into the desired factors of Ô´Þ µ using the following formu- be employed, see [10] for instance. las Having Æ at hand, construct a two-sided polynomial Ô´Þ µ Ô Ô Ô· Ô· as where the star stands for discrete-time conjugate, Ô´Þ µ Ô´Þ µÞ Æ Ô¼ Þ Æ · ¡ ¡ ¡ · Ô Þ Æ Þ Þ ½ . ¥ Ô Æ Þ Æ · ¡ ¡ ¡ · Ô¼ · ¡ ¡ ¡ · Ô Æ Þ Æ Note that one obtains Ê coefﬁcients of Ô · and Ô in the step Step 3 - Direct FFT (I): 9. However, Ô · ´Þ µ being the plus factor of Ô´Þ µ is known to Using the FFT algorithm, perform direct DFT, deﬁned by be of degree Æ only and only the ﬁrst Æ · ½ coefﬁcients of (2), on the vector Ô· ´Þ µ should be signiﬁcant as a result while the remaining ones Ô Ô¼ Ô½ Ô Æ ¼ ¼ ßÞ ¼ Ô Æ Ô ½ should be negligible. As the number Ê increases, these values theoretically converge to zero indeed since the formulas of DFT ¾Ê·½ become better approximations to the Z-transform deﬁnitions. In this way, the set È È ¼ È½ È¾Ê of the values of Ô´Þ µ at the Fourier points is obtained. 5 Radix-2 Modiﬁcation of the Algorithm Step 4 - Logarithmization: The basic version of the routine proposed above is based on the Compute the logarithms Æ ÐÒ´È µ of all particular È ’s and form the vector ÆÆ ¼ Æ½ Æ¾Ê of shifted dicrete Fourier transform. This modiﬁcation of DFT appears useful during the derivation of the Algorithm 1 due to them. Æ ’s thus obtained are the values of the function its more transparent relationship to the spectral theory. It can Ò´Þ µ ÐÒ´Ô´Þ µµ at related Fourier points on the unit com- be easily transformed to the standard DFT as it is deﬁned in plex circle. the section 3, simply by reordering related vector entries (see Step 5 - Inverse FFT (I): the steps 2 and 4 of Algorithm 1). However, ¾Ê · ½ interpo- To get the vector ÒÒ ¼ Ò½ ÒÊ Ò Ê Ò ½ , lation points are used for the FFT algorithm and unfortunately this number is always odd and cannot equal any power of two. containing the coefﬁcients of the two-sided polynomial Ò´Þ µ Ò Ê Þ Ê · ¡ ¡ ¡ · Ò ½ Þ ½ · Ò¼ · Ò½ Þ · ¡ ¡ ¡ · ÒÊ Þ Ê Therefore the radix-2 fast version of the FFT routine cannot approximating the power series of ÐÒ´Ñ´Þ µµ for the given be addressed. Nevertheless, this slight drawback can be easily Ê, perform inverse DFT, deﬁned by (3), on the vector Æ avoided if the periodicity of direct and inverse DFT formulas is taken into account. Basically, one can construct the initial set using the FFT algorithm. as Step 6 - Decomposition: Ô¼ Ô½ Ô Æ ¼ ¼ ¼ Ô Æ Ô ½ Take the ”causal part” Ü· of Ò: ßÞ Ü· Ò¼ ¾ Ò½ ÒÊ . Similarly, construct Ü as ¾Ê Ü Ò¼ ¾ Ò ½ Ò Ê . which has a power-of-two entries in total. The Algorithm 1 re- mains valid also in this case with ¾Ê · ½ replaced by ¾ Ê and Step 7 - Direct FFT (II): Ê · ½ by ¾Ê ½ respectively, up to one point: in the Step 6, Evaluate Ü· ´Þ µ Ò¼ ¾ · Ò½ Þ ½ · · ÒÊ Þ Ê at the the decomposition reads · Ü Ò¼ ¾ Ò½ ÒÊ ¾ instead Fourier points by applying direct FFT on the set · and Ü Ü of · Ò¼ ¾ Ò½ ÒÊ . This minor modiﬁcation of the get · Ê . Proceed with Ü ´Þ µ in obvi- · · ¼ proposed method further increases its efﬁciency since the pow- ous way. erful radix-2 FFT can be called. 6 Computational Complexity model ´Þ µ Ý ´Øµ Þ Ù´Øµ Thanks to the fact that the fast Fourier transform algorithm is ´Þ µ extensively used during the computation, the overall routine features an expedient computational complexity. Since the impulse response is rather long for a high sampling frequency (CD-quality standard of 44 kHz was used), both the Provided that the above modiﬁcations of the computational numerator and denominator of the model are of high orders, procedure are considered, namely if the resulting number of say one to ﬁve hundred. interpolation points is taken as a power of two, the fast radix-2 FFT can be employed. In this case, ´Ê ÐÓ ¾ Êµ ¾ multiplica- The model has an unstable inverse in general since some of its tions and Ê ÐÓ ¾ Ê additions are needed to evaluate either di- zeros may lie outside the unit disc. Hence a stable approxima- rect or inverse DFT of a vector of length Ê [4]. Let us suppose tion has to be calculated to be used in the feedforward structure. in addition that computing the logarithm or exponential of a The authors recall the LQG theory and seek for a compensating scalar constant takes at most multiplications and Ð additions. ﬁlter É´Þ µ Then the particular steps of the modiﬁed Algorithm 1 involve Ù´Øµ Û´Øµ ´Ê ÐÓ ¾ Êµ ¾ multiplications and Ê ÐÓ ¾ Ê additions (Steps 3, È ´Þ µ 5, 7, 9), and Ê multiplications and ÐÊ additions (Steps 4, 8) such that the criterion Â Ý ´Øµ Û´Ø µ ¾ · Ù´Øµ ¾ is respectively. Hence the overall procedure consumes minimized. Ê ÐÓ Ê For broadband audio signals, the optimal ﬁlter is given in the ·¾ Ê ¾Ê ÐÓ Ê · ¾ÐÊ ¾ form É½ ´Þ µ ´Þ µ complex multiplications, and Ù´Øµ Û´Øµ ¬ ´Þ µ Ê ÐÓ Ê · ¾ÐÊ where ¬ results from the spectral factorization complex additions. By inspecting the above formulas one ¬¬ £ £· £ can see that asymptotically the proposed method features Ç´Ê ÐÓ Êµ complex multiplications and additions. and É½ is the solution of a subsequent Diophantine equation Þ £ ´Þ µ Ö¬ £ ´Þ µÉ½ ´Þ ½ µ · ÞÄ£ ´Þ µ 7 Upgrading Loudspeakers Dynamics An original approach has been published by Sternad et al. in see [16]. [16] how to improve performance of an audio equipment at ÓÑÔ Ò× ØÓÖ ÅÓ Ð Ó ÐÓÙ ×Ô Ö low additional costs. The authors use the LQG optimal feed- forward compensator technique to receive an inverse dynamic Û´Øµ¹ É´Õ ½ µ È ´Õ ½ µ Ù´Øµ ¹ Õ ´Õ ½ µ ´Õ ½ µ Ý ´Øµ ¹ · ﬁlter for a moderate quality loudspeaker. By attaching a signal processor implementing this ﬁlter prior to the loudspeaker, the ÐÝ ¦ ¹ Ø ´ µ dynamical imperfections of the original device are eliminated and the overall equipment behaves as an aparatus of a much ¹ Õ higher class. To learn more about this research and to get some working examples, visit [17]. Figure 2: Optimal ﬁltering problem setup (adopted from [16]). As for the spectral factor computation, the authors employ the Newton-Raphson iterative scheme [11] in the cited work [16]. Inverse system Hi-fi system According to their results and our experience, this method has been probably the best available procedure for scalar polyno- mial spectral factorization so far [16, 8]. This method works Figure 1: Pre-ﬁlter compensation scheme (adopted from [16]). quite well also for high degrees of involved polynomials in con- trast to the straightforward way of computing and distributing Unlike their predecessors, the authors try to modify the sound the roots of £· £. over the whole range of frequencies. Such a complex compen- Let us perform a benchmark experiment to compare the exist- sation fully employs the increasing performance of signal hard- ing approach and our newly proposed algorithm for particu- ware dedicated to CD-quality audio signals, and at the same lar numerical data kindly provided by Mikael Sternad and col- time calls for fast and reliable factorization solvers [16]. We leagues from the University of Uppsala. Up to now, two models believe our new algorithm will signiﬁcantly contribute to this of the loudspeakers dynamics have been sent to us for testing goal. purposes and the results related to the more complex one are The loudspeaker dynamics is considered in the form of an ARX presented in the following. The data in concern are given as follows. The numerator proposed routine can be used to factor particular entries of an ´Þ µ ¼ · ½Þ ½ · ¡ ¡ ¡ · ¾ ¼Þ ¾ ¼ is an unstable equivalent diagonal matrix. Unfortunately, from the computa- polynomial of degree 250, ´Þ µ is stable of degree 90, and tional point of view it is not particularly useful as the Smith ½ ¼. Taking ¼, the spectral factorization of Ñ´Þ µ form transformation is known to be numerically fragile [8]. ´Þ µ £ ´Þ µ Ñ¾ ¼ Þ ¾ ¼ · · Ñ¼ · · Ñ¾ ¼ Þ ¾ ¼ is to A modiﬁcation for continuous time polynomials is under re- be performed. In this special case, the spectral factor Ü´Þ µ of search too. However, the suggested approach does not seem to Ñ´Þ µ can be effectively constructed as be as fruitful for continuous time systems as it is in the discrete Ü´Þ µ · ´Þ µ ´Þ µ¡£ Þ time case. The reason is that the relation between the Laplace transform, replacing the role of -transform in the continuous where · are the plus and minus factors of respectively time domain, and the DFT is not so close. The idea of a di- and is the degree of . rect bilinear transform converting the × variable to Þ cannot be All presented experiments were realized on a PC computer with recommended for numerical reasons. Pentium III/1.2GHz processor and 512 MB RAM, under MS Windows 2000 in MATLAB version 6.1. 9 Conclusion Results of this experiment for various values of the parameters A new method for the discrete-time plus-minus factorization Æ are summarized and related in the following table. Namely, problem in the scalar case has been proposed. The new method the computational time and accuracy of results are of interest. relies on numerically stable and efﬁcient FFT algorithm. Be- To obtain the former characteristic, the MATLAB abilities were sides its good numerical properties, the derivation of the routine employed (the built-in functions tic/toc ). The computa- also provides an interesting look into the related mathematics, tional error is deﬁned here as the largest coefﬁcient of the ex- combining the results of the theory of functions of complex pression · , evaluated in the MATLAB workspace, variable, the theory of sampled signals, and the discrete Fourier divided by the largest coefﬁcient of (all in absolute value). transform techniques. The suggested method is employed in a practical application of improving the quality of a hi-ﬁ system. Time [s] Accuracy FFT(14) 0.23 sec ¾ ¡ ½¼ ¿ Acknowledgements FFT(15) 0.45 sec ¾ ¾ ¡ ½¼ FFT(16) 0.89 sec ¡ ½¼ ½½ cı ˇ The work of M. Hromˇ´k and M. Sebek has been supported by FFT(17) 1.75 sec ¾ ¼ ¡ ½¼ ½¾ the Ministry of Education of the Czech Republic under contract No. LN00B096. TABLE 1: Accuracy and efﬁciency of compared algorithms. These tests prove the power of the new algorithm in such tough References examples. Neither of the two procedures described in para- c [1] Kuˇ era V., Analysis and Design of Discrete Linear Con- graph 2 can factor this large polynomial. Direct roots evalua- trol Systems, Academia Prague (1991). tion method, based on the standard MATLAB function roots, gives totally meaningless results (accuracy of ½¼ ¿ ) while the [2] M. Hromcik, J. Jezek, M. Sebek, New Algorithm routine based on spectral factorization fails due to numerical for Spectral Factorization and its Practical Applica- problems with greatest common polynomial divisor evaluation tion, Proceedings of the European Control Conference (Polynomial Toolbox function rdiv was used [8]). ECC’2001, Porto, Portugal, September 1-5, 2001. [3] Z. Hurak, M. Sebek, Algebraic Approach to l-1 Optimal 8 Further Research Control, to appear. [4] Bini D., Pan V., Polynomial and Matrix Computations, The success in modifying a selected numerical procedure, orig- a Volume 1: Fundamental algorithms, Birkh¨ user, Boston inally developed for polynomial spectral factorization, to han- (1994). dle the unsymmetric plus-minus decomposition suggests that other well known spectral factorization routines might work [5] Higham N. J., Accuracy and Stability of Numerical Al- well in the non-symmetric context as well. Actually, we gorithms, S.I.A.M., Philadelphia (1996). achieved some results of this kind recently and plan to incorpo- [6] Needham T., Visual Complex Analysis, Oxford Univer- rate them in the ﬁnal camera-ready version of our contribution. sity Press, 1997. Unfortunately, the described approach cannot be directly ex- c [7] Kuˇ era V., Analysis and Design of Discrete Linear Con- tended to the matrix case - the product of two matrices does not commute in general and since È · ´Þ µÈ ´Þ µ È ´Þ µÈ · ´Þ µ, trol Systems, Academia Prague (1991). one cannot write ÐÒ È ´Þ µ ÐÒ´È · ´Þ µµ · ÐÒ´È ´Þ µµ and per- ˇ [8] Kwakernaak H., Sebek M., PolyX Home Page, form the decomposition. Clearly, in combination with tech- http://www.polyx.cz/, niques for the Smith form of a polynomial matrix [7, 8], the http://www.polyx.com/. [9] Ljung L., System Identiﬁcation: Theory for the User, Prentice-Hall Information and Systems Sciences Series. Englewood Cliffs, Prentice-Hall (1987). [10] Barnett S., Polynomials and Linear Control. Marcel Dekker, New York and Basel (1983). z c [11] Jeˇ ek J. and Kuˇ era V., Efﬁcient Algorithm for Matrix Spectral Factorization,Automatica, vol. 29, pp. 663-669, 1985. cı ˇ [12] Hromˇ´k M., Sebek M., Numerical and Symbolic Com- putation of Polynomial Matrix Determinant, cı ˇ [13] Hromˇ´k M., Sebek M, Fast Fourier Transform and Ro- bustness Analysis with Respect to Parametric Uncertain- ties, Proceedings of the 3rd IFAC Symposium on Robust Control Design ROCOND 2000, Prague, CZ, June 21- 23, 2000. e [14] M. Sternad and A. Ahle´ n, Robust Filtering and Feed- forward Control Based on Probabilistic Descriptions of Model Errors, Automatica, 29, pp. 661-679. e [15] K. Ohrn, A. Ahle´ n and M. Sternad, A Probabilistic Ap- proach to Multivariable Robust Filtering and Open-loop Control, IEEE Transactions on Automatic Control, 40, pp. 405-417. [16] M. Sternad, M. Johansson, J. Rutstrom, Inversion of Loudspeaker Dynamics by Polynomial LQ Feedforward Control, Proceedings of the 3rd IFAC Symposium on Robust Control Design ROCOND 2000, Prague, CZ, June 21-23, 2000. [17] University of Uppsala, Signals and Systems De- partment, Adaptive Signal Processing, Course Homepage 2000, http://www.signal.uu.se/Courses/ CourseDirs/AdaptSignTF/Adapt00.html [18] The Mathworks, Using MATLAB 5.3, The Matrhworks, 1999.