PREPRINT VERSION. ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT

1

Best Linear Approximation of Wiener Systems
Using Multilevel Signals: Theory and Experiments

arXiv:1710.07067v1 [eess.SP] 19 Oct 2017

A. De Angelis, J. Schoukens, K. R. Godfrey, P. Carbone

Abstract—The problem of measuring the best linear approximation of a nonlinear system by means of multilevel excitation
sequences is analyzed. A comparison between different types of
sequences applied at the input of Wiener systems is provided
by numerical simulations and by experiments on a practical
circuit including an analog filter and a clipping nonlinearity.
The performance of the sequences is compared with a white
Gaussian noise signal for reference purposes. The theoretical
characterization of the best linear approximation when using
randomized constrained sequences is derived analytically for the
cubic nonlinearity case. Numerical and experimental results show
that the randomized constrained approach for designing ternary
sequences has a low sensitivity to both even and odd order
nonlinearities, resulting in a response close to the actual response
of the underlying linear system.
Index Terms—Nonlinear systems, best linear approximation,
ternary sequences, binary pseudorandom sequences.

I. I NTRODUCTION
The measurement of the frequency response function (FRF)
of a linear dynamical system is a fundamental step that is
typically performed in engineering applications for modeling
or control purposes. Such procedure is typically performed by
exciting the system under test with a properly designed input
signal and measuring the output. However, many practical
systems are affected by some degree of nonlinear distortion.
For such systems, the FRF will vary as a function of the
excitation that the system is subject to. To analyze these cases,
linearization is often a viable option. For this reason, the concept of best linear approximation (BLA) of nonlinear systems
is of great importance to practical applications. Therefore, the
BLA, which is obtained by solving a least squares problem
where the mean squared difference between the actual output
of the system and the output of a linear model is minimized,
has been thoroughly studied in the literature [1]–[4].
An analysis of the BLA for several cases of the amplitude
distribution of the input signal is developed in [5], considering
Gaussian signals and binary signals. The analysis is further
This work was supported in part by the Italian Ministry of Instruction,
University and Research (MIUR) through grant PRIN 2015C37B25.
A. De Angelis and P. Carbone are with the Engineering Department, University of Perugia, via G. Duranti 93, 06125 Perugia, Italy. {alessio.deangelis,
paolo.carbone}@unipg.it
J. Schoukens is with Department ELEC, Vrije Universiteit Brussel, Pleinlaan 2, B1050 Brussels, Belgium.
K. R. Godfrey is with the School of Engineering, University of Warwick,
CV4 7AL, Coventry, United Kingdom.
c 2017 IEEE. Personal use of this material is permitted. Permission from
IEEE must be obtained for all other uses, in any current or future media,
including reprinting/republishing this material for advertising or promotional
purposes, creating new collective works, for resale or redistribution to servers
or lists, or reuse of any copyrighted component of this work in other works.

extended in [6] to multilevel signals. A particularly useful class
of multilevel signals is that of ternary sequences, which are
easy to generate in practical applications where the number of
available levels is limited, and allow for defining a wider range
of harmonic behavior when compared to binary sequences.
In this context, a direct synthesis (DS) method for obtaining
ternary sequences with harmonic multiples of two and three
suppressed is proposed in [7]. Such sequences are well suited
for the characterization of even and odd nonlinearities in
dynamical systems and for the mitigation of the impact of
nonlinear distortion on the FRF measurement. Furthermore,
in [8], [9], randomized constrained sequences (RCS) are introduced, allowing for a more flexible approach to design ternary
sequences with defined spectral properties. The RCS approach
is particularly suited for practical measurement applications,
since it provides robustness in the presence of non-idealities
such as nonuniform digital-to-analog converter (DAC) levels.
In this paper, the study of the properties of RCS and DS
ternary sequences is developed by analyzing their behavior
when used as excitation signals to measure the BLA of a
nonlinear system. The performance is evaluated on a blockstructured nonlinear system and compared to commonly-used
binary sequences and Gaussian noise. The latter is assumed
as a reference, due to its well-studied properties for BLA
measurement [10].
The BLA is characterized theoretically, extending the numerical simulation results in [11]. Specifically, the analytical
expression of the cross-correlation between the input and the
output of Wiener systems with cubic nonlinearity and RCS
input is derived. Furthermore, the behavior of the BLA is
validated by experimental measurements performed on the
prototype of a Wiener system implemented by an analog filter
followed by a diode-based clipper. It is shown that the RCS
approach is relatively more robust to nonlinearities compared
to the other considered approaches, thus providing a BLA
closer to the underlying linear system.
II. T HE BLA OF A W IENER N ONLINEAR S YSTEM
The BLA of a nonlinear dynamical system is defined as
a linear system such that its impulse response minimizes the
mean square value of an error sequence e [·] defined as [12]:
e [k] = y [k] − (gBLA ∗ u) [k]

(1)

where gBLA [·] is the impulse response of the BLA, u [·] is
the input of the nonlinear system, y [k] is its output at time
instant k, and ∗ denotes the convolution operation.
When the system under test is subject to different excitation signals, a different BLA is obtained, depending on the

PREPRINT VERSION. ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT

Fig. 1. Diagram of a Wiener system (top) and its BLA (bottom). The symbol
g[k] denotes the impulse response of a linear system and the symbol f (x)
denotes a static nonlinear transformation acting on the sequence x.

distribution of the amplitude of the excitation and its power
spectrum. For a Wiener system, defined as the cascade of a
linear dynamical system and a static nonlinearity as shown in
Fig. 1, Gaussian excitations result in a BLA that is a scaled
version of the underlying linear system [5].
This property represents a significant advantage in measurement applications, especially those where the dynamics of the
underlying linear system must be measured while mitigating
the impact of nonlinear distortions [13]. However, the usage of
Gaussian signals is not always advisable in practical scenarios.
In fact, other classes of excitations are characterized by useful
properties, such as a lower crest factor, resulting in more power
into the system under test for the same maximum amplitude.
Furthermore, in some applications, there is a limitation on
the number of levels in the input signal. Such constraints
may be related to the particular measurement process, ease of
implementation considerations, or design of the actuator [14].
Therefore, in those applications, one may resort to the use of
binary sequences or multilevel signals. The latter are defined as
discrete sequences with a number of amplitude levels greater
than two [6]. Ternary signals, which use three amplitude levels,
are an example of a multilevel signal.
The use of signals with a limited number of levels implies
that the distribution of the input deviates from Gaussianity.
This causes the BLA to be biased with respect to the Gaussian
BLA. This, in turn, means that the BLA will not be a scaled
version of the linear system response. An in-depth theoretical
study of such modifications is presented in [5], which provides analytical expressions for the BLA using common nonGaussian signals for general Wiener–Hammerstein systems.
III. E XCITATION SIGNALS FOR MEASURING THE BLA
In this section, we provide a description of several signals
that will then be used in the next section for numerical
comparison of the BLA of a Wiener system. The aim is
to prove the applicability of such signals to the purpose of
measuring the BLA, and to highlight the advantages of ternary
sequences with respect to commonly-used binary sequences.
The excitation signals considered in this paper are ternary
sequences, specifically DS and RCS, pseudorandom binary
sequences, specifically maximum length binary sequences
(MLBS) and inverse-repeat MLBS (IRMLBS), together with
white Gaussian noise (WGN).
In particular, DS sequences are constructed according to
the analytical method introduced in [7]. A DS sequence

2

is obtained starting from a basic sequence with specified
harmonic properties, which can be a MLBS among other
classes of sequences. The basic sequence is repeated n times
and multiplied by a special fixed sequence. For suppression
of harmonic multiples of 2 and 3, n = 6. Due to this
synthesis mechanism, the length N of a MLBS-based DS must
satisfy the constraint N = 6(22k+1 − 1), with k integer. This
requirement implies a lack of flexibility in terms of available
periods. On the other hand, the advantage of the DS approach
is that harmonic multiples of two and three are suppressed,
which removes the effects of even order nonlinearities and
mitigates those of odd order ones [7], [15].
Recently, another approach, i.e. Randomized Constrained
Sequences (RCS), for generating ternary sequences with harmonic multiples of two and three suppressed has been proposed in [8], [9]. Such approach is based on random sequences
that are constrained to satisfy the harmonic suppression conditions and are iteratively modified in order to approach a design
goal, such as uniformity of the levels of excited frequencies.
The RCS method allows for a greater flexibility in terms of
achievable sequence lengths, since it can provide sequences of
length 6k, with k integer.
For comparison purposes, MLBS are considered in this
paper, since these binary pseudorandom sequences are commonly used for measuring the response of dynamical systems.
These sequences are generated using linear-feedback shift
registers with properly chosen feedback terms [16]. As a
consequence, the available MLBS sequence lengths are 2k −1,
with k integer. Furthermore, there is a limited number of
unique sequences for a given length. Notice that it is not
possible to suppress odd-order harmonics using MLBS (given
that they have only two levels). All harmonics are excited
when using an MLBS, but even harmonic suppression can be
achieved by IRMLBS, obtained by inverting every other digit
in an MLBS of period N to obtain an IRMLBS of period
2N [17]. For this reason, IRMLBS are also added to the
comparison.
IV. N UMERICAL COMPARISON
For numerical simulation, four Wiener systems are considered, with the parameters shown in Table I. In particular, the
linear block of the Wiener system consists of a FIR filter
having impulse response coefficients given by [1 0.7 0.3] in
the case of memory length 3, and [1 0.7 0.3 0.2 0.1 0.05] in
the case of memory length 6. Furthermore, for the nonlinear
block of the system, two nonlinearities are considered: x3 and
x3 + x2 . For each system, each of the considered sequences
is applied at the input, u[k], and the output sequence y[k]
is computed at steady state, by discarding the transient. The
system is noise free, because we are interested in studying the
effect of the sequence distribution and not of measurement
noise. The sequence length is 762 for DS, RCS, and WGN,
511 for the MLBS case, and 510 for the case of IRMLBS.
Different lengths are used because of the sequence-specific
length limitations explained in the previous section. As an
example, the length of DS must be 6(22k+1 − 1), which is
an even number. On the other hand, the MLBS must be of

PREPRINT VERSION. ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT

DS;

RCS;

WGN;

IRMLBS;

linear;

10

(a)

ĜBLA [dB]

ĜBLA [dB]

10

MLBS;

0
−10
−20

(b)

0
−10
−20

0

0.1

0.2

0.3

0.4

0.5

0

0.1

(c)

0
−10
−20
−30
0

0.1

0.2

0.3

0.3

0.4

0.5

0.4

(d)

0
−10
−20
−30

0.5

0

0.1

Normalized frequency

0.2

0.3

0.4

0.5

Normalized frequency

10

10

(e)

ĜBLA [dB]

ĜBLA [dB]

0.2

Normalized frequency
standard dev. [dB]

standard dev. [dB]

Normalized frequency

0
−10
−20

(f)

0
−10
−20

0

0.1

0.2

0.3

0.4

0.5

0

0.1

Normalized frequency

(g)

0
−10
−20
−30
0

0.1

0.2

0.3

0.2

0.3

0.4

0.5

Normalized frequency
standard dev. [dB]

standard dev. [dB]

3

0.4

0.5

(h)

0
−10
−20
−30
0

0.1

Normalized frequency

0.2

0.3

0.4

0.5

Normalized frequency

Fig. 2. Numerical simulation results. The absolute value and the standard deviation of the estimated BLA are shown for the systems defined in Table I: (a)
System 1; (b) System 2; (c) System 1 standard deviation; (d) System 2 standard deviation; (e) System 3; (f) System 4; (g) System 3 standard deviation; (h)
System 4 standard deviation.

the form 2k − 1, which is an odd number. Therefore, it is not
possible to generate a DS sequence and an MLBS sequence
having the same length. The sequence lengths were selected
to be as close as possible. First, a value of 762 was chosen
arbitrarily for the sequence length of DS. This choice allowed
the computational complexity of further processing operations
to be contained. Once a length of 762 was defined for DS,
the closest MLBS and IRMLBS lengths were chosen, namely
511 for MLBS and 2×255 = 510 for IRMLBS.
The BLA is estimated as follows [12]:
ĜBLA (jωk ) =

ŜY U [k]
ŜU U [k]

(2)

where ŜY U [k] is the cross-power spectrum of the input and
output sequences, and ŜU U [k] is the autopower spectrum of the
input. The estimates are obtained by first dividing the sequence
realizations into groups of 4 unique sequences each. Then,

TABLE I
S UMMARY OF CONSIDERED SYSTEMS .
Impulse response g[·]

Nonlinearity f (·)

System 1

[1 0.7 0.3]

x3

System 2

[1 0.7 0.3 0.2 0.1 0.05]

System 3

[1 0.7 0.3]

System 4

[1 0.7 0.3 0.2 0.1 0.05]

x3
x3

+ x2

x3 + x2

each group is used to provide an estimate of the BLA by
averaging before division. Such averaging significantly reduces
the variance increase due to the random behavior for Gaussian
signals [12]. The mean and standard deviation of the BLA
estimates provided by each group are then calculated. Due to
the limited availability of unique MLBSs for a given length,
as described in [18], 4 groups are formed for the DS and

PREPRINT VERSION. ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT

TABLE II
K ULLBACK -L EIBLER DIVERGENCE BETWEEN THE EMPIRICAL
PROBABILITY DISTRIBUTIONS , ESTIMATED BY HISTOGRAMS , OF THE
CONSIDERED SEQUENCES AND THE IDEAL NORMAL DISTRIBUTION , AT
THE OUTPUT OF THE LINEAR BLOCK OF S YSTEM 2.
DS

RCS

WGN

MLBS

IRMLBS

0.2346

0.0654

0.0255

0.1837

0.1897

IRMLBS cases, 12 groups for the MLBS case, and 100 groups
for the other cases.
The input sequences are normalized so that they have the
same RMS value. The normalization is performed such that
the impulse response of the Gaussian BLA is equivalent to
the true underlying linear impulse response, according to the
procedure in [5]. Additionally, the
p ternary sequences RCS and
DS are multiplied by the term 3/2, in order to obtain the
same power for all excitations.
Numerical results for all considered systems and input
sequence types are shown in Fig. 2, where only excited
harmonics are displayed. In the figure, scaled versions of the
BLA are plotted. The scaling factor is chosen such that the
least square errors between each BLA and the FRF of the
linear filter are minimized. This choice has been made because
we are interested in the dissimilarity between the dynamics.
The standard deviation has been scaled by the same factor.
In Fig. 3, the ratio between the linear system response
and the BLA is plotted, with the same scaling as in Fig.
2, from which it can be seen that the RCS input signal
provides the closest response to the true underlying linear
system response, except for Gaussian noise. This is because
the amplitude distribution of the RCS at the output of the linear
system is the closest to a Gaussian distribution, among the
considered signals. The latter behavior is also highlighted by
the Kullback-Leibler divergence [19] values listed in Table II.
V. D ISCUSSION OF THE RESULTS
Results show that odd ternary sequences are less sensitive
to the effect of even-order nonlinearities and memory length
of the linear system. In fact, they show less variability when
compared to white Gaussian noise excitation and to binary
signals that excite all harmonics (such as MLBS). This can be
noticed by observing that the standard deviation for MLBS in
Fig. 2(g), where the x2 term is present, is considerably larger
than that of Fig. 2(c).
Moreover, from the results of Fig. 3, it can be noticed
that DS is very sensitive to nonlinear distortions, due to its
particular construction. On the other hand, RCS, thanks to
its randomization properties, is less sensitive to nonlinearity.
Still, RCS maintains the ability to discriminate even and
odd distortions, due to harmonic multiples of two and three
suppressed. This is not possible with MLBS and IRMLBS.
VI. T HEORY FOR CUBIC NONLINEARITY AND RCS INPUT
Motivated by the numerical simulation results presented
in the previous sections, we analyze the performance of
RCS for Wiener systems also from a theoretical point of
view. Specifically, in this section, we provide a theoretical

4

characterization of the behavior of RCS when applied to the
input of a Wiener system with a cubic nonlinearity. Notice
that a cubic nonlinearity may be used to model numerous
real-world nonlinear systems often encountered in practice,
which can be represented using odd functions. The theoretical
characterization is performed by deriving the cross-correlation
between the input and the output, as stated in the following
Theorem.
Theorem 1. Consider an RCS of length N as defined in [9]
at the input of a causal Wiener system. Let the Wiener system
consist of a FIR filter of order H, with H < N/6, and a cubic
nonlinearity. Denote the input of the Wiener system by u[·], its
output by y[·], the autocorrelation of the input by Ru [·], and
the impulse response of the FIR filter by g[·].
Then, the cross-correlation between the output and the input
is given by

Ryu [r] = 2α2 (g ∗ Ru ) [r] − g 3 ∗ Ru [r]
(3)

2
(4)
= g [rm ] Ru [rq ] 2α2 − g [rm ] ,
r = 0...N − 1,
H
P

g 2 [k], rq denotes the quotient of the division
k
N
r
r
,
i.e.
r
,
q
N/6
N/6
6 , and rm denotes the reminder, with
rm = r − rq .
where α2 ,

k=0j

Proof: See Appendix I.
Note that the expressions for the cross-correlation in (3)
and (4) apply also to Wiener systems where the nonlinearity
is y = x3 + x2 . This is due to the general result that zeromean excitations with symmetric distribution, combined with
even-order nonlinearities (with even-degree Volterra kernels),
have a zero BLA [12].
It is interesting to compare the derived expression of the
cross-correlation to that for white Gaussian input sequences
of variance one, which is given by [5]:
Ryu(Gaussian) [r] = 3α2 g[r] .

(5)

If the cross-correlation in (3) is evaluated for |r| < N/6 and
for an RCS input
p of unitary variance, i.e. the RCS sequence
is scaled by 3/2, it reduces to
3
Ryu [r] = 3α2 g[r] − g 3 [r].
(6)
2
By comparing (6) to (5), it can be noticed that the bias term of
the RCS with respect to the Gaussian input case is −3/2g 3 [r],
whereas for random binary sequences it is −2g 3 [r], which
applies also to MLBS in the limit when N is large [5].
Therefore, the BLA of the RCS is closer to the underlying
linear system than that of the MLBS.
The theoretical results in Theorem 1 have been validated by
numerical simulations. In particular, 103 different instances
of an RCS of length N = 762 were generated and fed to
a Wiener system with H = 3, g = [1, 0.7, 0.3], and a
cubic nonlinearity. For each input sequence, the corresponding
output was stored. Subsequently, the cross-correlation between
the input and the output was calculated for each sequence.
Finally, the average of all cross-correlation sequences was

PREPRINT VERSION. ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT

RCS;

WGN;

15

MLBS;

(a)

10
5
0

IRMLBS;

G/ĜBLA [dB]

G/ĜBLA [dB]

DS;

−5

Linear;

15

(b)

10
5
0
−5

0

0.1

0.2

0.3

0.4

0.5

0

0.1

0.2

0.3

0.4

0.5

Normalized frequency

15

(c)

10
5
0
−5

G/ĜBLA [dB]

Normalized frequency
G/ĜBLA [dB]

5

15

(d)

10
5
0
−5

0

0.1

0.2

0.3

0.4

0.5

0

0.1

0.2

Normalized frequency

0.3

0.4

0.5

Normalized frequency

Fig. 3. Numerical simulation results. Ratio between linear system response and BLA. (a) System 1; (b) System 2; (c) System 3; (d) System 4.

simulation

theory, eq. (3)

(a)

1

Ryu (r)

100 nF

theory, eq. (4)

1 kΩ

0

−
1 kΩ
VI

−1
0

200

400

180 Ω

+
VO

600

r

Ryu (r)

1.5

(b)
Fig. 5. Schematic of the realized Wiener system, consisting of a low pass
filter and a soft clipper.

1
0.5
0
0

2

4

6

8

10

r
Fig. 4. Comparison between numerical simulations and theoretical equations
(3) and (4). (a): full cross-correlation sequence. (b): magnification of the first
peak of the cross-correlation sequence.

calculated. Such average cross-correlation is depicted in Fig.
4, together with the theoretical values obtained from (3) and
(4). As shown in Fig. 4, there is a very good agreement
between theoretic results and simulations, thus validating the
derivations of Theorem 1. In the following section, we show
that the behavior analyzed theoretically for cubic nonlinearity
may be observed also in a practical situation, with a realistic
type of nonlinearity.
VII. E XPERIMENTAL R ESULTS ON A R EAL -W ORLD
W IENER S YSTEM
To validate the theoretical results presented in Sections IV VI, experiments were performed on a real-world Wiener

system, implemented as an electronic circuit. The considered
Wiener system is comprised of a first-order low-pass active
filter, having a cutoff frequency of approximately 1600 Hz,
followed by a diode-based clipper with a series resistor, thus
performing soft clipping, as depicted in the schematic of
Fig. 5. This circuit was first analyzed by numerical simulations, then experimentally tested.
A. Simulation
Numerical simulations were carried out by modeling the
lowpass filter as a first-order iir filter and performing hard
clipping numerically. Hard clipping is an approximation of
the soft clipping operation of the actual circuit. The same
simulation procedure as in Section IV was used. Specifically,
due to a limited number of sequences available, 16 instances of
DS and IRMLBS, 48 instances of MLBS, and 1000 instances
of RCS and WGN were used. Oversampling by a factor of 10
was employed, to mitigate the effect of frequency warping
close to Nyquist frequency due to the iir filter. Therefore,
the sampling frequency of the data acquisition system was
200 kSa/s, whereas the signal generator frequency was 20 kHz.

PREPRINT VERSION. ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT

Simulation results are shown in Fig. 6, where plots are scaled
so that the mean squared error between each BLA and the
FRF of the linear filter is minimized.

A PPENDIX A
P ROOF OF T HEOREM 1
The cross-correlation between the output and the input is
Ryu [r] = E [y[t]u[t − r]]


!3
H−1
X
=E
g[k]u[t − k] u[t − r]

B. Experimental Characterization
A prototype of the considered Wiener system shown in
Fig. 5 was built using two 1N4148 diodes and a OP27
operational amplifier. Experiments were conducted using the
“Analog discovery 2” acquisition board by Digilent [20]. In
particular, the same sequences as those used in simulations
were applied at the input, scaled so they had the same desired
RMS amplitude, and the input and output of the system
were measured simultaneously, discarding the initial transient.
The experiment was repeated for several values of the RMS
amplitude. The sampling frequency of the on-board DAC
was set to 20 kSa/s, with a resolution of 14 bits. Moreover,
the sampling frequency of the ADC was set to 200 kSa/s
with a resolution of 14 bits. Note that oversampling was
performed also in the case of the experiments. The purpose
of oversampling in this context was to mitigate the effect of
sharp transitions of the input signal.
As a reference, the frequency response of the underlying
linear system was also measured, using a stepped sine excitation with a low amplitude of 0.1 V that was chosen in order
to avoid exciting the nonlinearity.
Experimental results are shown in Fig. 7, where the averaging process and number of sequence instances used is the same
as in the simulations. By comparing such results with numerical simulations in Fig. 6, a strong agreement in the behavior
of the considered sequences can be noticed. Specifically, the
RCS is closer to the underlying linear system, both in terms
of magnitude and phase, with respect to the other multilevel
sequences compared. This is particularly emphasized when the
nonlinearity is strongly excited.

6

(7)

k=0

By expanding (7), we obtain

H−1
X
Ryu [r] = E 
g 3 [k]u3 [t − k]
k=0

+3

H−1
X H−1
X

g 2 [k]g[h]u2 [t − k]u[t − h]

k=0 h=k+1

+6

H−1
X H−1
X

H−1
X


g[k]g[h]g[j]u[t − k]u[t − h]u[t − j]

k=0 h=k+1 j=h+1


×u[t − r] ,

r = 0...N − 1

(8)

where the terms of the expansion of a polynomial raised to
the third power are given by the multinomial coefficients [21].
Recall from [9] that E [u[n]] = 0 , ∀n and that the autocorrelation of the RCS is
 2
n=0
 3


2

−
 3 n = N2
1
(9)
Ru [n] =
n = N6 , 56 N .
3

1
N 2

 −3 n = 3 , 3N


0
elsewhere
By using (9), and considering that H < N/6, we obtain, for
k, h, j = 0 . . . H − 1, h 6= k, j 6= k, j 6= h, that the random
variables u [t − k], u[t−h], u[t−j], and u[t−r] are statistically
independent and
E [u [t − k] u[t − h]u[t − j]u[t − r]] = 0 .

(10)
3

VIII. C ONCLUSION
The use of ternary sequences for measuring the BLA of a
Wiener nonlinear system was investigated by numerical simulations, theoretical derivations, and experiments. The theoretical characterization of the randomized constrained sequences
was provided. Moreover, the performance of ternary sequences
was compared to that of commonly used binary sequences
and to white Gaussian noise. The comparison was carried out
both numerically and by experiments on a realized Wiener
circuit performing analog filtering and clipping. Results show
that the randomized constrained sequence approach provides
low sensitivity to nonlinearities, resulting in a BLA that is
close to the underlying linear system response. Compared to
binary sequences, it has the advantage of suppressing harmonic
multiples of two and three, allowing a nonparametric nonlinear
analysis to be made without any user interaction. Compared to
other signals that provide such suppression, such as multisines,
RCS is also robust to DAC nonlinearities, while a multisine
requires highly linear DAC for proper operation.

Then, by substituting (10) in (8) and noting that u [·] ≡
u[·] for any ternary sequence assuming values [1 , −1 , 0], the
expression in (8) for the cross-correlation simplifies to:
"H−1
X
Ryu [r] = E
g 3 [k]u[t − k]u[t − r]
k=0

+3

H−1
X H−1
X



g 2 [k]g[h]u2 [t − k]u[t − h]u[t − r]

k=0 h=0
h6=k

=

H−1
X

g 3 [k]Ru [r − k]

k=0
H−1
X H−1
X

+3



g 2 [k]g[h]E u2 [t − k] Ru [r − h]

k=0 h=0
h6=k

=

H−1
X
k=0

g 3 [k]Ru [r − k] + 2

H−1
X H−1
X

g 2 [k]g[h]Ru [r − h]

k=0 h=0
h6=k

(11)

PREPRINT VERSION. ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT

DS;

RCS;

WGN;

MLBS;

IRMLBS;

7

linear

G/ĜBLA [dB]

0.2
0
−0.2
0

0.2

0.4

0.6

0.8



∠ (G) − ∠ ĜBLA [deg]

Frequency [Hz]

0.4

−2
0

0.2

0.4

0.6

0.8

Frequency [Hz]

0.6

0.8

Frequency [Hz]

−2
0.2

0

1

0

0

2

·104


∠ (G) − ∠ ĜBLA [deg]

G/ĜBLA [dB]

4
0.4

1

1
·104

20
0
−20
0

0.2

0.4

0.6

0.8

Frequency [Hz]

·104

1
·104

Fig. 6. Wiener system of Fig. 5, numerical simulation results. The two graphs on the left are obtained with an input RMS level of 1 V, whereas those on
the right with an input RMS level of 2 V, thus providing a stronger excitation of the nonlinearity.
DS;

RCS;

WGN;

MLBS;

IRMLBS;

G/ĜBLA [dB]

1
0
−1
0

0.2

0.4

0.6

0.8



∠ (G) − ∠ ĜBLA [deg]

Frequency [Hz]

−5
−10
0.4

0.6

Frequency [Hz]
Fig. 7.

−2

0

0.2

0.8

0.4

0.6

0.8

Frequency [Hz]

0

0.2

0

1

5

0

2

·104

1



∠ (G) − ∠ ĜBLA [deg]

G/ĜBLA [dB]

2

linear (stepped sine)

1
·104

20

0

0

0.2

0.4

0.6

0.8

Frequency [Hz]

·104

1
·104

Wiener system of Fig. 5, experimental test results for an input RMS level of 1 V (left) and 2 V (right).

where
we

 used the fact that the variance of the RCS is
E u2 [n] = 2/3 , ∀n, and the result that functions of independent random variables are independent.
By using simple algebraic manipulations, (11) can be written in terms of the convolution operator as follows:
H−1
X H−1
X

Ryu [r] = g 3 ∗ Ru [r] + 2
g 2 [k]g[h]Ru [r − h]
k=0 h=0
h6=k
3



= g ∗ Ru [r] + 2

H−1
X
k=0

2

g [k]

! H−1
X
h=0

g[h]Ru [r − h]

−2

H−1
X

g 2 [k]g[k]Ru [r − k]

k=0

= g 3 ∗ Ru [r] + 2


H−1
X

!
g 2 [k] (g ∗ Ru ) [r]

k=0


− 2 g 3 ∗ Ru [r]

= 2α2 (g ∗ Ru ) [r] − g 3 ∗ Ru [r]
Therefore, (3) is proven. Finally, the expression in (4)
follows from (3) by noticing that only one term of the
convolution sums in (3) is nonzero for a given r.

PREPRINT VERSION. ACCEPTED FOR PUBLICATION IN IEEE TRANSACTIONS ON INSTRUMENTATION AND MEASUREMENT

R EFERENCES
[1] C. Evans and D. Rees, “Nonlinear distortions and multisine signals.
I. Measuring the best linear approximation,” in IEEE Trans. Instrum.
Meas., vol. 49, no. 3, pp. 602-609, Jun 2000.
[2] L. Ljung, “Estimating linear time-invariant models of nonlinear timevarying systems.” Eur. J. Control vol. 7, no. 2, pp. 203-219, 2001.
[3] P.M. Mäkilä, “LTI approximation of nonlinear systems via signal
distribution theory,” Automatica, vol. 42, no. 6, pp. 917-928, June 2006.
[4] A. Fakhrizadeh Esfahani, J. Schoukens and L. Vanbeylen, “Using
the Best Linear Approximation With Varying Excitation Signals for
Nonlinear System Characterization,” in IEEE Trans. Instrum. Meas.,
vol. 65, no. 5, pp. 1271-1280, May 2016.
[5] H. K. Wong, J. Schoukens and K. R. Godfrey, “Analysis of Best
Linear Approximation of a Wiener-Hammerstein System for Arbitrary
Amplitude Distributions,” in IEEE Trans. Instrum. Meas., vol. 61, no.
3, pp. 645-654, March 2012.
[6] H. K. Wong, J. Schoukens and K. R. Godfrey, “Design of Multilevel
Signals for Identifying the Best Linear Approximation of Nonlinear
Systems,” in IEEE Trans. Instrum. Meas., vol. 62, no. 2, pp. 519-524,
Feb. 2013.
[7] A. H. Tan, “Direct synthesis of pseudo-random ternary perturbation
signals with harmonic multiples of two and three suppressed,” Automatica, vol. 49, no. 10, pp. 2975-2981, October 2013.
[8] A. De Angelis, J. Schoukens, K. R. Godfrey and P. Carbone, “Practical
synthesis of ternary sequences for system identification,” IEEE Int.
Instrum. Meas. Technology Conference Proceedings, pp. 423-428,
Taipei, May 2016.
[9] A. De Angelis, J. Schoukens, K. R. Godfrey and P. Carbone, ”Practical
Issues in the Synthesis of Ternary Sequences,” in IEEE Trans. Instrum.
Meas., vol. 66, no. 2, pp. 212-222, Feb. 2017.
[10] M. Enqvist and L. Ljung, “Linear approximations of nonlinear FIR
systems for separable input processes,” Automatica, vol. 41, no. 3, pp.
459-473, March 2005.
[11] A. De Angelis, J. Schoukens, K. R. Godfrey and P. Carbone, “Measuring the Best Linear Approximation of Wiener Systems Using Multilevel Sequences”, IEEE Int. Instrum. Meas. Technology Conference
Proceedings, Torino, May 2017.

8

[12] R. Pintelon and J. Schoukens, System identification: a frequency
domain approach, John Wiley and Sons, Inc., Hoboken, NJ, 2nd
edition, 2012.
[13] C. Evans and D. Rees, “Nonlinear distortions and multisine signals.
II. Minimizing the distortion,” in IEEE Trans. Instrum. Meas., vol. 49,
no. 3, pp. 610-616, Jun 2000.
[14] H.A. Barker and K.R. Godfrey, “System identification with multi-level
periodic perturbation signals,” Control Engineering Practice, vol. 7,
no. 6, pp. 717-726, 1999.
[15] J. Schoukens, M. Vaes and R. Pintelon, “Linear System Identification
in a Nonlinear Setting: Nonparametric Analysis of the Nonlinear
Distortions and Their Impact on the Best Linear Approximation,” in
IEEE Control Systems, vol. 36, no. 3, pp. 38-69, June 2016.
[16] K.R. Godfrey, A.H. Tan, H.A. Barker and B. Chong, “A survey of
readily accessible perturbation signals for system identification in the
frequency domain,” Control Engineering Practice, vol. 13, no. 11,
November 2005, pp. 1391-1402.
[17] A. H. Tan, K. R. Godfrey and H. A. Barker, “Design of computeroptimized pseudorandom maximum length signals for linear identification in the presence of nonlinear distortions,” in IEEE Transactions
on Instrumentation and Measurement, vol. 54, no. 6, pp. 2513-2519,
Dec. 2005.
[18] New Wave Instruments (2010): Linear feedback shift registers implementation, m-sequence properties, feedback tables. [Online].
Available: http://www.newwaveinstruments.com/resources/articles/m
sequence linear feedback shift register lfsr.htm
[19] S. Kullback and R. A. Leibler, “On information and sufficiency,” The
Annals of Mathematical Statistics, 22(1), pp. 79-86, 1951.
[20] Digilent Inc., “Analog Discovery 2 Reference Manual,” [Online].
Available: https://reference.digilentinc.com/reference/instrumentation/
analog-discovery-2/reference-manual?redirect=1. Retrieved: April
2017.
[21] NIST Digital Library of Mathematical Functions, Section 26.4. http:
//dlmf.nist.gov/26.4, Release 1.0.14 of 2016-12-21. F. W. J. Olver, A.
B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W.
Clark, B. R. Miller and B. V. Saunders, eds.

