FFT based algorithm for polynomial plus-minus factorization by dffhrtcv3

VIEWS: 10 PAGES: 6

									   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 first order factors or, alternatively, emply a more efficient
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 efficiency 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 filter 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 efficient 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 defined 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, effi-
                                                                      Ý       ݼ ݽ      Ý is given, then its inverse DFT recovers
                                                                      the original vector Ô
cient algebraic design methods for time-optimal controllers [7],                            Æ

quadratically optimal filters 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 fields. 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
identification 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 efficient 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 beneficial compu-                    roots of Ô´Þ µ from infinity 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 (finite) 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
                                                                               infinite power series. Similar reasoning proves the Ô   factor
                   Ô´Þ µ      Ô¼ · Ô½ Þ · ¡ ¡ ¡ · Ô Þ                          desired properties.
nonzero for Þ     ½, we first 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 coefficients Ô
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 Æ
Ô½
              ¼     ½                                    ¼
                  Æ                                                            finite set of its coefficients or by its values in a finite 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.                       ficients), 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 defined as
Ò´Þ µ, obtained from the given Ô´Þ µ, is a Laurent infinite power                              Ê                                Ê
                                                                                                    Ü  
                                                                                                                       ½
series                                                                                                       Ü
                                                                                                                    ¾Ê · ½
                                                                                                Ê                                Ê
          Ò´Þ µ       ¡¡¡·Ò Þ ·Ò½        ¼         ½ · ¡ ¡ ¡
                                             · Ò ½Þ
                                                                               which approximates the Z-transform by dealing with  Ê
It can be directly decomposed,                                                     ·Ê instead of infinite  ½           ·½, 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 coefficients     Ô·
                                                                                                                Ô·
                                                                                                                 ¼        Ô· 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 modification 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 Ê coefficients 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, defined by              be of degree Æ only and only the first Æ · ½ coefficients of
    (2), on the vector                                                  Ô· ´Þ µ should be significant 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 definitions.
     In this way, the set   È     È ¼ Ƚ       È¾Ê of the values
     of Ô´Þ µ at the Fourier points is obtained.                        5 Radix-2 Modification 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 modification 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 defined 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 coefficients 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, defined 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 modification of the
    get ·                   Ê . Proceed with Ü ´Þ µ in obvi-
                   ·        ·
                   ¼                                                    proposed method further increases its efficiency 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 modifications 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 five 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.      filter
                                                                                                             É´Þ µ
Then the particular steps of the modified 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 filter 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
                                                                            ۴ص¹   É´Õ ½ µ
                                                                                    È ´Õ ½ µ
                                                                                                  ٴص   ¹    Õ 
                                                                                                                     ´Õ ½ µ
                                                                                                                     ´Õ ½ µ
                                                                                                                                          Ý ´Øµ   ¹
                                                                                                                                      ·
filter for a moderate quality loudspeaker. By attaching a signal
processor implementing this filter 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 filtering 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-filter 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 significantly 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 modification 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 efficient 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 defined here as the largest coefficient 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 coefficient of (all in absolute value).               transform techniques. The suggested method is employed in a
                                                                            practical application of improving the quality of a hi-fi 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 efficiency 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 final 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 Identification: 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., Efficient 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.

								
To top