City, University of London Institutional Repository
Citation: Siesakul, B. T., Gkoktsi, K. and Giaralis, A. (2015). Compressive power
spectrum sensing for vibration-based output-only system identification of structural systems
in the presence of noise. Compressive Sensing IV, 9484, 94840K. doi: 10.1117/12.2177162
This is the accepted version of the paper.
This version of the publication may differ from the final published
version.
Permanent repository link: https://openaccess.city.ac.uk/id/eprint/13617/
Link to published version: http://dx.doi.org/10.1117/12.2177162
Copyright: City Research Online aims to make research outputs of City,
University of London available to a wider audience. Copyright and Moral
Rights remain with the author(s) and/or copyright holders. URLs from
City Research Online may be freely distributed and linked to.
Reuse: Copies of full items can be used for personal research or study,
educational, or not-for-profit purposes without prior permission or
charge. Provided that the authors, title and full bibliographic details are
credited, a hyperlink and/or URL is given for the original metadata page
and the content is not changed in any way.
City Research Online: http://openaccess.city.ac.uk/ publications@city.ac.uk
City Research Online
*Agathoklis.Giaralis.1@city.ac.uk
Compressive power spectrum sensing for vibration-based output-only
system identification of structural systems in the presence of noise
Bamrung Tau Siesakul
a
, Kyriaki Gkoktsi
a
, Agathoklis Giaralis*
a
a
Department of Civil Engineering, City University London, Northampton Square, EC1V 0HB,
London, UK
ABSTRACT
Motivated by the need to reduce monetary and energy consumption costs of wireless sensor networks in undertaking
output-only/operational modal analysis of engineering structures, this paper considers a multi-coset analog-to-
information converter for structural system identification from acceleration response signals of white noise excited linear
damped structures sampled at sub-Nyquist rates. The underlying natural frequencies, peak gains in the frequency
domain, and critical damping ratios of the vibrating structures are estimated directly from the sub-Nyquist measurements
and, therefore, the computationally demanding signal reconstruction step is by-passed. This is accomplished by first
employing a power spectrum blind sampling (PSBS) technique for multi-band wide sense stationary stochastic processes
in conjunction with deterministic non-uniform multi-coset sampling patterns derived from solving a weighted least
square optimization problem. Next, modal properties are derived by the standard frequency domain peak picking
algorithm. Special attention is focused on assessing the potential of the adopted PSBS technique, which poses no sparsity
requirements to the sensed signals, to derive accurate estimates of modal structural system properties from noisy sub-
Nyquist measurements. To this aim, sub-Nyquist sampled acceleration response signals corrupted by various levels of
additive white noise pertaining to a benchmark space truss structure with closely spaced natural frequencies are obtained
within an efficient Monte Carlo simulation-based framework. Accurate estimates of natural frequencies and reasonable
estimates of local peak spectral ordinates and critical damping ratios are derived from measurements sampled at about
70% below the Nyquist rate and for SNR as low as 0db demonstrating that the adopted approach enjoys noise immunity.
Keywords: Sub-Nyquist multi-coset sampling, Multi-band stationary stochastic processes, Power spectrum blind
estimation, Output-only system identification, ARMA filter, analog-to-information converter, additive white noise.
1. INTRODUCTION
Operational modal analysis (OMA), also termed as output-only system identification, is historically, the first, and,
arguably, the most widely used technique for global condition assessment, design verification, and health monitoring of
civil engineering structures and structural components
1,2
. OMA relies on linear system identification and modal testing
concepts and techniques
3
to derive structural dynamic/modal properties (i.e., natural frequencies, damping ratios, and
mode shapes) by acquiring and processing low-amplitude response (output) acceleration signals from linearly vibrating
structures due to ambient noise excitation (input). In practice, the input/excitation is not measured but is assumed to be
sufficiently white (i.e., have a flat spectrum over a range of frequencies within which the structural natural frequencies of
interest lie) by allowing for a long enough observation time. Consequently, ignoring any additive noise, the (analog)
acceleration signals considered in OMA are naturally sparse in the frequency domain since their energy is concentrated
around the natural frequencies of the structural system at which (for lightly damped structures) a prominent local peak in
the spectrum is attained. From a theoretical viewpoint, a plethora of system/modal identification algorithms and signal
processing techniques have been considered
2,4,5
to address the two major challenges of OMA: (i) the modal coupling
effect
3
(i.e., the identification of closely spaced natural frequencies and the separation of the corresponding mode shapes)
commonly encountered in higher modes of vibration
1
, and (ii) bias errors due to environmental conditions, measurement
noise, and model uncertainties
4
. From a technological viewpoint, the consideration of wireless sensor networks (WSNs)
has been an important development in OMA in the past 15 years
6,7
. This is because WSNs allow for less obtrusive and
more economical and rapid implementation of OMA compared to arrays of tethered sensors.
To this end, significant research effort has been devoted to develop low energy-consuming wireless sensors while
maintaining their cost and sophistication at a reasonable level
6
. Specifically, in a typical WSN deployment for OMA,
sensors acquire acceleration measurements, perform local data processing at on-board micro-processors, and transmit
data to a base station for further processing. The local data processing step is undertaken to: (i) reduce the amount of
required data to be transmitted within the WSN, and to (ii) minimize the sensor energy consumption (and, therefore, the
requirements for energy harvesting and/or frequency of battery replacement), since data transmission is the most power
consuming operation. This step commonly involves off-line lossy or lossless data compression of stored discrete-time
signals acquired at an analog-to-digital converter (ADC) sampling by, at least, twice the Nyquist rate. The compressed
data are transmitted to a base station where they are de-compressed to reconstruct the originally acquired signals (or an
estimate of them in case of lossy compression). These signals are then further processed by standard OMA algorithms to
derive structural dynamical properties. Despite these efforts, current wireless sensors require battery replacement at
intervals of few weeks to few months, depending on various factors such as the sampling frequency, the duration of each
monitoring interval, the on-board hardware and software to be executed.
In this respect, quite recently, certain compressive sensing (CS) or, in general, sub-Nyquist sampling techniques have
been considered by various researchers (including the authors) in undertaking OMA operations
8-11
, recognizing that such
techniques can potentially minimize monetary and energy consumption costs in WSNs
12
. Specifically, O’Connor et al.
10
reported significantly reduced energy consumption in a long-term field deployment of wireless sensors acquiring
randomly sampled sub-Nyquist measurements compared to conventional (Nyquist sampling) sensors. In this application,
the acquired compressed sensed signals were transmitted to a base station and reconstructed in the time-domain (at
Nyquist rate) by means of a standard CS reconstruction algorithm. The reconstructed signals from each channel/sensors
were finally Fourier-transformed to obtain frequency response functions (FRFs) and the mode shapes were derived using
the standard frequency domain decomposition algorithm of OMA
13
. Following a similar CS-based strategy, that is,
considering reconstructed Nyquist sampled signals in the time-domain from randomly sampled sub-Nyquist
measurements, Yang and Nagarajaiah
11
explored the potential of CS-based OMA in conjunction with blind source
separation for mode shape and natural frequency estimation. In a different study, Park et al.
8
recognized that for the
purposes of modal system identification of linear systems, signal reconstruction in time-domain from CS measurements
is not necessary. In this regard, they considered a singular value decomposition based algorithm to retrieve mode shapes
directly from sub-Nyquist non-uniform random measurements assumed to be acquired by means of a particular CS-based
analog-to-information (AIC) sampling device, namely the random demodulator (RD)
14
. Note that the RD assumes a
multi-tone analog signal model attaining a spectrum with a number of “spikes” (harmonics) and zero elsewhere. In the
theory of structural dynamics, such a signal model is consistent with the (deterministic) response of free oscillating non-
decaying (i.e., undamped) structures. Indeed, the analytical work of Park et al.
8
rely on the above assumption, which is
not in alignment with the typical analog OMA signals (multi-band wide sense stationary random signals/ stochastic
processes). Still, reasonably accurate results in terms of mode shapes were derived from noisy field recorded data
pertaining to a specific bridge structure using the proposed algorithm for system identification.
Notably, the above discussed CS-based approaches for OMA are sensitive to added noise in the CS measurements since
they rely on signal sparsity considerations
15
. It has been shown that a signal reconstruction error of approximately 3 dB
of signal-to-noise ratio (SNR) for each halving of the number of measurements (per octave of subsampling) occurs
which implies signal noise amplification at higher compression rates. In fact, a trade-off between noise level and
sampling rate exists for CS-based AICs (i.e., random sampling assuming analog signal sparsity): within noisy
environments (low SNRs), faster sampling rates and, therefore, an increased number of measurements are required to
extract the underlying signal information
16
.
Inspired by recent developments in power spectrum blind sampling (PSBS) driven mainly by applications in
telecommunications and especially by cognitive radio
17
, the authors
9
considered the use of an AIC device
18
implementing deterministic non-uniform periodic multi-coset sampling
19
to acquire sub-Nyquist acceleration
measurements from linear white-noise excited structural systems. Next, a PSBS approach that does not require any
signal sparsity assumption
18,20
has been applied to approximate the underlying FRFs directly from the sub-Nyquist
measurements. Finally, the standard “peak picking” in the frequency domain OMA algorithm is applied to extract the
useful information from these FRFs. That is, the central frequency of each spiky “occupied” band in the spectrum
corresponding to a natural frequency of the system, the width of the occupied bands associated with structural damping,
and the local peak values of the occupied bands associated with mode shapes which can be retrieved in a multi-channel
setting. It is noted in passing that in cognitive radio applications the goal is to efficiently detect unoccupied bands in a
wide spectral range of telecommunication signals; a relatively simple task compared to the retrieval of modal properties
in OMA. Importantly, the above approach, adopted in this paper as well, does not require signal reconstruction and does
not pose any sparsity requirements to the analog signal. Therefore, it is expected to be insensitive to added noise to the
sub-Nyquist measurements offering an important advantage over the CS-based random sampling considered in the
previously reviewed approaches
8,10,11
. In this regard, this paper furnishes novel numerical results to assess the potential
of the adopted PSBS approach
18,20
to retrieve useful data for OMA from a single channel/sensor acquiring noisy sub-
Nyquist multi-coset samples. To this aim, a benchmark space truss structure with closely-spaced natural frequencies is
considered as a testbed
21
together with a novel frequency domain simulation-based framework which allows for the
efficient generation of wide sense stationary structural response acceleration signals with additive noise consistent with
the assumptions of OMA.
The remainder of the paper is organized as follows. Section 2 introduces the adopted AIC device and multi-coset sub-
Nyquist sampling strategy for stochastic processes (random signals). Section 3 reviews the mathematical details to
accomplish PSBS directly from the coset sub-Nyquist measurements, while section 4 outlines the optimization problem
that needs to be solved to design efficient deterministic multi-coset sampling patterns. Section 5 provides pertinent
numerical results to assess the robustness of the adopted approach to noise and, lastly, section 6 summarizes conclusions
and points to directions for future work.
2. MULTI-COSET SUB-NYQUIST SAMPLING
Multi-coset Sampling Strategy and Device
Let x(t) be a continuous in time t real-valued wide-sense-stationary random signal (or stochastic process) characterized in
the frequency domain by the power spectrum P
x
(ω) band-limited by 2π/T. It is desired to sample x(t) at a rate lower than
the Nyquist sampling rate 1/Τ (in Hz), and still be able to obtain a useful estimate of the power spectrum P
x
(ω). To this
aim, the multi-coset sampling strategy is herein adopted
19
, according to which the grid of Nyquist samples x(nT) is
divided into blocks of N consecutive samples and from each block M (<N) Nyquist-rate samples are selected. The
resulting sampling is periodic with period N; non-uniform since any subset of M samples may be selected from a total of
N Nyquist-rate samples within each block; and deterministic since the position of the M samples on the Nyquist grid of
samples x(nT) is defined a priori and applies to all considered blocks. The above sampling strategy can be implemented
by utilizing M interleaved ADC units operating at a sampling rate 1/(NT) as discussed in Ariananda and Leus
18
. A
discrete-time model of such a sampling device is shown in Figure 1 in which the discrete-time signal x[n]= x(n/T) enters
M branches and at each m branch (m= 0,1,…,M-1), the signal is convolved (filtered) by an N-length sequence c
m
[n] and
down-sampled by N.
Figure 1. Discrete-time model of the considered AIC multi-coset sampling device
18
.
The selection of M samples (sampling pattern) within each block is governed by the coefficients c
m
[n] of the filter
written as
1, ,
[]
0, ,
m
m
m
nn
cn
nn
(1)
where there is no repetition in n
m
, i.e.,
12
12
,.
mm
n n m m
(2)
The output of the m-th branch is given by
[ ] [ ],
mm
y k z kN
(3)
where
[]
m
z
is expressed as
0
1
[ ] [ ] [ ].
mm
kN
z n c k x n k
(4)
Relation of the Input and Output Correlation Functions
Consider the cross-correlation function between the output sequences of the different branches of the device in Figure 1
and the autocorrelation function of the input signal to the device given by
12
12
,y
[ ] E [ ] [ ] ,
mm
y y m m
r k y l y l k
(5)
and
x
[ ] E [ ] [ ] ,
x
r n x m x m n
(6)
respectively, where E
a
{·} is the mathematical expectation operator with respect to a and the “*” superscript denotes
complex conjugation. Further, consider the pattern cross-correlation function between the different sampling patterns of
each branch of the same device expressed as
12
12
0
,
1
[ ] [ ] [ ].
mm
c c m m
kN
r n c k c k n
(7)
Substituting Eq. (1) into Eq. (7) yields
21
12
,
[ ] [ ( )].
mm
c c m m
r n n n n
(8)
where δ[n]=1 for n=0 and δ[n]=0 for n≠0.
It can be shown that the following relationship holds
18
1
0
[ ] [ ] [ ],
y c x
l
k l k l
r R r
(9)
where
y
r
[k] is the M
2
-by-1 vector collecting the output cross-correlation functions
12
,
mm
yy
r
between the M branches of the
considered sampling device evaluated at index k, that is,
0 0 0 1 1 0 1 1
T
, , , ,
[ ] [ ] [ ] [ ] [ ] ,
M M M
y y y y y y y y y
k r k r k r k r k
r
(10)
where the “T” superscript denotes matrix transposition.
x
r
[n] is the N-by-1 vector collecting the input autocorrelation function evaluated at certain indices as in
T
[ ] [ ] [ 1] [( 1) 1] ,
x x x x
n r nN r nN r n N r
(11)
and R
c
[l] is the M
2
-by-N matrix defined as
0 0 0 1 1 0 1 1
T
, , , ,
[ ] [ ] [ ] [ ] [ ] ,
M M M
c c c c c c c c c
l l l l l
R r r r r
(12)
where
