Analysis of fast-sampled non-linear systems: Generalised frequency response functions for δ -operator models

ARTICLE IN PRESS

Signal Processing 86 (2006) 3246–3257 www.elsevier.com/locate/sigpro

Analysis of fast-sampled non-linear systems: Generalised frequency response functions for d-operator models M.A. Chadwick, V. Kadirkamanathan, S.A. Billings Department of Automatic Control and Systems Engineering, The University of Sheffield, Mappin Street, Sheffield S1 3JD, UK Received 13 April 2005; received in revised form 4 November 2005; accepted 15 January 2006 Available online 9 February 2006

Abstract It is widely known that non-linear systems can respond at frequencies other than the excitation frequency. This harmonic distortion must be accommodated when sampling system data for use in identification and can result in high sample frequencies to ensure information about the underlying continuous-time system is not lost. Unfortunately, high sampling rates can have adverse effects on algorithms when used with the absolute-positional operator q. However, ddomain models do not suffer from such problems and are more suited to fast-sampled systems. Therefore, it would appear that the d-domain provides a more robust and natural basis for identification and analysis of non-linear systems when fastsampled data are involved. In this paper a new framework is proposed for analysing polynomial non-linear d-domain models in the higher-order transform and frequency domains. The use of the NARMAX approach is advocated and is adapted to cater for non-linear models of the d-domain description. r 2006 Elsevier B.V. All rights reserved. Keywords: Non-linear system identification; NARMAX model; d-operator; Fast-sampled systems; Generalised frequency response functions (GFRFs)

1. Introduction In linear system analysis a major tool in understanding the dynamics of a system is the study of the frequency response of the system transfer function, to reveal the changes to the gain and phase of the inputs as they pass independently through the system. This implies the range of frequencies over which the output has power is governed by the range of frequencies in the input. This underlying Corresponding author. Tel.: +44 (0)114 222 5553;

fax: +44 (0)114 222 5661. E-mail addresses: m.chadwick@sheffield.ac.uk, [email protected] (M.A. Chadwick).

principle led to the guidelines given by Nyquist [1] and formalised later by Shannon [2] that a signal should be sampled at least twice the cut-off frequency (or highest frequency of interest). Non-linear systems, however, can generate power both at frequencies over the input power range and at new frequency locations to produce outputs which are rich in spectral content. This phenomena is well known and arises from the effects of harmonics and inter-modulation generation [3]. Sampling of non-linear systems must therefore take account of this harmonic distortion, for otherwise information about the non-linear dynamics at higher frequencies may be lost. In the non-linear case therefore, the sampling frequency should be

0165-1684/$ - see front matter r 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.sigpro.2006.01.004

ARTICLE IN PRESS M.A. Chadwick et al. / Signal Processing 86 (2006) 3246–3257

linked to the frequency content of the output taking into account this distortion [4]. The problematic case occurs when the non-linear system generates new power at higher frequencies which in turn requires fast sampling. Fast sampling induces several problems when analysing systems, for example, (i) analysing vast amounts of data collected, (ii) large wordlengths and numerical ill-conditioning within algorithms [5], primarily caused by the discrete time-shift operator q, and (iii) the encroachment of the sample frequency onto the noise bandwidth causing aliasing/frequency folding. These drawbacks have led to a renewed interest in the gradient-based discrete-time operator d, used to parameterise models of fast-sampled discrete-time systems [6]. These methods have two main benefits for fast-sampled data, namely (i) less numerical illconditioning compared to the absolute-positional shift operator, q [7] and (ii) an inherent ability to unify continuous- and discrete-time systems [8]. Currently the application of d-domain techniques includes: identification and modelling [9–11], predictive control [12,13], H1 control [14] and adaptive control in [15]. In [6] the authors provide a thorough treatment of frequency response analysis and design for linear d-domain models. They discuss the issue of qdomain models falling foul of fast-sampled systems and the lack of connection between the q-domain and continuous-time poles of the equivalent discrete- and continuous-time systems. This lack of connection is further reinforced by the clustering of the dynamic information to the point (0,1) in the zplane when the sampling interval becomes very small [16]. Whereas for d-domain models, the locations of the poles and zeros are preserved, and as the sample time tends to zero, the d-domain poles and zeros converge to the continuous-time pole/zero locations. One method of investigating non-linear phenomena involves analysing the generalised frequency response functions (GFRFs) of the non-linear system. The GFRFs were, up until a few decades ago, found by taking multi-dimensional Fourier transforms of the Volterra kernels. When applied in practise, this involves taking multi-dimensional fast Fourier transforms of data collected from the nonlinear system [17]. These methods have severe limitations for application to real systems as multidimensional windows are needed and large amounts of data may need to be processed to obtain smooth

3247

estimates, making these methods computationally burdensome [18]. An alternative parametric method was developed by Billings et al. that avoids most of these restrictions and is based on the NARMAX approach [18–20]. This method involves estimating a NARMAX (Non-linear AutoRegressive Moving Average model with eXogenous inputs) model and computing the GFRFs directly from the model. Work done thus far on frequency response methods for d-operator parameterised models is limited to linear systems. There is currently no framework from which to analyse polynomial nonlinear d-domain models in the higher-order transform and frequency domains. This paper aims to address these issues and to introduce the reader to the d-domain discrete-time framework for fastsampled non-linear systems. The paper is structured as follows: Section 2 gives a brief review of the d-operator, Section 3 summarises continuous Volterra series theory with multi-dimensional Laplace and Fourier transforms. Section 4 focuses on the discrete Volterra series combined with d-operators, multi-dimensional Delta transforms and how to map such models into the higher-order frequency domain. Section 5 introduces the d-NARMAX model and describes the procedure by which the GFRFs may be extracted from models using the probing method. Section 6 presents the results of a d-NARMAX Duffing equation example, and finally conclusions are drawn in Section 7. 2. A review of the d-operator A linear discrete time-invariant system can be represented as an Auto Regressive Moving Average model with eXogenous inputs (ARMAX). In the forward-shifted sense, the model takes the form yðk þ nÞ ¼
an
1 yðk þ n
1Þ
  
a0 yðkÞ þ bn
1 uðk þ n
1Þ þ    þ b1 uðk þ 1Þ þ b0 uðkÞ þ eðk þ nÞ þ cn
1 eðk þ n
1Þ þ    þ c1 eðk
1Þ þ c0 eðkÞ,

ð1Þ

where n is the order of the system and uðkÞ, yðkÞ and eðkÞ are taken as the values of the input, output and prediction error sequences, respectively, at time kT. The model can be seen to be constructed using weighted (scaled) sums of past inputs, outputs and prediction error sequences through the model parameters ai , bi and ci . This description is often shortened by the use of an absolute-positional shift

ARTICLE IN PRESS M.A. Chadwick et al. / Signal Processing 86 (2006) 3246–3257

3248

operator, q where

3. Continuous non-linear system representations

qyðkÞ ¼ yðk þ 1Þ.

(2)

Using the q operator, Eq. (1) can be shortened to yðk þ nÞ ¼


n
1 X

ai qi yðkÞ þ

i¼0

þ

n
1 X

n
1 X

bi qi uðkÞ

i¼0

ci qi eðkÞ þ eðk þ nÞ.

ð3Þ

i¼0

One major drawback in representing models of discrete-time systems in the above format is the lack of connection to the continuous-time model equivalent. This is due to the inherent differences between the operators q and d=dt used to parameterise each model. An alternative approach can be adopted by using the d-operator defined as a first-forward difference (4)

where q is the forward shift operator and T is the sampling interval. The d-operator, due to its definition, has a convenient property, given a differentiable signal, yðkÞ

T!0

d yðkÞ. dt

(5)

Given the ARMAX model of Eq. (3) in the qdomain, a d-domain equivalent may be constructed by substituting qn ¼ ð1 þ TdÞn

(6)

to give dn yðkÞ ¼


n
1 X i¼0

þ

yðtÞ ¼

1 X

yn ðtÞ,

(8)

n¼1

where each of the higher-order outputs are given by Z 1 Z 1  hn ðt1 ; . . . ; tn Þ yn ðtÞ ¼
1 n Y




1

uðt
ti Þ dti ;

n40.

ð9Þ

i¼1

q
1 d¼ , T

lim dyðkÞ ¼

The Volterra series is the classical method of representing a wide class of non-linear systems [22]. The series is composed of a zero-memory (static) non-linearity (represented by a Taylor series) and the convolution integral, and comprises of a possibly infinite number of parallel sub-systems called Volterra functionals which are each described by the Volterra kernels hn ðt1 ; . . . ; tn Þ. The Volterra series can be expressed as

n
1 X

a¯ i di yðkÞ þ

n
1 X

b¯ i di uðkÞ

i¼0

c¯ i di eðkÞ þ dn eðkÞ.

ð7Þ

i¼0

Notice that the model order is preserved and the only change in model representation is the replacement of d for q and the adjustment in the parameters. This indicates that any polynomial in q of degree n will be exactly equivalent to some polynomial in d of degree n [21]. This representation then enables a more intuitive comparison between the discrete-time system and the underlying continuous-time system.

Each of the sub-systems represented by Eq. (9) is called degree-n homogeneous systems. In situations where the non-linearities are not too harsh or when a truncation is sufficient, a finite number of subsystems will suffice in representing the non-linear system. 3.1. Multi-dimensional Laplace transforms Laplace transforms have been pivotal in understanding the dynamics of linear systems and are a powerful and integral part of linear control theory. Since engineers are familiar with Laplace transforms of linear systems, it would be advantageous to have a similar construct for non-linear systems. Unfortunately, due to the dimensional properties, there is no direct way of doing this. An indirect way is adopted by amending Eq. (9) and introducing an auxiliary multi-dimensional time function [23] Z 1 Z 1 yn ðt1 ; . . . ; tn Þ ¼  hn ðt1 ; . . . ; tn Þ
1 n Y




1

uðti
ti Þ dti ,

ð10Þ

i¼1

where the real output yn ðtÞ is recovered by the restriction yn ðtÞ ¼ yn ðt1 ; . . . ; tn Þjt1 ¼¼tn ¼t .

(11)

Taking a multi-dimensional Laplace transform of Eq. (10) and assuming causality i.e., hn ðt1 ; . . . ; tn Þ ¼ 0 8ti o0 i ¼ 1; 2; . . . ; n yields

ARTICLE IN PRESS M.A. Chadwick et al. / Signal Processing 86 (2006) 3246–3257

Y n ðs1 ; . . . ; sn Þ ¼ H n ðs1 ; . . . ; sn ÞUðs1 Þ . . . Uðsn Þ

(12)

which, for the case when n ¼ 1, reduces to the familiar linear transfer function definition. This result implies that hn ðt1 ; . . . ; tn Þ and H n ðs1 ; . . . ; sn Þ provide two equivalent representations in the timeand transform-domain, respectively, for the nthorder homogenous sub-system [3].

d-domain procedure and is different from the traditional q-domain discretisation process used throughout the existing literature. 4.1. Multi-dimensional Delta ðDÞ transforms Consider the linear (or one-dimensional) Dtransform definition as follows [6]:

3.2. Multi-dimensional Fourier transforms F ðgÞ ¼ Dff ðkÞg ¼ The transfer function of the system was introduced using the Laplace transform. This theory is often complemented by the frequency-domain counterpart namely, the Fourier transform. For a homogenous system of nth degree taking the multidimensional Fourier transform of the kernel yields Z 1 Z 1 H n ðjo1 ; . . . ; jon Þ ¼  hn ðt1 ; . . . ; tn Þ
1
1
jðo1 t1 þþon tn Þ

e

4. Discrete-delta non-linear system representations

f ðkÞ½1 þ gT
k T.

(17)

A multi-dimensional form of the D-transform can therefore be defined as a direct extension of Eq. (17) to higher dimensions as follows: F n ðg1 ; . . . ; gn Þ ¼ Dn ff n ðk1 ; . . . ; kn Þg 1 1 X X ¼  f n ðk1 ; . . . ; kn Þ k1 ¼0 n Y



kn ¼0

½1 þ gm T
km T,

ð18Þ

m¼1

where g1 ; . . . ; gn are complex numbers and are assumed to be within the region of convergence to ensure the existence of F n ðg1 ; . . . ; gn Þ. Amending Eq. (16) and introducing a multidimensional sample time function

H n ðjo1 ; . . . ; jon Þ ¼ H n ðs1 ; . . . ; sn Þjs1 ¼jo1 ;...;sn ¼jon . (14)

1 X k¼0

dt1 . . . dtn ð13Þ

which are called the GFRFs of the system transfer function [24]. There is an obvious link between the transforms of Laplace and Fourier: if a Laplace transfer function H n ðs1 ; . . . ; sn Þ exists for Re½si  ¼ 0, i ¼ 1; . . . ; n, the Fourier transfer function is given by the relationship

3249

yn ðk1 ; . . . ; kn Þ ¼

1 X



i1 ¼0 n Y



1 X

hdn ði1 ; . . . ; in Þ

in ¼0

uðkm
im ÞT,

ð19Þ

m¼1

As shown previously, a non-linear system may be represented by the Volterra series in the continuoustime domain. The direct parallel to this series in the discrete time domain is given by the following: yðkÞ ¼

1 X

yn ðkÞ,

(15)

n¼1

where each of the higher-order outputs are given by yn ðkÞ ¼

1 X



i1 ¼0 n Y



1 X in ¼0

k ¼ 0; 1; 2 . . .

yn ðkÞ ¼ yn ðk1 ; . . . ; kn Þjk1 ¼¼kn ¼k

(20)

then the multi-dimensional D-transform can be applied directly to Eq. (19). Assuming causality i.e., hdn ðk1 ; . . . ; kn Þ ¼ 0; 8ki o0, i ¼ 1; . . . ; n gives Y n ðg1 ; . . . ; gn Þ ¼ H n ðg1 ; . . . ; gn ÞUðg1 Þ . . . Uðgn Þ,

(21)

where for the case when n ¼ 1, Eq. (21) reduces to the familiar linear input–output relationship Y ðgÞ ¼ HðgÞUðgÞ. For a full derivation of the above result see the Appendix.

hdn ði1 ; . . . ; in Þ

uðk
im ÞT;

where the real output is recovered by the restriction

ð16Þ

m¼1

and the kernel hdn ði1 ; . . . ; in Þ is a real sequence which is equal to zero if any argument is negative (onesided). The superscript d notation is used to indicate the kernel above is derived by discretising the continuous Volterra series in accordance with the

4.2. Mapping multi-dimensional D-transforms into the higher-order frequency domain Multi-dimensional D-transforms may be mapped directly into the higher-order frequency domain by considering the function H n ðg1 ; . . . ; gn Þ evaluated on

ARTICLE IN PRESS M.A. Chadwick et al. / Signal Processing 86 (2006) 3246–3257

3250

a hypersphere jg1 j ¼ jg2 j ¼    ¼ jgn j ¼ 1=T centred at
1=T with radius 1=T where T is the sampling interval. The Fourier transfer function is then given by the mathematical relationship H dn ðjo1 ; . . . ; jon Þ ¼ H n ðg1 ; . . . ; gn Þjg1 ¼ðejo1 T
1Þ=T;...;gn ¼ðejon T
1Þ=T .

ð22Þ

It can be seen that as the sampling interval T ! 0, the hypersphere expands to cover the region defined by Re½gi o0, i ¼ 1; . . . ; n of the multi-dimensional complex domain with the boundary becoming the imaginary axis. This implies for small sampling intervals, the multi-dimensional g-domain is very similar to the multi-dimensional s-domain. This leads to the higher-order frequency domain of ddomain models being approximately equal to continuous-domain models for small values of the sampling interval; and in the limit, exactly equal. 4.3. Properties of the d-operator and the multidimensional D-transform One inherent characteristic of the d-operator is the key property that it converges to the continuous-time derivative in the fast-sampling limit. A consequence of this property is that the onedimensional D-transform converges to the associated Laplace transform when T ! 0. This convergence property also holds for the multidimensional D-transform where

5. Non-linear d-domain analysis 5.1. The d-NARMAX model The NARMAX (Non-linear AutoRegressive Moving Average with eXogenous inputs) model has received much attention in recent decades due to the ability to represent a wide class of non-linear systems using a choice of suitable basis functions [25,26]. Since a NARMAX model includes regressors in both the inputs and outputs, it has a much smaller parameter set than functional series expansions such as the Volterra series. The NARMAX model can therefore be visualised as a Volterra series in the inputs and outputs [27]. Traditionally, the NARMAX model is written in terms of the backward-shift operator q
1 , but will be written here in terms of the forward-shift operator q; this is to provide a more amenable relationship between the forward-shifted d-operator used throughout this paper. The NARMAX model can be represented as yðk þ nÞ ¼ a þ F l ½yðk þ n
1Þ; . . . ; yðkÞ, uðk þ n
1Þ; . . . ; uðkÞ, eðk þ n
1Þ; . . . ; eðkÞ þ eðk þ nÞ,

lim Y n ðg1 ; . . . ; gn Þ ¼ Y n ðs1 ; . . . ; sn Þjs1 ¼gn ;...;sn ¼gn .

T!0

(23) Although this property is widely accepted, it is also understood that the fast-sampling limit can never practically be reached. Other interesting features of the D-transforms are the relation to the Z-transform commonly used in the literature. The discrete multi-dimensional D-transform can be related to the multi-dimensional Z-transform by noting that Y n ðg1 ; . . . ; gn Þ ¼ T n Y n ðz1 ; . . . ; zn Þjz1 ¼1þg1 T;...;zn ¼1þgn T , (24) n

where Y n ðz1 ; . . . ; zn Þ ¼ Z ½yn ðk1 ; . . . ; kn Þ. Conversely, Y n ðz1 ; . . . ; zn Þ ¼ T
n Y n ðg1 ; . . . ; gn Þjg1 ¼ðz1
1Þ=T;...;gn ¼ðzn
1Þ=T .

mains. This reinforces the ability to unify both domains and motivates the investigation of these operators in this study. For a more thorough investigation into the origins and properties of these operators, the reader is referred to [6].

ð25Þ

These properties show the direct connection the doperator and multi-dimensional D-transform have with both the continuous and discrete time do-

ð26Þ

where uðkÞ, yðkÞ and eðkÞ are the sampled input, output and prediction error sequences, respectively, F is some non-linear function (taken as polynomial in this case), l is the degree of non-linearity, n is the order of the system and a is a constant (dc-shift) term. A polynomial non-linear d-domain model, termed the d-NARMAX model, may be represented by the following equation: l dn yðkÞ ¼ a þ F¯ ½dn
1 yðkÞ; . . . ; d0 yðkÞ,

dn
1 uðkÞ; . . . ; d0 uðkÞ, dn
1 eðkÞ; . . . ; d0 eðkÞ þ dn eðkÞ,

ð27Þ

where Eq. (27) can be seen to be a direct mapping of the polynomial q-domain NARMAX model of Eq. (26) into the d-domain. The structure of Eq. (27) is the same as that of Eq. (26), except instead of using forward-shifts to parameterise the model, forwarddifferences are used. If the data sampled from the

ARTICLE IN PRESS M.A. Chadwick et al. / Signal Processing 86 (2006) 3246–3257

system are contaminated by wideband noise, taking differences will greatly amplify the noise and should be avoided. Measures have been taken to cope with this issue and are discussed in [6,28,7]. Although the d-domain model terms should have a direct correspondence to the derivative terms that naturally describe the continuous-time physical system, the parameter estimates will necessarily be different due to the discrete nature of the model. Therefore, the d-domain provides an exact discretisation framework that incorporates the key attribute of directly linking model terms to the continuous-time domain. It is important to make the distinction between the use of the d-operator as a discrete-time parameterisation and a continuoustime approximation; only when the sampling interval is very small compared to the dynamics of the system can a d-domain model be said to approximate the underlying continuous system.

5.2. The probing method for the discrete Volterra series The probing (or harmonic input) method [29] investigates the steady-state output of a non-linear system subjected to a sum of sinusoids. Consider a P-tone input signal consisting of an ensemble of sinusoids with different amplitude, frequency and phase represented by uðkÞ ¼

P X

jAp j cosðop kT þ ffAp Þ

p¼1

¼

P 1X ðAp ejop kT þ Ap e
jop kT Þ 2 p¼1

P 1 X Ap ejop kT , ¼ 2 p¼
P

ð28Þ

yn ðkÞ ¼



i1 ¼0



1 X in ¼0

n P Y 1 X m¼1



2 p¼
P

Ap ejop ðk
im ÞT T

1 X

n

¼

n Y

Apm ejopm ðk
im ÞT T

m¼1

P P n Y X 1 X    Apm ejopm kT 2n p ¼
P p ¼
P m¼1 n

1



1 X i1 ¼0

¼

1

hdn ði1 ; . . . ; in Þ

in ¼0



1 X

hdn ði1 ; . . . ; in Þ

in ¼0

n Y

e
jopm im T T

m¼1

P P X 1 X    Ap1 . . . Apn H dn 2n p ¼
P p ¼
P n

1

ðjop1 ; . . . ; jopn Þejðop1 þþopn ÞkT

ð29Þ

which coupled with Eq. (15) yields the overall nonlinear system output in the frequency domain. This equation can simply be viewed as the direct parallel to that of the continuous domain in [3] but written in the delta domain used throughout this paper. Note also that since the output of the system is real, the right-hand side of the equation contains conjugate pairs. The above analysis was carried out assuming the input was applied to a real system so it was constructed as an ensemble of K ‘real’ sinusoids. In an analytical setting where a model of the nonlinear system is available, a ‘convenience of the theoretical world’ may be drawn on where an ideal (or non-physical) sinusoidal input can be applied to the mathematical model, and the GFRFs obtained more readily. The form of these inputs are a restriction of Eq. (28) where the ideal sinusoid is taken as the conjugate with positive exponent and is represented as follows: P X

Ap ejop kT .

(30)

p¼1

Applying Eq. (30) to Eq. (16) and setting the number of sinusoids equal to the degree of nonlinearity i.e., P ¼ n and the magnitude of the input signals to unity i.e., Ap ¼ 1 for all p ¼ 1; 2; . . . ; n yields yn ðkÞ ¼

hdn ði1 ; . . . ; in Þ

P P X 1 X 1 X   n 2 p ¼
P p ¼
P i ¼0 1

uðkÞ ¼

where Ap is a complex number to give the amplitude and phase of the pth frequency in the input with A0 ¼ 0, A
p ¼ Ap and o
p ¼
op . Substitution of Eq. (28) into the higher-order Volterra functional model of Eq. (16) yields 1 X

¼

3251

n X p1 ¼1



n X

H dn ðjop1 ; . . . ; jopn Þejðop1 þþopn ÞkT

pn ¼1

(31) which describes the output response of the system to the ideal sinusoidal input. This equation is a subset of Eq. (29) and represents the conjugates with

ARTICLE IN PRESS 3252

M.A. Chadwick et al. / Signal Processing 86 (2006) 3246–3257

positive exponents contained in Eq. (29). Only half of the conjugate pairs are needed to characterise the function H dn ðjop1 ; . . . ; jopn Þ, as the conjugates with negative exponents are ‘mirror images’ of their positive counterparts. Eq. (31) contains all the necessary information to compute the function H dn ðjop1 ; . . . ; jopn Þ as in Eq. (29) yet, by exploiting the reflectional symmetry in the o1 ¼ o2 axis, requires less computation. This then enables the GFRFs to be extracted analytically by equating the coefficients of n!ejðo1 þþon ÞkT in the system output when the input is defined as in Eq. (30) with P ¼ n and Ap ¼ 1 for all p ¼ 1; 2; . . . ; n. 6. A d-NARMAX Duffing equation simulation example

md ¯ yðkÞ þ c¯ dyðkÞ þ k¯ 1 yðkÞ þ k¯ 2 y ðkÞ þ k¯ 3 y3 ðkÞ ¼ uðkÞ,

uðkÞ ¼ ejokT .

(34)

From Eqs. (15) and (31) with n ¼ 1 yðkÞ ¼ H 1 ðJðoÞÞejokT ,

(35)

where the forward-differences of inputs and outputs are computed using the definition q
1 T therefore

(36)

d¼

uðk þ 1Þ
uðkÞ T joT ¼ ½e
1T
1 ejokT ¼ JðoÞejokT

duðkÞ ¼

Mapping a d-NARMAX model into the GFRFs enables an analysis of the transfer function frequency response to be carried out revealing information about the non-linear systems dynamics. Since the frequency domain is the invariant descriptor of a system, irrespective of the model representation in the time domain, analysing these properties gives valuable insight into the system characteristics. For such an understanding to be gained, the probing method has to be adapted to accommodate models of the d-NARMAX description. The probing or harmonic input method introduced in [29], and later developed into a recursive algorithm for application to NARMAX models in [30], can readily be applied to d-NARMAX models. Consider the d-domain Duffing equation below, which is commonly used to represent mechanical oscillators such as beams or wires/strings under tension. The system contains a constant damping factor, c¯ and a non-linear displacement/stiffness term k¯ 1 yðtÞ þ k¯ 2 y2 ðtÞ þ k¯ 3 y3 ðtÞ representing a nonlinear mechanical spring, and can be expressed as 2

This will be used for mathematical convenience which will become more apparent shortly and will frequently be written as JðoÞ for simplicity. Defining the first probing input as

2

ð37Þ

and dyðkÞ ¼ JðoÞH 1 ðJðoÞÞejokT .

(38)

Substituting into (32) yields mJ ¯ 2 ðoÞH 1 ðJðoÞÞejokT þ c¯ JðoÞH 1 ðJðoÞÞejokT þ k¯ 1 H 1 ðJðoÞÞejokT þ k¯ 2 ½H 1 ðJðoÞÞejokT 2 þ k¯ 3 ½H 1 ðJðoÞÞejokT 3 ¼ ejokT

ð39Þ

and equating coefficients of ejokT on both sides yields H 1 ðJðoÞÞ ¼

1 . mJ ¯ ðoÞ þ c¯ JðoÞ þ k¯ 1 2

(40)

Probing with two exponentials uðkÞ ¼ ejo1 kT þ ejo2 kT .

(41)

From Eq. (31) the output is yðkÞ ¼ H 1 ðJðo1 ÞÞejo1 kT þ H 1 ðJðo2 ÞÞejo1 kT þ 2!H 2 ðJðo1 Þ; Jðo2 ÞÞejðo1 þo2 ÞkT þ H 2 ðJðo1 Þ; Jðo1 ÞÞejð2o1 ÞkT

ð32Þ

where m, ¯ c¯ and k¯ i , i ¼ 1; 2; 3 are constants Before the probing method for d-operator models is introduced, consider the operator proposed in [8] and defined as ( ðejoT
1Þ=T : T40; Jðo; TÞ ¼ (33) jo : T ¼ 0:

þ H 2 ðJðo2 Þ; Jðo2 ÞÞejð2o2 ÞkT

ð42Þ

and the forward difference of the output is dyðkÞ ¼ JðoÞH 1 ðJðo1 ÞÞejo1 kT þ JðoÞH 1 ðJðo2 ÞÞejo1 kT þ 2!Jðo1 þ o2 ÞH 2 ðJðo1 Þ; Jðo2 ÞÞejðo1 þo2 ÞkT þ Jð2o1 ÞH 2 ðJðo1 Þ; Jðo1 ÞÞejð2o1 ÞkT þ Jð2o2 ÞH 2 ðJðo2 Þ; Jðo2 ÞÞejð2o2 ÞkT .

ð43Þ

ARTICLE IN PRESS M.A. Chadwick et al. / Signal Processing 86 (2006) 3246–3257

Substituting into (32) and equating coefficients of 2!ejðo1 þo2 ÞkT yields mJ ¯ 2 ðo1 þ o2 ÞH 2 ðJðo1 Þ; Jðo2 ÞÞ þ c¯ Jðo1 þ o2 ÞH 2 ðJðo1 ÞÞ

þ k¯ 1 H 2 ðJðo1 ÞÞ þ k¯ 2 H 1 ðJðo1 ÞÞH 1 ðJðo2 ÞÞ ¼ 0 ð44Þ therefore
k¯ 2 H 1 ðJðo1 ÞÞH 1 ðJðo2 ÞÞ mJ ¯ ðo1 þ o2 Þ þ c¯ Jðo1 þ o2 Þ þ k¯ 1 2

(45) and since H 1 ðJðo1 þ o2 ÞÞ ¼

1 mJ ¯ ðo1 þ o2 Þ þ c¯ Jðo1 þ o2 Þ þ k¯ 1 (46) 2

then H 2 ðJðo1 Þ; Jðo2 ÞÞ ¼
k¯ 2 H 1 ðJðo1 ÞÞH 1 ðJðo2 ÞÞH 1 ðJðo1 þ o2 ÞÞ. ð47Þ Continuing the process further probing with three exponentials and equating coefficients of 3!ejðo1 þo2 þo3 ÞkT gives H 3 ðJðo1 Þ; Jðo2 Þ; Jðo3 ÞÞ ¼
16H 1 ðJðo1 þ o2 þ o3 ÞÞ f4k¯ 2 ðH 1 ðJðo1 ÞÞÞH 2 ðJðo2 Þ; Jðo3 ÞÞ þ H 1 ðJðo2 ÞÞH 2 ðJðo3 Þ; Jðo1 ÞÞ þ H 1 ðJðo3 ÞÞH 2 ðJðo1 Þ; Jðo2 ÞÞ þ 6k¯ 3 H 1 ðJðo1 ÞÞH 1 ðJðo2 ÞÞH 1 ðJðo3 ÞÞg.

ð48Þ

Since the non-linearity is on the output, the Volterra series will be infinite and so an infinite number of GFRFs will be required to completely characterise the system. The successive higher-order GFRFs may be obtained recursively by adopting the routine demonstrated above and constructed using combinations of the lower-order GFRFs. Evaluating the GFRFs above in the fast-sampling limit to give lim H 2 ðJðo1 Þ; Jðo2 ÞÞ ¼ H 2 ðjo1 ; jo2 Þ

T!0

implies that for small sampling intervals, the GFRFs of d-domain models are very similar to that of the continuous-time model. Notice that any d-NARMAX models identified from data may be mapped into the GFRFs quite readily, as shown previously. For example, for a sample frequency os , then a Nyquist frequency oNq (defined as os =2) is the limit for which the discrete and continuous linear FRFs are approximately equal to each other. In the second-order GFRF case, this region is defined by o1 poNq , o2 poNq , o1 þ o2 poNq . These restrictions also hold for the higher dimension cases i.e., oi poNq for all i ¼ 1 . . . n, oi þ oj poNq for all i; j ¼ 1; . . . ; n, iaj, oi þ oj þ ol poNq for all i; j; l ¼ 1; . . . ; n, iajal, etc. Beyond such limits, the GFRFs become repeated or aliased versions at higher frequencies which are meaningless. This repetitive nature of the discrete-time GFRFs is due to the periodicity of ½ejoT
1T
1 in the frequency response function. These limits are best viewed graphically as depicted in Fig. 1 for the second-order case. In Fig. 1 it can be seen that the limits on the frequency range for which the discrete GFRFs are meaningful form a prism or ‘wedge’ shape. As the sampling interval T s is reduced, the wedge expands to cover more and more of the three-dimensional space; depicted in Fig. 1 by the shaded cross-sections. In the simulation example, the following coefficients m ¯ ¼ 1, c¯ ¼ 20, k¯ 1 ¼ 104 , k¯ 2 ¼ 107 and k¯ 3 ¼ 9 5  10 for the d-domain model of Eq. (32) were

(49)

Ts

Gain

H 2 ðJðo1 Þ; Jðo2 ÞÞ ¼

3253


2


1 +
2 ≤
Nq


2 ≤
Nq

or more generally lim H n ðJðo1 Þ; . . . ; Jðon ÞÞ ¼ H n ðjo1 ; . . . ; jon Þ

T!0


1 ≤
Nq

(50)

which indicates that as the sampling interval tends towards zero, the linear and non-linear frequency response functions of d-domain models converge to the continuous-time equivalent. In practise, this


1 Fig. 1. Meaningful frequency range for discrete second-order GFRFs.

ARTICLE IN PRESS M.A. Chadwick et al. / Signal Processing 86 (2006) 3246–3257

3254

substituted to produce the following form:

600

d2 yðkÞ þ 20dyðkÞ þ 104 yðkÞ þ 107 y2 ðkÞ

400

ð51Þ

which may be viewed as a direct analogue of that found in [24] in the d-domain but with the continuous-time coefficients. The same system in the continuous domain with the same parameters, produces a resonant frequency of or ¼ 99 rad s
1 or f r ¼ 15:75 Hz. Comparing the GFRFs of the continuous model and the discrete delta model generated using the algorithm in Section 6 (see Figs. 2, 3 and 4), it can be seen that they become more and more alike as the sample interval is reduced. This indicates the ability of the d-operator to approximate a derivative also extends to the non-linear case and transform domains demonstrating the usefulness and capability of this method. The features in Fig. 3 indicate significant nonlinear effects along the ridges in H 2 ðÞ. The lines of the ridges for large magnitudes indicate that if an input signal has two components at frequencies which satisfy the relationship, ao1 þ bo2 ¼ or , then the non-linear model terms transfer energy to excite the system at the resonant frequency. The values of a and b are determined by the position in o1 –o2 space of the resonant frequency ridge which in turn is governed by the properties of the system. This implies that a resonance can be excited even though the input signal may not have energy at the resonance frequency itself. This is contrary to linear theory where a resonance can only be excited by an −60 −70

Gain (dBs)

−80 −90 −100 −110 −120 −130

200
2

þ 5  10 y ðkÞ ¼ uðkÞ

0 −200 −400 −600 0

100

(a)

Gain

9 3

200

300
1

400

500

600

Contour plot

−80 −100 −120 −140 −160 −180 −200 0 200
1

(b)

500 400

600 −500 Actual magnitude plot

0


2

Fig. 3. Magnitude plots of H 2 ðjo1 ; jo2 Þ. (a) Contour plot. (b) Actual magnitude plot.

input signal with energy at the frequency of the resonance. The third-order GFRFs cannot be plotted directly but methods have been developed to visualise the effects of H 3 ðÞ [31]. One approach is to plot two-dimensional submanifolds of H 3 ðÞ by fixing one frequency to a resonance frequency, a constant or even the same as one of the other frequencies. Additional methods are discussed in [24]. These plots reveal interesting and valuable information about the non-linear systems dynamics and can be directly produced from models of the ddomain description. 7. Conclusions

0

200

400

600 800
(rad/s)

1000 1200 1400

Fig. 2. Magnitude plot of H 1 ðjoÞ and H d1 ðjoÞ for F s ¼ 100; 200 and 3200 Hz.

A new framework for identifying and analysing fast-sampled non-linear systems has been presented. Classical methods employing the absolute-positional

ARTICLE IN PRESS M.A. Chadwick et al. / Signal Processing 86 (2006) 3246–3257

600

Appendix

400

For a sampled-data degree-n homogeneous system, the output (in the d-domain) can be calculated at regularly spaced sampling instants by means of a weighted sum of input values.


2

200 0 −200

yn ðkÞ ¼

1 X



i1 ¼0

−400

1 X

hdn ði1 ; . . . ; in Þ

in ¼0

n Y

uðk
im ÞT,

m¼1

k ¼ 0; 1; 2; . . . ,

−600 0

100

(a)

200

300
1

400

500

600

Contour plot

−80 −100 Gain

3255

−120 −140 −160 0 200
1

(b)

500 400

−500 600 Actual magnitude plot

0


2

ðA1Þ

where the kernel hdn ði1 ; . . . ; in Þ is a real sequence and equal to zero if any argument is negative (onesided). The superscript d notation is used to indicate the kernel above is derived by discretising the continuous Volterra series in accordance with the d-domain procedure and is different from the traditional q-domain process used throughout existing literature. As with the continuous case, it is desirable to find a similar form of input–output relationship such as Y ðgÞ ¼ HðgÞUðgÞ which exists in the linear case. However, due to the dimensional properties there is no direct way to do this. An indirect way is adopted by amending Eq. (A1) and introducing an auxiliary multi-dimensional time function yn ðk1 ; . . . ; kn Þ ¼

1 X



i1 ¼0 n Y

H d2 ðjo1 ; jo2 Þ

Fig. 4. Magnitude plots of with F s ¼ 200 Hz. (a) Contour plot. (b) Actual magnitude plot.



1 X

hdn ði1 ; . . . ; in Þ

in ¼0

uðkm
im ÞT,

ðA2Þ

m¼1

q operator were found to be less suited to fastsampled systems than the d-operator. A promising method marrying together the theories of NARMAX and d-operators, to exploit both their benefits, has been illustrated. The method is simple in principle and easily integrates into current methodologies. Current theories based on the linear D-transform were extended to the nonlinear case and applied to discrete-time Volterra models. An example using a Duffing equation was used to motivate and illustrate the new method.

where the real output is recovered by the restriction yn ðkÞ ¼ yn ðk1 ; . . . ; kn Þjk1 ¼¼kn ¼k .

If we take a multiple D transform of both sides of the above equation we get Y n ðg1 ; . . . ; gn Þ ( ) 1 1 n X X Y d  hn ði1 ; . . . ; in Þ uðkm
im ÞT ¼ Dn i1 ¼0

¼

1 X



k1 ¼0 n Y

Acknowledgements The authors gratefully acknowledge that this work was partly supported by EPSRC (UK).

(A3)



m¼1

i1 ¼0

1 X 1 X

m¼1



kn ¼0 i1 ¼0

uðkm
im ÞT

1 X in ¼0 n Y

hdn ði1 ; . . . ; in Þ

½1 þ gm T
km T.

ðA4Þ

m¼1

Noting the common domain of the product terms Qn Qn
km Eq. (A4) m¼1 uðkm
im Þ and m¼1 ½1 þ gm T

ARTICLE IN PRESS M.A. Chadwick et al. / Signal Processing 86 (2006) 3246–3257

3256

reduces to

time Fourier transformation definition.

Y n ðg1 ; . . . ; gn Þ ¼

1 X



1 X 1 X



kn ¼0 i1 ¼0

k1 ¼0 n Y

1 X

hdn ði1 ; . . . ; in Þ

in ¼0

uðkm
im Þ½1 þ gm T
km T 2



m¼1

ðA5Þ In order to achieve the desired form of Eq. (18) we make a substitution of variables k^m ¼ km
im where km ¼ k^m þ im produces the following: Y n ðg1 ; . . . ; gn Þ ¼

1 X



1 X 1 X kn ¼0 i1 ¼0

k1 ¼0 n Y



1 X

hdn ði1 ; . . . ; in Þ

in ¼0 ^

uðk^m Þ½1 þ gm T
ðkm þim Þ T 2



m¼1

ðA6Þ since the terms involving km are independent of the terms involving im this can be rearranged and grouped as follows: Y n ðg1 ; . . . ; gn Þ ¼

1 X



i1 ¼0 n Y



1 X

½1 þ gm T
im T

1 X

1 X



k^1 ¼0 n Y

^

uðk^m Þ½1 þ gm T
km T ðA7Þ

which, assuming causality (i.e., hn ðk1 ; . . . ; kn Þ ¼ uðki Þ ¼ 0 8ki o0Þ and the definition of H n ðg1 ; . . . ; gn Þ gives the following: Y n ðg1 ; . . . ; gn Þ ¼ H n ðg1 ; . . . ; gn ÞUðg1 Þ . . . Uðgn Þ. (A8) Notice that by using Dirac delta functions (not to be confused with the d-operator) the sampled form of the sequence hdn ði1 ; . . . ; in Þ can be expressed as a time function as 1 X i1 ¼0


1 i1 ¼0

in ¼0

dðt1
i1 TÞ    dðtn
in TÞ e
jðo1 t1 þþon tn Þ dt1 . . . dtn 1 1 X X ¼  hdn ði1 T; . . . ; in TÞ i1 ¼0

e

in ¼0


jðo1 i1 þþon in ÞT

.

ðA10Þ

Therefore the Fourier transfer function can be viewed as H n ðg1 ; . . . ; gn Þ evaluated on hypersphere jg1 j ¼    ¼ jgn j ¼ 1=T centred at
1=T where T is the sampling interval. The Fourier transfer function is then given by the mathematical relationship H dn ðjo1 ; . . . ; jon Þ ¼ H n ðg1 ; . . . ; gn Þjg1 ¼ðejo1 T
1Þ=T;...;gn ¼ðejon T
1Þ=T . ðA11Þ References

k^n ¼0 m¼1

hdn ðt1 ; . . . ; tn Þ ¼


1

in ¼0

m¼1



hdn ði1 ; . . . ; in Þ

H dn ðjo1 ; . . . ; jon Þ Z 1X Z 1 1 1 X   hdn ði1 T; . . . ; in TÞ ¼



1 X

hdn ði1 T; . . . ; in TÞ

in ¼0

dðt1
i1 TÞ    dðtn
in TÞ,

ðA9Þ

where T is the sampling interval for the discrete system and the dðti
ii TÞ are Dirac delta functions. The Fourier transform of hdn ðt1 ; . . . ; tn Þ is then obtained from the continuous-

[1] H. Nyquist, Certain topics in telegraph transmission theory, Trans. AIEE 47 (1928) 617–644. [2] C.E. Shannon, Communication in the presence of noise, Proc. Inst. Radio Eng. 37 (1) (1949) 617–644. [3] H. Zhang, S.A. Billings, Analysing nonlinear systems in the frequency domain. I: the transfer function, Mech. Syst. Signal Process. 7 (6) (1993) 531–550. [4] S.A. Billings, L.A. Aguirre, Effects of the sampling time on the dynamics and identification of nonlinear models, Internat. J. Bifurcation Chaos 5 (6) (1995) 1541–1556. [5] R.M. Goodhall, B.J. Donoghue, Very high sample rate digital filters using the d operator, IEE Proc.-G 140 (3) (1993) 199–206. [6] R.H. Middleton, G.C. Goodwin, Digital Control and Estimation: A Unified Approach, Prentice-Hall, Englewood Cliffs, NJ, 1990. [7] H.H. Fan, P. De, High speed adaptive signal processing using the delta operator, Digital Signal Process. 11 (1) (2001) 3–34. [8] G.C. Goodwin, R.H. Middleton, High-speed digital signal processing and control, Proc. IEEE 80 (2) (1992) 240–259. [9] H. Fan, T. So¨derstro¨m, M. Mossberg, B. Carlsson, Y. Zou, Estimation of continuous-time AR process parameters from discrete-time data, IEEE Trans. Signal Process. 47 (5) (1999) 1232–1244. [10] T. So¨derstro¨m, H. Fan, B. Carlsson, S. Bigi, Least squares parameter estimation of continuous-time ARX models from discrete-time data, IEEE Trans. Automat. Control 42 (5) (1997) 659–673. [11] A. Khodabakhshian, K. Jamshidi, A new identification model for power system transfer functions, Electric Power Components Systems 31 (5) (2003) 513–524.

ARTICLE IN PRESS M.A. Chadwick et al. / Signal Processing 86 (2006) 3246–3257 [12] M. Rostgaard, M.B. Lauritsen, N.K. Poulsen, A state-space approach to the emulator-based gpc design, Systems Control Lett. 28 (5) (1996) 291–301. [13] A.G. Kuznetsov, R.O. Bowyer, D.W. Clarke, Estimation of multiple order models in the d-domain, Internat. J. Control 72 (7/8) (1999) 629–642. [14] E.G. Collins Jr., T. Song, A delta operator approach to discrete-time H1 control, Internat. J. Control 72 (4) (1999) 315–320. [15] G.C. Goodwin, R.L. Leal, D.Q. Mayne, R.H. Middleton, Rapprochment between continuous and discrete model reference adaptive control, Automatica 22 (2) (1986) 199–207. [16] C.P. Neuman, Properties of the delta operator model of dynamic systems, IEEE Trans. Systems Man Cybernet. 23 (1) (1993) 296–301. [17] K.I. Kim, E.J. Powers, A digital method of modelling quadratically nonlinear systems with a general random input, IEEE Trans. Acoust. Speech Signal Process. 36 (1988) 1758–1769. [18] S.A. Billings, K.M. Tsang, Spectral analysis for nonlinear systems. Part I: parametric nonlinear spectral analysis, Mech. Systems Signal Process. 3 (4) (1989) 319–339. [19] S.A. Billings, K.M. Tsang, Spectral analysis for nonlinear systems. Part II: interpretation of nonlinear frequency response functions, Mech. Systems Signal Process. 3 (4) (1989) 341–359. [20] S.A. Billings, K.M. Tsang, G.R. Tomlinson, Spectral analysis for nonlinear systems. Part III: case study examples, Mech. Systems Signal Process. 4 (1) (1990) 3–21. [21] R.H. Middleton, G.C. Goodwin, Improved finite word length characteristics in digital control using delta operators, IEEE Trans. Automat. Control 31 (11) (1986) 1015–1021.

3257

[22] W.J. Rugh, Nonlinear System Theory: The Volterra/Weiner Approach, John Hopkins University Press, Baltimore, MD, 1981. [23] M. Schetzen, The Volterra and Weiner Theories of Nonlinear Systems, Wiley, New York, 1980. [24] K. Worden, G.R. Tomlinson, Nonlinearity in Structural Dynamics: Detection, Identification and Modelling, Institute of Physics, 2001. [25] I.J. Leontaritis, S.A. Billings, Input–output parametric models for nonlinear systems. Part I: deterministic nonlinear systems, Internat. J. Control 41 (2) (1985) 303–328. [26] I.J. Leontaritis, S.A. Billings, Input–output parametric models for nonlinear systems. Part II: stochastic nonlinear systems, Internat. J. Control 41 (1985). [27] S.A. Billings, K.M. Tsang, G.R. Tomlinson, Application of the narmax method to nonlinear frequency response estimation, in: Sixth International Modal Analysis Conference, 1988. [28] S.A. Billings, K.L. Lee, A smoothing algorithm for nonlinear time series, Internat. J. Bifurcation Chaos 14 (3) (2004) 1037–1051. [29] E. Bedrosian, S.O. Rice, The output properties of volterra systems (nonlinear systems with memory) driven by harmonic and gaussian inputs, IEEE Proc. 59 (1971) 1688–1707. [30] J.C. Peyton-Jones, S.A. Billings, Recursive algorithm for computing the frequency response of a class of non-linear difference equation models, Internat. J. Control 50 (5) (1989) 1925–1940. [31] R. Yue, S.A. Billings, Z.-Q. Lang, An investigation into the characteristics of nonlinear frequency response functions, part 1: understanding the higher dimensional frequency spaces, Internat. J. Control (78) (13) (2005) 1031–1044.