System Identification via Polynomial Transformation Method
Prof Pradip Sircar
Department of Electrical Engineering
Indian Institute of Technology Kanpur
Email address: sircar@iitk.ac.in
Abstract
We propose a method based on minimum-variance polynomial approximation to extract system poles
from a data set of samples of the impulse response of a linear system. The method is capable of handling
the problem under general conditions of sampling and noise characteristics. The superiority of the
proposed method is demonstrated by statistical comparison of its performance with the performances of
two exiting methods in the special case of uniform sampling.
Keywords - Minimum-variance polynomial approximation; System identification; Impulse response;
Complex exponential signal
Introduction
In the all-pole model for linear systems, the characterization of the impulse response by a sum of
weighted complex exponential signals, and then estimating the complex frequency- and amplitudeparameters of the signal is equivalent to determining the poles and residues at the poles, respectively, of
the system. This particular problem arises in the diversified fields of control, electromagnetics,
geophysics, NMR spectroscopy, etc., and is, therefore, of eminent importance.
The parameter estimation of complex exponential signals, when the noise corrupted signal samples are
provided at non-uniform spacing, has been attempted in [1], [2]. This problem has special significance
because non-uniform sampling can be preferred over uniform sampling for various considerations like
sampling efficiency, SNR enhancement on the sampled set of signal values, etc.
To deal with the non-uniform sampling case, the idea of applying the orthogonal polynomial
approximation together with the minimum error-variance criterion to find a closed-form expression for
the approximating signal has been introduced in [1], [2]. The closed-form expression can be utilized to
reconstruct the signal values at uniform spacing.
In an extension of that work, it is presented in [3], [4] that reconstruction of signal values at uniform
spacing is through a linear causal transformation based on minimum-variance polynomial approximation.
Therefore, the statistical properties of the approximating error sequence in the reconstructed signal values
can be obtained from the known statistics of the original noise process. We assume time-invariance of the
transformation here.
The maximum likelihood estimator combined with a rank-reduced SVD algorithm can then be employed
to obtain estimation with high accuracy in the general case of non-uniform sampling. The most desirable
feature of the proposed method is that it provides estimation in the maximum likelihood sense even when
the corrupting noise process is not white or Gaussian.

It is the purpose of this article to provide a statistical comparison of the performance of the proposed
method with the performances of two existing SVD-based techniques. Our study shows the proposed
method provides better accuracy of estimation over the existing methods even for the case where the
noise process is white and Gaussian.
Non-Uniform Sampling and Estimation of Parameters of Complex Exponential Signals
The si -parameters of the signal modeled as
M

g (t ) =
∑ exp ( si t ) , t ≥ 0 ; si =−α i + j 2π fi (1)
i =1

g ( tk ) + w ( tk ) ; tk =
t1 , t2 , , t K }
are to be estimated by utilizing the noise corrupted signal samples { x ( tk ) =
at nonuniform spacing, where {w ( tk )} is the sequence sampled from a zero-mean noise process w(t )
whose normalized autocorrelation functions are known.
By employing the orthogonal polynomial approximation and minimum error-variance criterion, the
reconstructed signal vector, YL =
 y ( 0 ) , y (T ) , , y ( ( L − 1)T )  is computed as [4],
T

YL = PQ −1 P1T X K (2)
where X K =
 x ( t1 ) , x ( t2 ) , , x ( t K )  ,
T

1,=
=
P1  p=
 , K ; j 1, , N  ,
j −1 ( ti ) ; i
K
K

Q = diag  ∑ p02 ( tk ), , ∑ pN2 −1 ( tk )  ,
=
 k 1 =k 1


P =  p j −1 ( ( i − 1) T ) ; i = 1, , L ; j = 1, , N  ,
and the polynomials p j (t ) are given by
pj =
( t − a j ) p j −1 − b j p j −2 , j ≥ 1 (3)

with=
p0 1,=
p−1 0 ,
and a=
j

(1 Φ ) ∑ t
K

j −1

k =1

K



Φ j Φ j −1 , Φ j =∑  p j ( tk )  .
k  p j −1 ( t k )  , b j =
2

2

k =1

The approximation order N is determined so that the error variance
K

σ N2 =
∑  x ( tk ) − y ( tk )

2

(K − N )

(4)

k =1

is minimum. The sampling interval T is utilized to convert the z -plane poles into the s -plane poles as,

( 0.5ln z

=
si

2
i

)

+ j arg ( zi ) T (5)

In order to apply the statistics of the approximating error sequence {e=
( iT ) y ( iT ) − g ( iT )} into the
estimation procedure, it is observed that the transformation ( PQ −1 P1T ) in (2) is linear, causal, and its
impulse response sequence {h ( iT )} is obtained in the first column of the transformation matrix. We
assume that the system is time-invariant.
It is then easy to visualize that the noise sequence {w ( tk )} is converted into the error sequence {e ( iT )}
by a linear system {h ( iT )} . As a consequence, since the noise sequence is zero-mean, the error sequence
is also zero-mean; furthermore, the autocorrelation functions are related by
ree =
[ k ] h [ k ] ∗∗ h∗ [ −k ] ∗∗ rww [ k ] (6)

where [ k ] = ( kT ) ,

∗

stands for complex conjugation, and ∗ ∗ denotes discrete convolution.

The probability density function of {e ( iT )} will be jointly Gaussian when the noise process w ( t ) is
Gaussian. In fact, by invoking the central limit theorem, it can be shown that the error sequence will be
close to Gaussian when w ( t ) is white, but may not be Gaussian [4].
By employing the maximum likelihood estimator (MLE) in the linear prediction model well known for
the form of signal expressed in (1), the prediction coefficients {− ai } are obtained by solving the
following equation [4], [5],

(G

H

−1
−1
R ee
G ) A = G H R ee
Y (7)

where Y =  y [ J ] , y [ J + 1] , , y [ L − 1] ,
T

A =−
[ a1 , −a2 ,, −aJ ] ,
T

=
G  g [ J + i − j − 1]=
; i 1, , L − J=
; j 1, , J  ,
R ee=  ree [i − j ]; i= 1, , L − J ; j= 1, , L − J  ,
and, J is the extended model order, J ≤ L − J , and preferably, J is much larger than M .
To compute A from (7), the pseudoinverse of ( G H R ee−1 G ) is needed. The characteristic equation is then
formed as

z J + a1 z J −1 +  + aJ =
0 (8)
which is solved to find M signal and ( J − M ) noise zi -poles.

In practice, the matrix G cannot be formed because of unavailability of { g [i ]} . Hence, (7) is formed by
Y replacing G ,

(Y

H

−1
−1
R ee
Y ) A = Y H R ee
Y (9)

Y  y [ J + i − j − 1]=
; i 1, , L − J =
; j 1, , J  .
where =
In this case, the rank- M pseudoinverse of

(Y

H

−1
R ee
Y ) is employed to compute A , and the SVD

technique is conveniently utilized for the purpose [6]. To get accurate results at very low SNR (say, 5
dB), corrections of the principal singular values are necessary [4].
Uniform Sampling and Other SVD-based Methods
Uniform sampling is the special case of nonuniform sampling. As such, the performance of the developed
method which is based on minimum-variance polynomial approximation and rank-reduced SVD
algorithm can be effectively compared to the performances of other SVD-based methods in uniform
sampling case. We will consider here two such methods which can be employed for parameter estimation

T ) ; k 0, , L − 1} are provided at uniform spacing.
when noise-corrupted signal samples =
{x [ k ] x ( k=

Autocorrelation-like Matrix Method
The autocorrelation-like matrix (ALM) method is an extension of the SVD-Prony method in the secondorder statistics domain [7], [8]. The developed method provides estimation of signal parameters with less
sensitivity to the effect of noise. A brief description of the ALM matrix method is given below.
The linear prediction equations satisfied for the signal given by (1) are expressed in matrix form as,
X A = X (10)

where X and X are defined similar to Y and Y respectively, as in (7) and (9). Then, instead of utilizing
the rank- M pseudoinverse of X to compute A , first both sides of (10) are processed to obtain,
R A = R (11)
L − J − i −1


∗
where
=
R (1 L ) ∑ x [ J + l + i − j ] x=
, I ; j 1, , J  ,
[l ]; i 1,=
l =0


L − J − i −1


and RT (1 L ) ∑ x [ J =
=
+ l + i ] x∗ [l ]; i 1, , I  ,
l =0



and the actual value of I should be considerably less than its maximum value I max = L − J − 1 . The rankM pseudoinverse of R is employed now in determining A , and once the coefficients are known, the
characteristic equation is solved for its roots.

The noise desensitization of the developed method is achieved because first the matrix [ R | R ] can be
separated into the autocorrelation-like signal matrix and autocorrelation noise matrix. Furthermore, the
autocorrelation noise matrix tends to be the null matrix as the data length L approaches infinity. We
assume here that the additive noise is white.
Matrix Pencil Method
The matrix pencil (MP) method has its root in the pencil-of function method which has been extensively
used in system identification [9], [10]. The improved method while extracting the signal pole information
from a generalized singular-value problem, utilizes extended order modeling. A truncated-SVD algorithm
reduces the dimension of the final singular-value problem back to the true model order. As a result, the
noise sensitivity of estimation is substantially reduced, and since only the signal poles are estimated,
separation of signal and noise poles becomes unnecessary.
In

the

MP

method,

two

matrices

are

formed,

X

as

in

(10),

and

X1

given

by

 x [ J + i − j ]=
X
=
; i 1, , L − J ;=
j 1, , J  . The generalized singular-value equation is then written as
1
X1q = z Xq (12)

which is solved for the generalized singular values. In the no-noise case, the M non-zero generalized
singular-values can be shown to be the M zi -poles of the signal.
The generalized singular-value problem can be reduced to a singular-value problem by premultiplying
both sides of (12) with X + which is the rank- M pseudoinverse of X ,
+
X
=
X1 q z=
X + Xq z q (13)

Therefore, the singular-values of the matrix X + X1 are same as the M signal and ( J − M ) noise zi poles. The rank- M pseudoinverse of X is defined as,
=
X+

M

(1 σ ) v u
∑=
i =1

i

i

H
i

V Σ −M1 U H (14)

where {σ i ; i = 1, , M } are the M largest singular values, v i 's and ui 's are the corresponding singular
vectors, V = [ v1 ,, v M ] , U = [u1 , , u M ] , and Σ M = diag [σ 1 , ,σ M ] .
Substituting (14) into (13) and simplifying, it can be shown that the estimates of signal zi -poles can be
found by computing the eigenvalues of the M × M nonsymmetrical matrix, ( Σ −M1 U H X1V ) .

Comparison of Three SVD-based Methods
The comparison is valid only in the uniformly sampling case. Since in this case, the question of
reconstructing the signal values at uniform spacing does not arise, the purpose of processing of data by
minimum-variance polynomial approximation is to be explained.

It is not difficult to show that the autocorrelation-like matrix method will have the asymptotic properties
of optimal estimation, provided the embedded noise is known to be white and Gaussian. For the method
based on minimum-variance polynomial approximation, no such assumptions are needed. The fact that
the preprocessing of data will make the remnant error closely Gaussian is the most desirable feature of the
proposed method, which makes the designed estimator perform like an optimal estimator even in a
realistic situation where the noise process is neither Gaussian nor white. This explains why preprocessing
of data by the polynomial approximation is suggested even in the case where the signal samples are given
at uniform spacing.
Among other advantages of the proposed method are significant enhancement of SNR level in the
processed data, which reduces the effect of ill-conditioning of model equation in estimation [1], [4], and
capability of handling the problem when the original noise process is known to be correlated and colored.
Compared to the first two SVD-based methods, it is observed that the MP method essentially leads to a
suboptimal estimator. The main advantage of the MP method is the reduced dimension in the final
eigenvalue problem, which leads to less computation and as a result, less computational error. The
method possesses high efficiency in implementation because only the signal poles are computed here.
Simulation Results
Example 1: Nonuniform sampling case
The real-valued signal of desired form,
g ( t ) = exp ( s1 t ) + exp ( s1∗ t ) + exp ( s2 t ) + exp ( s2∗ t ) ,

where s1 =
−1 90 + j 2π 0.11 , corrupted with zero-mean white Gaussian noise
−1 75 + j 2π 0.08 and s2 =
is sampled at 50 nonuniformly spaced points which are chosen arbitrarily except that no sampling
interval exceeds 1.1 unit. By applying minimum-variance polynomial approximation, the reconstructed
signal values y ( nT ) are computed at uniform interval of T = 0.5 , as shown in Figure 1. The polynomial
order is chosen to be 19 from the error-variance, σ N2 plot of Figure 2. Observe how the correlated error
process e ( t ) differs from the uncorrelated noise process w ( t ) , both plotted in Figure 3.
By setting the SNR level at 10 , 20 and 40 dB, the magnitude of bias and the variance of estimates of
each parameter are numerically computed from 100 independent realizations of noise sequences. The
results are shown in Table 1. The extended model order is chosen to be 16 for optimal performance in all
the cases.

Figure 1 Orthogonal polynomial approximation: Peak SNR = 5 dB

Figure 2 Minimum error-variance criterion: Approximation order = 19

Figure 3 Preprocessing by minimum-variance polynomial approximation

Table 1 Bias and Variance of Parameter Estimates
Estimated
Parameters
α1
f1
α2
f2

SNR = 40 dB
Bias
Variance
2.576 × 10−4
1.174 × 10−7
7.659 × 10−6
2.675 × 10−8
4.118 × 10−5
2.859 × 10−7
5.677 × 10−5
5.772 × 10−9

SNR = 20 dB
Bias
Variance
7.880 × 10−4
8.303 × 10−6
8.240 × 10−5
2.360 × 10−7
6.330 × 10−5
2.235 × 10−6
1.492 × 10−4
2.649 × 10−7

SNR = 10 dB
Bias
Variance
1.418 × 10−3
1.409 × 10−5
1.431 × 10−4
2.353 × 10−6
5.521 × 10−4
7.067 × 10−6
2.621 × 10−4
1.027 × 10−6

Example 2: Uniform sampling case
In this example, the transient signal
g ( t ) = β1 exp ( s1 t ) + β1 exp ( s1∗ t ) + β 2 exp ( s2 t ) + β 2 exp ( s2∗ t ) ,

with s1 =
−0.00555 + j 0.08 , β1 = 1.5 , s2 =
−0.00666 + j 0.11 , β 2 = 3.5 , is sampled at 50 equispaced
points with interval 5.6 units. The SNR value is set at 5 dB by mixing zero-mean white Gaussian noise
sequence generated by computer.
The magnitude of bias and the variance of estimation for each parameter are computed from 200
independent realizations, by applying in turn the proposed method, the ALM method and the MP method.
The comparative results, as shown in Table 2, agree well with the discussions presented in Section 3
where it is argued that the proposed method will be superior in accuracy of estimation.

Table 2 Comparison of Three SVD-based Methods
Estimated
Parameters
α1
f1
α2
f2

Our Method
Bias
Variance
1.006 × 10−4
9.298 × 10−5
8.469 × 10−4
7.983 × 10−5
3.993 × 10−4
1.820 × 10−5
3.212 × 10−4
1.298 × 10−5

ALM Method
Bias
Variance
4.704 × 10−3
7.333 × 10−4
3.582 × 10−3
1.301 × 10−4
6.589 × 10−3
2.518 × 10−5
4.114 × 10−3
2.507 × 10−5

MP Method
Bias
Variance
6.583 × 10−3
9.973 × 10−3
8.790 × 10−3
8.049 × 10−3
6.934 × 10−3
3.346 × 10−4
7.063 × 10−3
2.889 × 10−4

For the MP method, the extended model order is set at 16 to provide best results. For the other two
methods, the order is optimally chosen to be 20.
Concluding Remarks
We have proposed a complete approach based on minimum-variance polynomial approximation,
maximum likelihood estimation and SVD algorithm. The proposed method provides very accurate
estimates of parameters of complex exponentials under general conditions, viz., when the signal samples
are not necessarily at uniform spacing, and/or the superimposed noise process may not be white or
Gaussian.
The superiority of the proposed method is demonstrated by comparing its performance with the
performance of the SVD-based autocorrelation-like matrix method and matrix pencil method in the
uniform sampling case. Better results are obtained by the proposed method even when the noise process is
white and Gaussian. Therefore, we may conclude that the proposed method with pre-processing of data
should be the choice for better accuracy in all cases, although it needs some extra computation.
References
[1] Sircar, P.: Accurate Parameter Estimation of Damped Sinusoidal Signals Sampled at Nonuniform
Spacings. PhD Dissertation, Elec. Eng., Syracuse University (1987)
[2] Sircar, P., Sarkar, T.K.: System Identification from Nonuniformly Spaced Signal Measurements.
Signal Processing. 14(3), 253–268 (1988)
[3] Ranade, A.C.: System Identification from Nonuniformly Sampled Data. MTech Thesis, Elec. Eng.
Dept., IIT Kanpur (1990)
[4] Sircar, P., Ranade, A.C.: Nonuniform Sampling and Study of Transient Systems Response. IEE Proc.
F. 139(1), 49–55 (1992)
[5] Kay, S.M.: Modern Spectral Estimation: Theory and Application. Prentice Hall, Englewood Cliffs,
N.J. (1988)
[6] Deprettere, E.F. (ed.): SVD and Signal Processing: Algorithms, Applications and Architectures.
Elsevier Science, North-Holland (1988)

[7] Kumaresan, R., Tufts, D.W.: Estimating the Parameters of Exponentially Damped Sinusoids and PoleZero Modeling in Noise. IEEE Trans. on Acoust., Speech, Signal Processing. 30(6), 833–840 (1982)
[8] Cadzow, J.A., Wu, M.M.: Analysis of Transient Data in Noise. IEE Proc. F. 134(1), 69–78 (1987)
[9] Jain, V.K., Sarkar, T.K., Weiner, D.D.: Rational Modeling by Pencil-of-Functions Methods. IEEE
Trans. on Acoust., Speech, Signal Processing. 31(3), 564–573 (1983)
[10] Hua, Y., Sarkar, T.K.: Matrix Pencil Method for Estimating Parameters of Exponentially
Damped/Undamped Sinusoids in Noise. IEEE Trans. on Acoust., Speech, Signal Processing. 38(5), 814–
824 (1990)
Further reading
•

Ravi Shankar, M., Sircar, P.: Nonuniform sampling and polynomial transformation method. In
Proc. IEEE International Conference on Communications, ICC 2002, 3, 1721–1725 (2002)

•

Banoth, J.K., Sircar, P.: Polynomial transformation method for non-Gaussian noise
environment. arXiv preprint arXiv:1401.5580 (2014)

•

Sircar, P.: Mathematical Aspects of Signal Processing. Cambridge University Press (2016)

