FFT BASED ALGORITHM FOR POLYNOMIAL PLUS-MINUS
Martin Hromˇ´k, Michael Sebek
Centre for Applied Cybernetics
Czech Technical University, Prague, Czech Republic
e-mail: firstname.lastname@example.org, email@example.com
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 .
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  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 . 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
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 , Æ
quadratically optimal ﬁlters for mobile phones [14, 15], and Ð ½ Ô Ô ¼ ½Ô , where Æ
optimal regulators , 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 . 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  and to treat robustness analysis problems of certain
Cauchy’s theorem of argument , the curve Ô´Þ µ for Þ ½
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 . 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 . 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 Ê Ê
¾Ê · ½
Ò´Þ µ ¡¡¡·Ò Þ ·Ò½ ¼ ½ · ¡ ¡ ¡
· Ò ½Þ
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
Ü ´Þ µ Ü · Ü Þ · ¡ ¡ ¡
¼ ½ · Ò½ Þ · ¡¡¡ (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  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.
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 Ê . 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
·¾ Ê ¾Ê ÐÓ Ê · ¾ÐÊ
É½ ´Þ µ ´Þ µ
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 .
 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 .
Figure 2: Optimal ﬁltering problem setup (adopted from ).
As for the spectral factor computation, the authors employ the
Newton-Raphson iterative scheme  in the cited work .
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 ). 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 . 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
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 .
´Þ µ £ ´Þ µ Ñ¾ ¼ Þ ¾ ¼ · · Ñ¼ · · Ñ¾ ¼ Þ ¾ ¼ 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 ¾ ¡ ½¼ ¿
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
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
 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  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 ). ECC’2001, Porto, Portugal, September 1-5, 2001.
 Z. Hurak, M. Sebek, Algebraic Approach to l-1 Optimal
8 Further Research Control, to appear.
 Bini D., Pan V., Polynomial and Matrix Computations,
The success in modifying a selected numerical procedure, orig-
Volume 1: Fundamental algorithms, Birkh¨ user, Boston
inally developed for polynomial spectral factorization, to han-
dle the unsymmetric plus-minus decomposition suggests that
other well known spectral factorization routines might work  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-
 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-
 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- ˇ
 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/.
 Ljung L., System Identiﬁcation: Theory for the User,
Prentice-Hall Information and Systems Sciences Series.
Englewood Cliffs, Prentice-Hall (1987).
 Barnett S., Polynomials and Linear Control. Marcel
Dekker, New York and Basel (1983).
 Jeˇ ek J. and Kuˇ era V., Efﬁcient Algorithm for Matrix
Spectral Factorization,Automatica, vol. 29, pp. 663-669,
 Hromˇ´k M., Sebek M., Numerical and Symbolic Com-
putation of Polynomial Matrix Determinant,
 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-
 M. Sternad and A. Ahle´ n, Robust Filtering and Feed-
forward Control Based on Probabilistic Descriptions of
Model Errors, Automatica, 29, pp. 661-679.
 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,
 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.
 University of Uppsala, Signals and Systems De-
partment, Adaptive Signal Processing, Course
Homepage 2000, http://www.signal.uu.se/Courses/
 The Mathworks, Using MATLAB 5.3, The Matrhworks,