Bradford Thesis Twoside PDF
Bradford Thesis Twoside PDF
Thesis by
2006
(Defended July 19, 2006)
ii
c 2006
Samuel Case Bradford v
All Rights Reserved
iii
Acknowledgements
Many thanks to my advisor, Tom Heaton, for his valuable guidance as I worked
through this research. Always enthusiastic, I am very grateful to him for helping me
expand on aspects of my thesis that I found intriguing, and for applying his insight
to my studies.
Large portions of my thesis are dependent on digitized data, and as such I would
like to acknowledge the many organizations who have assisted with the collection
and digitization of data: NSMP for digitization of Millikan Library film records,
CGS for Imperial County Services Building records, SCEDC for Southern Califor-
nia seismic data, USC for Millikan Library records (I would like to thank Professor
Todorovska and Professor Trifunac for other appreciated discussions), USGS for the
Millikan Library dense instrumentation network, and the PBIC for the loan of instru-
ments/recorders for the Broad Center ambient vibration test.
My thesis has benefitted from the support of the entire Thomas Building – in-
cluding Carolina Oseguera, Julie Wolf, Swami Krishnan, Jeff Scruggs, Andy Guyader,
Judy Mitrani-Reiser, Sam Feakins-Daly, Alex Taflanidis, and my officemate Jing
Yang. In particular, the experimental guidance of Javier Favela, the thoughtful edit-
ing of Anna Olsen, and the many contributions of Matt Muto made this a stronger
work.
Special thanks to Georgia Cua and John Clinton, for the good times and compan-
ionship, and I am grateful as well to John for the technical foundation upon which
much of my research depends.
iv
In a larger sense, the support and understanding of my family has helped me
through the years, and I would like to thank my parents (Clete and Teri and Ken
and Jessica) and grandparents (Amber and Carl and Julia) for their encouragement
to pursue my goals.
Finally, my deepest thanks are reserved for my patient wife, Charlotte.
v
Abstract
Contents
Acknowledgements iii
Abstract v
1 Introduction 1
2 Time-Frequency Representations 7
2.1 Frequency Content of Non-Stationary Signals . . . . . . . . . . . . . 9
2.2 Uncertainty Principle in Temporal and Frequency Resolution . . . . . 13
2.3 Support & Marginal Conditions . . . . . . . . . . . . . . . . . . . . . 17
2.4 Short-Time Fourier Transform & Spectrogram . . . . . . . . . . . . . 18
2.5 Continuous Wavelet Transform & Scalogram . . . . . . . . . . . . . . 20
2.6 Gabor Transform . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.7 Wigner-Ville Distribution . . . . . . . . . . . . . . . . . . . . . . . . 23
2.8 Classes of Time-Frequency Representations . . . . . . . . . . . . . . . 25
2.8.1 Cohen’s Class . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.8.2 Affine Class . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.9 Reduced Interference Distribution . . . . . . . . . . . . . . . . . . . . 28
2.10 Smoothed Pseudo Wigner-Ville Distribution . . . . . . . . . . . . . . 29
2.11 Other Time-Frequency Representations . . . . . . . . . . . . . . . . . 29
2.11.1 Choi-Williams . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.11.2 Born-Jordan . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.11.3 Zhao-Atlas-Marks . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.12 Sample Signal with Transient Components . . . . . . . . . . . . . . . 32
viii
2.12.1 Spectrogram . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.12.2 Continuous Wavelet Transform . . . . . . . . . . . . . . . . . 40
2.12.3 Gabor Transform . . . . . . . . . . . . . . . . . . . . . . . . . 42
2.12.4 Wigner-Ville Distribution . . . . . . . . . . . . . . . . . . . . 42
2.12.5 Smoothed Pseudo Wigner-Ville Distribution . . . . . . . . . . 45
2.12.6 Reduced Interference Distribution . . . . . . . . . . . . . . . . 45
2.12.7 Choi-Williams Distribution . . . . . . . . . . . . . . . . . . . . 52
2.12.8 Born-Jordan Distribution . . . . . . . . . . . . . . . . . . . . 52
2.12.9 Zhao-Atlas-Marks Distribution . . . . . . . . . . . . . . . . . 52
2.12.10 Masked Wigner-Ville Methods . . . . . . . . . . . . . . . . . . 52
2.12.11 Simple Masking Representation . . . . . . . . . . . . . . . . . 61
2.13 Further Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
2.13.1 Chirp Function . . . . . . . . . . . . . . . . . . . . . . . . . . 66
2.13.2 Double Chirp Function . . . . . . . . . . . . . . . . . . . . . . 66
2.13.3 Millikan Library: 2005 Parkfield Earthquake . . . . . . . . . . 73
2.14 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
3 Time-Frequency Representations:
Response of Simulated Nonlinear Systems 85
3.1 Distributed Element Model . . . . . . . . . . . . . . . . . . . . . . . 86
3.2 Nonlinear Finite Element Building . . . . . . . . . . . . . . . . . . . 96
3.2.1 Response to Varying Levels of Input Motion . . . . . . . . . . 102
3.3 Conclusions and Discussion . . . . . . . . . . . . . . . . . . . . . . . 112
4 Time-Frequency Representations:
Response of Instrumented Structures 113
4.1 Millikan Library (Caltech), Pasadena, CA . . . . . . . . . . . . . . . 114
4.1.1 Millikan Library: Sensitivity to Weather Patterns . . . . . . . 116
4.1.2 Millikan Library: 1971 San Fernando Earthquake . . . . . . . 123
4.1.3 Millikan Library: 1987 Whittier Narrows Earthquake . . . . . 138
4.1.4 Millikan Library: 1991 Sierra Madre Earthquake . . . . . . . . 140
ix
4.1.5 Millikan Library: 1994 Northridge Earthquake & Aftershock . 152
4.1.6 Millikan Library: Amplitude & Frequency Nonlinearity . . . . 175
4.2 52-Story Office Building, Los Angeles, CA
(1994 Northridge Earthquake) . . . . . . . . . . . . . . . . . . . . . . 184
4.3 Imperial County Services Building, El Centro, CA
(1979 Imperial Valley Earthquake) . . . . . . . . . . . . . . . . . . . 190
4.4 Conclusions and Discussion . . . . . . . . . . . . . . . . . . . . . . . 193
5 Conclusions 195
Bibliography 251
xi
List of Figures
4.1 Deviation from the mean natural frequency for the 3 fundamental fre-
quencies at Millikan Library station MIK, May 2001 – Nov 2003 . . . . 117
4.3 Deviation from the mean for the natural frequencies of Millikan Library
during February 2003, which includes a major rainstorm . . . . . . . . 120
4.4 Deviation from the mean for the natural frequencies of Millikan Library
in January 2003, as in Figure 4.3, for strong winds . . . . . . . . . . . 121
4.5 Deviation from the mean for the natural frequencies of Millikan Library
in August and September of 2002, as in Figure 4.3, for high temperatures122
B.1 Robert A. Millikan Memorial Library: View from the northeast . . . . 212
B.2 Millikan Library structural diagrams. . . . . . . . . . . . . . . . . . . . 213
B.3 Graphical depiction of Table B.1, the historical behavior of Millikan
Library . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 216
B.4 Kinemetrics VG-1 Synchronized Vibration Generator (Shaker) . . . . . 217
B.5 Lin-Log normalized peak displacement curves for the frequency sweep
performed on July 10, 2002 . . . . . . . . . . . . . . . . . . . . . . . . 220
B.6 Resonance curves and mode shapes for the E-W fundamental mode un-
der two loading conditions . . . . . . . . . . . . . . . . . . . . . . . . . 224
B.7 Resonance curves and mode shapes for the N-S fundamental mode under
two loading conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 226
xix
B.8 Resonance curves and mode shapes for the Torsional fundamental mode 227
B.9 Second and third E-W modes (first and second E-W overtones) . . . . 229
B.10 Resonance curves and mode shapes for the second NS mode (first NS
overtone) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 230
B.11 Resonance curves and mode shapes for the second Torsional mode (first
Torsional overtone) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 232
B.12 Least squares curve fitting for E-W and N-S modes. Linear tilt and
translation removed when calculating the best fit. . . . . . . . . . . . . 234
B.13 From left to right, theoretical mode shapes for the fundamental mode
(1st mode) and the first two overtones (2nd and 3rd modes) for a can-
tilevered bending beam . . . . . . . . . . . . . . . . . . . . . . . . . . 236
B.14 From left to right, theoretical mode shapes for the fundamental mode
(1st mode) and the first two overtones (2nd and 3rd modes) for a can-
tilevered shear beam . . . . . . . . . . . . . . . . . . . . . . . . . . . . 237
xx
xxi
List of Tables
Chapter 1
Introduction
Modern digital seismology has advanced considerably in recent years. As seismic in-
strumentation systems have become more sensitive and accurate, stations have been
deployed in civil engineering structures such as buildings, bridges, and dams, allow-
ing for an unprecedented investigation into their dynamic properties. The fields of
civil/structural engineering and earthquake engineering have benefitted greatly from
this increased instrumentation of structures. 24-bit digitization systems, now com-
monplace, allow for the analysis of strong earthquake motions and ambient vibrations
on the same instrument, with a signal-to-noise ratio of over 100 for typical ambient
vibration levels. When these instruments and digitization systems are combined with
a continuous telemetry and storage system, such as that found in the Southern Cali-
fornia Seismic Network, building behavior can be observed over a timespan of seconds,
hours, days, weeks, months, or even years; this provides a rich data set for investiga-
tions into changes in the dynamic properties of these buildings.
Transform methods are fundamental in understanding signals from a wide variety
of fields. Scientists and engineers are often restricted to studying systems through
signals associated with these systems (vs. physically inspecting the components of the
system), and the properties of the signal are typically difficult to identify in the time-
domain. Transforming the signal into the frequency domain via the Fourier Transform
identifies the frequency content of the signal, which is often more useful than time-
domain information for investigation of dynamic properties. A Fourier Transform of
a time series (such as an earthquake record from an instrumented structure) contains
2
information regarding the frequency content of the signal, but it cannot resolve the
exact onset of changes in natural frequency, as temporal information is contained only
in the phase of the transform. Comparing a Fourier Transform of the first portion of
the record with a later portion reveals information about the evolution of the dynamic
properties, but decreases the frequency resolution of each portion of the record – this
idea is at the heart of the spectrogram. By applying the Fourier Transform with a
running window, a properly constructed spectrogram is better able to resolve temporal
evolution of frequency content. However, there is a trade-off in time resolution versus
frequency resolution in accordance with the uncertainty principle.
I present a mathematical foundation for analyzing the evolution of frequency con-
tent in a signal and apply these techniques to synthetic records from linear and non-
linear finite element models. These analysis techniques are then applied to records
from instrumented structures, with examples drawn from forced vibration tests, am-
bient vibrations, and buildings damaged during historical earthquakes. Many of these
techniques are adaptations of the Wigner-Ville Distribution, which is one of the old-
est and most fundamental time-frequency representations. As investigations into in-
strumented structures are dependent on digital signal processing, applying advanced
signal processing tools to this data allows for a deeper understanding of the system’s
evolution with time.
2
34
LC M5.3 !=57
WN M6.1 !=19
SM M5.8 !=18
BB M5.4 !=119
NR M6.7 !=34
BH M4.2 !=26
Frequency, Hz
SS M6.5!=323
351
1.4
Figure 1.1: Historical behavior of Millikan Library - EW (red) and NS (blue) fun-
damental frequencies since construction. Crosses indicate frequencies from forced
vibration tests. Circles indicate estimates from recorded earthquakes, and the asso-
ciated numbers in italics are peak acceleration values recorded for the event (cm/s2 ).
Dashed lines represent the observed natural frequencies of the library, and the shaded
region is the likely variance from such factors as weather conditions, weight config-
uration of the shaker used for forced vibration tests, and experimental error. Note
the significant drop in observed frequencies during the events and the permanent
decrease in both the EW and NS natural frequencies. [Earthquake Abbreviations:
LC: Lytle Creek, SF: San Fernando, WN: Whittier Narrows, SM: Santa Monica, NR:
Northridge, BH: Beverly Hills, BB: Big Bear, SS: San Simeon] (Adapted from Clinton
et al. (2006).)
5
the effects of these factors generally requires that the operator have some knowledge
about the system and its properties.
Time-frequency analysis, exploring system behavior through the exploration of the
behavior in the time-frequency plane, provides important insight into the dynamic
properties of a system. Wandering natural frequencies can be extracted from the
ridges of a 2D representation, and the onset of changes in the dynamic properties
can be further correlated with physical properties of the system such as deformation,
damage, peak velocity, etc. While these methods can be used as a starting point for
parameterization in the general system identification sense, they also by themselves
provide interesting information regarding the character and extent of changes in the
apparent natural frequency of systems. Changes in physical properties are typically
represented in the time-frequency plane as changes in the observed natural frequency
of some signal of interest.
6
7
Chapter 2
Time-Frequency Representations
1
The Hilbert Transform is sometimes referred to as a “quadrature filter” and the transformed
signal as the “quadrature signal.” Running the transform four times will return the signal to the
original signal, as each transformation shifts real frequencies by π/2 and negative frequencies by
−π/2.
9
distribution of pseudo-energy in the time-frequency plane is a useful approximation
to the distribution of energy at each instant, and time-frequency representations can
be interpreted as approximate energy even in cases where the representation is not
entirely positive.
The Fourier Transform (FT) of a signal decomposes the original signal into harmonic
components, identifying the spectral content of the signal. This process allows for
system identification in terms of the natural frequencies and corresponding mode
shapes, which are directly related to the physical properties of the system. However,
when analyzing signals with evolving frequency content, the amplitude of the FT
does not give the full information regarding the behavior of the system (Bracewell,
2000). A sample signal with evolving frequency is shown in Figure 2.1(a), along with
the amplitude of the FT of the signal. While the FT correctly identifies the main
components of the signal, it allows for neither a straightforward identification of the
onset of each signal component nor an identification of which frequency component
arrived first. All such temporal information is contained in the phase of the transform,
but cannot be easily extracted for system identification purposes.
In Figure 2.2, I introduce a Time-Frequency Representation (TFR) of the synthetic
signal. This representation creates a frequency estimate at each instant in the signal
and a more thorough understanding of the evolution of the natural frequency. If the
signal represents the behavior of a system, then changing frequency content can be
correlated with the evolving properties of a system (e.g., changing stiffness, damage
to a structural component).
To create the time-frequency representations in Figure 2.2, a sliding window is
used to split the signal into segments, and the Fourier Transform of each segment is
then assembled into the final time-frequency matrix. Two windowing widths are pre-
sented to emphasize the trade-off in temporal resolution versus frequency resolution.
In the top figure, a wide window is used, and the coarse temporal resolution of the
10
0.5
Amplitude
0.025
0
−0.5
Frequency (Hz)
0.02
−1
128 256 384 512 640 768 896 1024 1152 1280 1408 1536 1664 1792 1920 2048
Time
Fourier Transform of Test Signal
400 0.015
Fourier Amplitude
300
200 0.01
100
0 0.005
0.005 0.01 0.015 0.02 0.025 0.03 200 400 600 800 1000 1200 1400 1600 1800 2000
Frequency (Hz) Time
(a) Sample signal with evolving frequency, and (b) Spectrogram. this distribution has limited
Amplitude of the Fourier Transform – note resolution and takes on a “blocky” structure
that the Fourier Transform, while correctly due to the trade-off between temporal reso-
identifying the components of the signal, lution and frequency resolution in the time-
does not allow for a straightforward inter- frequency plane. However, this representa-
pretation of frequency evolution. tion clearly captures information about pro-
gression of the mean frequency during the sig-
nal.
Scalogram −− Continuous Morlet Wavelet Wigner−Ville Distribution
0.03 0.03
Approximate Frequency [scale = 5/(2*π*frequency)]
0.025 0.025
Frequency (Hz)
0.02 0.02
0.015 0.015
0.01 0.01
0.005 0.005
128 256 384 512 640 768 896 1024 1152 1280 1408 1536 1664 1792 1920 2048 128 256 384 512 640 768 896 1024 1152 1280 1408 1536 1664 1792 1920 2048
Time Time
(c) Scalogram of sample signal (Continuous Mor- (d) Wigner-Ville Distribution for the sample sig-
let Wavelet). Wavelet transformations create nal. Note that resolution of the signal is quite
a time-scale representation, where scale has crisp, though the strong interference terms
an approximate equivalency to frequency. make this method unsuitable for general sig-
nal analysis.
Reduced Interference Distribution Smoothed Pseudo Wigner−Ville Distribution
0.03 0.03
0.025 0.025
Frequency (Hz)
Frequency (Hz)
0.02 0.02
0.015 0.015
0.01 0.01
0.005 0.005
128 256 384 512 640 768 896 1024 1152 1280 1408 1536 1664 1792 1920 2048 128 256 384 512 640 768 896 1024 1152 1280 1408 1536 1664 1792 1920 2048
Time Time
(e) Reduced Interference Distribution for the (f) Smoothed Pseudo Wigner-Ville Distribution
sample signal. for the sample signal
Figure 2.1: Selected Time-Frequency Representations: (a) Sample signal, (b) Spec-
trogram (Section 2.4), (c) Scalogram (Section 2.5), (d) Wigner-Ville Distribution
(Section 2.7), (e) Reduced Interference Distribution (Section 2.9), and (f) Smoothed
Pseudo Wigner-Ville Distribution (Section 2.10).
11
representation is immediately obvious . In the second figure, a narrower window is
used which increases the temporal resolution, and it is now easier to identify the evo-
lution of the frequency content contained in the signal; however, the narrow window
decreases the maximum frequency resolution for each slice. Both of these windowing
choices smear information in the time-frequency plane along the time and frequency
axes. The third spectrogram plot in Figure 2.2 is an instantaneous frequency repre-
sentation, an improvement over the methods in the first two plots. This method still
leaks energy along both axes, but closely matches the theoretical frequency content,
shown in the fourth plot. This instantaneous spectrogram, also presented in Fig-
ure 2.1(b), still has a noticeable block-like structure, a result of the inability of this
method to create a true instantaneous energy estimation. The uncertainty principle
(Section 2.2) precludes a TFR with perfect resolution in both time and frequency,
but various methods have been proposed that can create useful representations for
instantaneous frequency content.
In time-frequency representations, the goal is to create a distribution which cor-
rectly identifies energy in the time-frequency plane. For the spectrogram example,
the information in the time-frequency plane closely matches that of the theoretical
components of the signal, shown in the last plot of Figure 2.2. However, the spec-
trogram method introduces some further complications for analyzing the temporal
evolution of the frequency content of the signal. A long time window will smear
the time-frequency information across the time axis, changing the perceived duration
and onset of a signal component. A shorter time window, while improving temporal
resolution, decreases the maximum resolution along the frequency axis. Temporal
resolution and frequency resolution are inversely proportional in accordance with the
uncertainty principle, which limits the effective resolution of all time-frequency rep-
resentations.
Another method for analyzing the instantaneous spectrum of a signal is to de-
compose it into different bases: one such choice is the wavelet transformation. The
Continuous Wavelet Transform (CWT) creates a scalogram – a time-scale represen-
tation of the signal, where the scale of the wavelet bases has an intuitive association
12
0.025
0.02
0.015
0.01
0.005
0 128 256 384 512 640 768 896 1024 1152 1280 1408 1536 1664 1792 1920 2048
Spectrogram, Short Time Window
Frequency (Hz)
0.025
0.02
0.015
0.01
0.005
0 128 256 384 512 640 768 896 1024 1152 1280 1408 1536 1664 1792 1920 2048
Spectrogram, instantaneous frequency estimation
Frequency (Hz)
0.025
0.02
0.015
0.01
0.005
128 256 384 512 640 768 896 1024 1152 1280 1408 1536 1664 1792 1920 2048
Idealized frequency of test signal
Frequency (Hz)
0.025
0.02
0.015
0.01
0.005
128 256 384 512 640 768 896 1024 1152 1280 1408 1536 1664 1792 1920 2048
Time
The inverse relationship between temporal resolution and frequency resolution has
been formalized as the Heisenberg-Gabor Uncertainty Principle.
The uncertainty principle as applied to quantum physics (“Heisenberg Uncer-
tainty Principle” or “Heisenberg Indeterminacy Principle”) refers to the relationship
between certain canonically related quantum observables, such as position and mo-
mentum.2 Briefly, if you measure the momentum of a particle p to within an accuracy
of ∆p, then you cannot resolve the position x more accurately than ∆x = h/∆p, where
h is Plank’s constant, ∼ 6.63x10−34 Joule · seconds (Feynman et al., 1966), such that:
∆x · ∆p ≥ h. (2.1)
An analogous inequality can be formulated for time and frequency resolution. The
standard proof offered in many textbooks (Flandrin (1999), Qian (2002), and others)
2
This is a separate principle from the observer effect, which states that observation of a system
changes the phenomenon being observed.
14
is straightforward, depending on the definition of the Fourier Transform, Parseval’s
Formula (Equation 2.2), and the Cauchy-Schwarz Inequality (Equation 2.6).
For vectors a(t) and b(t) (with Fourier Transforms A(ω) and B(ω)), the definition
of the convolution operator leads to Parseval’s Formula:
Z+∞ Z+∞
∗ 1
a(t)b (t)dt = A(ω)B ∗ (ω)dω. (2.2)
2π
−∞ −∞
Z+∞ Z+∞
2 1
|b(t)| dt = |B(ω)|2 dω. (2.3)
2π
−∞ −∞
Parseval’s Formula implies a conservation of energy between the time domain and
frequency domain (the form of Equation 2.3 is also referred to as the Plancherel
Formula (Mallat, 1999)).
d
For the case of a derivative, dt
b(t) ↔ iωB(ω):
Z+∞ Z+∞
db 2
1
(t) dt =
dt ω 2 |B(ω)|2 dω. (2.4)
2π
−∞ −∞
I will also require the Cauchy-Schwarz Inequality, which for functions a(t) and b(t)
can be expressed as:
where h · , · i denotes an inner product. For the usual integral (L2 ) inner product for
complex-valued functions, Cauchy-Schwarz implies that:
Z 2 Z Z
a(t)b(t) dt ≤ |a(t)|2 dt · |b(t)|2 dt.
(2.6)
15
With these formulas, I can now derive the uncertainty in temporal and frequency
resolution for a signal, x(t), (Fourier Transform X(ω)). I begin with the assumption
of finite energy:
Z∞
Ex = |x(t)|2 dt < +∞. (2.7)
−∞
Without loss of generality I can further define the axes such that:
Z∞
t|x(t)|2 dt = 0, (2.8)
−∞
Z∞
ω|X(ω)|2 dω = 0. (2.9)
−∞
Taking the second moment as a measure of time and frequency support (temporal
duration and frequency bandwidth):
Z∞
2
∆t = t2 |x(t)|2 dt, (2.10)
−∞
Z∞
1 2
∆ω 2 = ω |X(ω)|2 dω, (2.11)
2π
−∞
Z∞
dx 2
= (t) dt.
dt (Equation 2.4.) (2.12)
−∞
Z∞ Z∞
dx 2
∆t2 ∆ω 2 = 2 2
t |x(t)| dt · dt (t) dt.
(2.13)
−∞ −∞
16
Applying the Cauchy-Schwarz inequality of Equation 2.6 to the right-hand side of
Equation 2.13 gives:
∞ 2
Z∞ Z∞ 2 Z
2 2
dx ∗ dx
t |x(t)| dt · dt (t) dt ≥ t x (t) dt (t)dt
(2.14)
−∞ −∞ −∞
Take |I|2 to be the right hand side of Equation 2.14 such that:
Z∞
dx
I = t x∗ (t) (t)dt, (2.15)
dt
−∞
Z∞
dx∗
∗ +∞ ∗
I = [t x (t)x(t)]−∞ − t x(t) + x (t)x(t) dt, (2.16)
dt
−∞
Z∞ Z∞
+∞ dx∗
t |x(t)|2 −∞ − |x(t)|2 dt −
I = t x(t)
(t)dt, (2.17)
dt
−∞ −∞
=⇒ I + I ∗ = −1, (2.19)
The restriction in Equation 2.18 is necessary for this form of the proof, though the
final result is valid for all x ∈ L2 (R) (Mallat, 1999). From Equation 2.13, 2.14 ,and
2.22:
1
∆t2 ∆ω 2 ≥ |I|2 ≥ (Real{I})2 = . (2.23)
4
Finally, the time-frequency analog to the Heisenberg Principle, the Heisenberg-Gabor
Uncertainty Principle:
1
∆t · ∆ω ≥ . (2.24)
2
where X(ω) is the Fourier Transform of x(t), and Ex is the total energy of signal x(t).
Satisfying the marginal conditions states that the integral of the TFR along the time
axis is equal to the power spectrum of the signal (the squared Fourier Transform),
while the integral along the frequency axis gives the squared envelope of the original
time series – the energy condition states that the total area of the TFR, the double
integral across time and frequency, is the energy contained in the original signal.
These conditions have an intuitive appeal, as they imply a limitation on the extent
of the signal in the time-frequency plane. In an ideal representation, a signal of short
duration would have a narrow representation along the time axis, and the frequency
content would be localized to the frequency of the signal. The total energy condition
is also a useful property, particularly for the situation where the TFR is being used
18
as an instantaneous energy estimate – bounded energy makes for a more meaningful
physical interpretation.
The marginal conditions are important criteria for constructing a TFR, though
these conditions are neither necessary nor sufficient for the construction of useful rep-
resentations.3 For a graphical example of the time and frequency marginal conditions,
please see Section 2.13.3, and Figures 2.43 and 2.44.
Z∞
X(ω) = x(t)e−iωt dt. (2.28)
−∞
For a signal with evolving frequency content, the Fourier Transform does not give
information regarding the temporal evolution of changes in character of the signal.
However, the Fourier Transform of the first portion of the record can be easily com-
pared with the Fourier Transform of the last portion of the record, revealing the
frequency content in each portion. A logical extension of this concept is to window
the original signal and compute a running Fourier Transform, or Short-Time Fourier
Transform (STFT):
Z∞
ST F T (t, ω) = x(u)h∗ (u − t)e−iωu du, (2.29)
−∞
3
For example, the spectrogram and scalogram (Sections 2.4 and 2.5) are useful representations
that do not satisfy these marginal conditions. The WVD (Section 2.7) satisfies the marginal con-
ditions, but it contains interference terms which make it unsuitable as an instantaneous energy
estimate.
19
of the STFT: 2
Z∞
M
X −1
X̂[n] = x̂[m]e−2πmn/N n = 0, 1, 2, . . . , N − 1, (2.31)
m=0
where x̂[m] is the discrete signal (Â will represent the discrete version of A) and X̂[n]
is the corresponding Discrete Fourier Transform. Equation 2.29 is then replaced with
the Discrete Short-Time Fourier Transform:
L−1
X
ST
\ F T [m, n] = x̂[m]ĥ∗ [k − m]e−2πkn/N
k=0 (2.32)
= ST F T (t, ω)|t=m∆t, ω=2πn/(L∆t) ,
When discussing these discrete representations, some basic relations govern the
amount of independent values present in a signal and its transform. For a discrete
vector x[n] with length N (even), the Discrete Fourier Transform results in a complex
vector of length N/2 (properly, the Fourier Transform is of length N , with a redundant
mirroring about ω = 0 which reduces the effective length of the vector to N/2). This
complex vector has N/2 amplitude values and N/2 phase values for a total of N pieces
of unique information – it is a linear transform with no redundancy. When considering
the amplitude of the Fourier Transform, as in the STFT and spectrogram, there is an
apparent loss of 50% of the information in the signal. The phase of the complex value
is typically discarded, leaving an amplitude vector of length N/2. A non-redundant
20
discrete STFT will be restricted to J-by-K matrices where J × K = N/2. A vector of
length N = 4096, for example, could be composed of an 8-by-256 matrix, a 32-by-64
matrix, a 128-by-16 matrix, etc., all of total size 2048 (N/2). The 8-by-256 matrix
has 8 segments along the time domain, each with a frequency resolution of 256 points
across the frequency axis; all information in the first segment is smeared across the
first 1/8th of the TFR. The 128-by-16 matrix would have 128 segments along the
time domain, representing much finer temporal information, but with a resolution
along the frequency axis of only 16 points across the entire frequency range.
This limitation on resolution in time and frequency is another way to interpret the
fundamental restrictions of the uncertainty principle, as an increase in resolution along
the time axis is only possible by decreasing the resolution along the frequency axis. In
practice, spectrograms are created with redundant information, and time-windowing
parameters are selected to localize the expression of energy in the time-frequency
plane. As the window impinges upon a portion of the record, the frequency content
of the portion is reflected in the time-frequency plane, which will smear information
forwards and backwards along the time axis. Taking a shorter window localizes
information in time at the cost of a wider representation along the frequency axis.
where a is the scale (dilation factor) for the analyzing wavelet, Ψ. To create the
wavelet transformation, the original wavelet is smoothly shifted over the full signal
(dilating and contracting as necessary), creating a scale estimate for each instant of
the original signal. Different wavelet families have been proposed (e.g., Harr, Mexican
Hat, Meyer, Daubechies) with various mathematical properties (Mallat, 1999). For
21
the continuation of this thesis, I will make use of the Morlet Wavelet for wavelet
transform methods, defined as a Gaussian tapered cosine.
t2
Ψmorlet (t) = e− 2 cos(5t). (2.34)
Wavelets and the Continuous Wavelet Transform are well suited to detect changes
in a signal that occur over a short timespan, and these methods have been used
with great success for detecting anomalies in signals. Using wavelet techniques, it
is possible to identify, with great precision, the exact temporal onset of “spikes” or
“knocks” in a signal, as well as correct for long-period drifting of a signal (Boashash,
2003).
Figure 2.3: Morlet wavelet and equivalent frequency content at different “scale”
values. The Morlet wavelet is one family of wavelets based on a Gaussian-tapered
cosine. The scale of the wavelet represents the dilation/expansion, which has an
approximate inverse relationship to the frequency content of the wavelet (frequency
5/2π
≈ scale ). By definition, wavelets are not monochromatic, and wavelet transformations
create a time-scale representation (hence “scalogram”) rather than a time-frequency
representation. For higher frequencies (lower scale values), the bandwidth of the
frequency content grows proportionally to the center frequency of the wavelet. This
creates a diffuse frequency approximation at high frequencies.
23
2.6 Gabor Transform
∞
X
GABORx (λ, µ) = x(n)G∗ (n − λL)e−i(µ/W )(n−λL) , (2.36)
n=−∞
µ
= ST F Tx (λL, ), [as in Equation 2.32]
W
Recent signal processing advances have led to a new class of TFR tools. These new
methods allow for distributions that function more accurately than the spectrogram
method, though the uncertainty principle of Section 2.2 does preclude a perfect in-
stantaneous frequency estimation at each point in the time-frequency plane.
24
The most fundamental of these TFR methods is the Wigner-Ville Distribution
(WVD). Many other TFR methods (including the spectrogram and the scalogram)
can be derived from the WVD, with a suitable choice of smoothing factors. For a
signal, s(t), with analytic associate x(t), the Wigner-Ville Distribution, W V Dx (t, ω),
is defined as: Z +∞
τ τ
W V Dx (t, ω) = x(t + )x∗ (t − )e−iωτ dτ. (2.37)
−∞ 2 2
E. Wigner, who first proposed this transformation (Wigner, 1932), was applying
this equation to derive a joint representation of position and momentum in quantum
mechanics – position and momentum are the standard examples of a canonically
related pair. The Heisenberg Uncertainty Principle (Section 2.2) maintains that at
the quantum level the accuracy of a measurement in position is inversely related to
the accuracy of the measurement in momentum. In the field of signal processing, the
inverse relationship between temporal and frequency resolution is analogous to the
inverse relationship between quantum measurements of momentum and position.
The WVD is an entirely real-valued function which satisfies the marginal condi-
tions (Section 2.3), and for certain classes of signals it has ideal resolution in time
and frequency. Interference terms and negative values, however, are significant obsta-
cles to using the WVD as a system identification tool for general signals (Bracewell,
2000). The function x(t) appears twice in the integral, as x and x∗ , which makes the
distribution “bilinear” or “quadratic.” Cross-term interference from the quadratic
x(t + τ2 )x∗ (t − τ2 ) term in the WVD generates highly oscillatory interference when
there are multiple frequency components in a signal. This interference can have an
amplitude larger than the expression of the auto terms, and typically ranges from
25
extreme negative to extreme positive values. Negative values in the time-frequency
plane, under an instantaneous energy interpretation of the WVD, would correspond to
negative energy; this negates the physical interpretation of the WVD as representing
the energy content of a signal. However, the satisfaction of the marginal conditions
and many other desirable properties have led to an interest in improving the accuracy
of the WVD as a time-frequency representation, primarily by reducing the effects of
the interference terms on the final result.
Note that the WVD is similar to the Fourier Transform, though instead of trans-
forming the original signal, the kernel of the WVD contains a type of autocorrela-
tion term (in this case, the phase lag of the ambiguity function, or “properly sym-
metrized covariance function” (Flandrin and Martin, 1997)). The Ambiguity Func-
tion, AFx (τ, ν), is the 2-D Fourier Transform of the WVD and has the expression:
Z +∞
τ τ
AFx (τ, ν) = x(t + )x∗ (t − )e−iντ dt, (2.38)
−∞ 2 2
which is the same form as the WVD in Equation 2.37, but with integration with
respect to t instead of τ . This transformation is not always real valued and is often
used as a cross-AF to determine correlation functions between different states.
(Sτ X)(ω) = eiωτ X(ω)
Time Shift Covariance : (2.39)
⇒ TSτ X (t, ω) = TX (t − τ, ω),
(Mν X)(ω) = X(ω − ν)
Frequency Shift Covariance : (2.40)
⇒ TMν X (t, ω) = TX (t, ω − ν),
(Ca X)(ω) = √1 X(ω/a)
|a|
Scale Shift Covariance : (2.41)
⇒ TCa X (t, ω) = TX (at, ω/a).
The spectrogram and Wigner-Ville Distribution are members of the general class of
time-frequency representations known as Cohen’s Class. It the most general form, a
time-frequency representation in Cohen’s Class, T F RC (t, ω), can be represented by:
+∞
R +∞
R
T F RC (t, ω) = W V Dx (u, ν)fC (t − u, ω − ν) du dν,
−∞ −∞
+∞
R +∞ (2.43)
−iωu iνt
R
fC (t, ω) = ΦC (u, ν)e e du dν,
−∞ −∞
where fC (t, ω) and ΦC (t, ω) are kernel functions (either fC or ΦC uniquely defines the
representation). For fC = δ(t)δ(ω) ↔ ΦC = 1, Equation 2.43 reduces to the WVD.
For fC = W V Dx (−t, −ω) ↔ ΦC = AFx (−t, −ω) (the WVD or Ambiguity Function
(Equation 2.38) for the time/frequency reversed signal), Equation 2.43 reduces to
the spectrogram (c.f. Equation 2.30). This smoothed WVD formulation, for different
kernel properties, makes it straightforward to design time-frequency representations
satisfying useful properties.
Cohen’s class is equivalently defined as representations satisfying the time-shift
and frequency-shift covariances of Equations 2.39 and 2.40, since representations that
satisfy these properties can always be expressed in the form of Equation 2.42 or 2.43.
The scalogram and Wigner-Ville Distribution are members of the Affine Class, distri-
butions that satisfy time-shift and scale-shift covariance (Equations 2.39 and 2.41).
As in the formulation of Cohen’s Class, any time-frequency representation in the
Affine Class, T F RA (t, ω), can be expressed as:
+∞
R +∞
R
T F RA (t, ω) = W V Dx (u, ν)fA (ω(t − u), ν/ω) du dν,
−∞ −∞
+∞
R +∞ (2.44)
−iβλ iαγ
R
fA (γ, λ) = ΦA (α, β)e e dα dβ.
−∞ −∞
For fA (γ, λ) = δ(γ)δ(λ + 1), T F RA is the WVD of the signal; for fA (γ, λ) =
W V Dφ (−γ, −λ), the WVD of a time/frequency reversed wavelet, T F RA reduces
to the scalogram (c.f. Equation 2.35).
28
2.9 Reduced Interference Distribution
The Reduced Interference Distribution (RID) has some advantages over the standard
Wigner-Ville Distribution. The RID and WVD are both time-frequency distribu-
tions in Cohen’s general class, but reduced interference methods are better suited to
transient, non-stationary signals, as well as reducing the quadratic interference that
complicates interpretation of the Wigner-Ville Distribution. In general, a “Reduced
Interference Distribution” refers to any distribution that reduces the expression of the
cross-terms relative to the auto-terms in a quadratic TFR representation; the RID
given here is one particular member of a broad class of such representations.
The Reduced Interference Distribution, RID(t, ω), with kernel Rx (t, τ ), based on
a time series s(t) with analytic associate x(t), is defined as:
R∞
RID(t, ω) = h(τ )Rx (t, τ )e−iωτ dτ,
−∞
+
|τ |
2
(2.45)
g(ν)
1 + cos 2πν x t + ν + τ2 x∗ t + ν − τ2 dν.
R
Rx (t, τ ) = |τ | τ
|τ |
− 2
Note the similarity to the WVD (Equation 2.37) in general form, with the ad-
dition of a smoothing kernel in Rx . In this formulation, h(τ ) is a time-smoothing
window, and g(ν) is a frequency-smoothing window. The smoothing windows can
be adjusted to fit specific needs, see Appendix C for selected software packages with
RID implementations. I have presented a RID formulation for a Hanning window of
[(1 + cos(2πν/τ ))/2]. Other common smoothing kernels include binomial, bessel, and
triangular (or “bartlett”) windows (http://tftb.nongnu.org/, 2006; Boashash, 2003).
For different choices of kernel, the formulation of Rx (t, τ ) will change, but it will
retain its general auto-correlation structure.
29
2.10 Smoothed Pseudo Wigner-Ville Distribution
Z∞ Z∞
SP W V Dx (t, ω) = h(τ ) g(u − t)x(u + τ /2)x∗ (u − τ /2) du e−iωτ dτ. (2.46)
−∞ −∞
Z∞ Z∞ r
α −ν 2 α/(16τ 2 ) τ ∗ τ −iωt
CWx (t, ω) = 2 e x t+ν+ x t+ν− e dν dτ.
16πτ 2 2 2
−∞ −∞
(2.47)
The Choi-Williams Distribution preserves interference horizontally and vertically,
causing a ripple effect when there are auto-terms with the same time or frequency
value. As α goes to zero, the Choi-Williams Distribution reduces to the Wigner-Ville
Distribution. Increasing values of α reduces the expression of cross-terms, but also
begins to interfere with the expression of the auto-terms. As in most time-frequency
representations, there is a trade-off in implementing this distribution.
As a member of Cohen’s Class, the Choi-Williams Distribution can be expressed
πωt 2
as in Equation 2.42 with ΦC = e−[ α ] , a Gaussian function of (tω).
2.11.2 Born-Jordan
The Born-Jordan Distribution (Equation 2.49) preserves time and frequency support,
sin(πωt)
with the kernel ΦC (t, ω) = πωt
, a sinc function of (tω):
Z∞ Z |/2
t+|τ
1
BORN JORDANx = x(u + τ /2)x∗ (u − τ /2) du e−iωt dτ. (2.48)
|τ |
−∞ t−|τ |/2
2.11.3 Zhao-Atlas-Marks
The Zhao-Atlas-Marks (ZAM) Distribution, Equation 2.49, also known as the “Cone-
Shaped Distribution” from the shape of the region |τ | ≥ 2|t| in the t − τ plane, again
reduces interference between arbitrary components in the time-frequency plane. It
is able to almost completely remove cross-term interference from auto-terms with
the same frequency value (“horizontal” interference in the time-frequency plane) but
31
preserves interference from auto-terms with the same time value (“vertical” interfer-
ence). This is a smoothed version of the Born-Jordan Distribution, filtered along the
frequency axis:
Z2t Z∞
−iωt τ ∗ τ
ZAMx (t, ω) = e x u+ x u− f (t − u, τ )du dτ, (2.49)
2 2
−2t −∞
1 −ατ 2
τ
e , |τ | ≥ 2|t|
f (t, τ ) =
0 , otherwise
sin(πωt) −αt2
The Cohen kernel for this distribution is ΦC (t, ω) = πωt
e , which suppresses the
−αt2
Born-Jordan sinc function with the e term. As α in f (t, τ ) increases, cross-terms
are suppressed at the expense of interference with the character of auto-terms (again,
only “vertical” interference, from auto-terms with different frequency components
occurring at the same time instant). The ZAM Distribution has been extensively
used in speech analysis (Qian, 2002) as an alternative to spectrogram methods.
32
2.12 Sample Signal with Transient Components
It is useful to apply these time-frequency analysis methods to a synthetic signal with
known frequency content. The test signal in this section (Figure 2.4, along with the
Fourier Transform) is composed of 5 main segments of 512 points each, for a total
length of N = 5 ∗ 512 = 2560:
3. Baseline signal,
I present linear and logarithmic plots of various TFR methods for this sample sig-
nal. First I present traditional time-frequency representations (spectrogram, wavelets,
Gabor transform) that may be more familiar to the engineering community, then
modern quadratic TFRs and masked TFR methods. The kernel-smoothing done by
modern TFR methods does a better job of removing the cross term interference than
masking/overlay filtering methods, though these techniques do still have problems
with positivity and time-frequency resolution. Each of these methods is subject to
the fundamental restrictions of the uncertainty principle (Section 2.2), which limits
the ultimate resolution in the time-frequency plane.
The Reduced Interference Distribution presented in Section 2.12.6 has certain
properties that make it nearly ideal for the types of signals often encountered in
seismology and civil engineering.
2.12.1 Spectrogram
The Spectrogram, with varying window lengths, shows the evolution of the frequency
content of the signal. The two transient pulses and the baseline signal can be identified
33
Synthetic Signal
(Transient 1) (Transient 2)
512 1024 1536 2048 2560
Time [Points]
Fourier Transform
Amplitude of Fourier Transform
Baseline,
ω = 1/(18π)
Transient 2, Transient 1,
ω = 1/(5π) ω = 2/(5π)
Figure 2.4: Test signal with transient components, five segments of different charac-
ter – 1: Baseline signal (noisy sine wave), 2: Baseline with a high-frequency transient,
which dies out by the end of the segment, 3: Baseline signal, 4: Baseline signal dou-
bles in amplitude, with a intermediate frequency transient, which dies out by the
end of the segment, 5: Baseline signal (doubled). The Fourier Transform correctly
identifies the three main components of the record but does not provide information
regarding the evolution of the signal’s frequency content.
34
in the time-frequency plane, subject to the uncertainty relationship of Section 2.2.
In these figures, the log plot alongside the linear plots allow for a more detailed
investigation of the location of the energy in the Time-Frequency plane.
Of note in Figures 2.5 – 2.9 is the smearing of the transient signal components
along the time axis. Energy in the time-frequency plane appears before the true tem-
poral onset of the transient signals. The spectrogram, as a running Fourier Transform,
extracts frequency content from the signal for a sliding window. However, the win-
dowing involved in the Fourier Transform smears information such that the transient
components appear in the time-frequency plane before their onset in the time se-
ries, as soon as the analysis window impinges on the transient. This is an artifact
of the window used for calculating the Fourier Transform and represents one of the
severe shortcomings in the spectrogram method for system identification. Identify-
ing the onset of changes in the natural frequencies of a structure is necessary for
the identification of changes in stiffness – and ultimately for correlating changes in
dynamic properties with changes in physical properties of the system. The spectro-
gram, regardless of windowing length, is not able to create a joint time-frequency
representation for the system without this uncertainty in resolution.
Figures 2.5 and 2.6 show sliding windows of 256 and 512 points, respectively.
Each window overlaps the previous window by 50%. Choosing discrete windows
reveals the “blocky” nature of the spectrogram method and emphasizes the trade-off
in temporal vs. frequency resolution. The remaining spectrograms, Figures 2.7 – 2.9,
are continuous spectrogram representations for different window sizes, obtained by
shifting the window only one point at a time as the window moves along the original
signal. This creates a smooth distribution, though the smearing of information in the
time-frequency plane is still apparent at all window lengths. It is particularly useful
to compare the linear plots in Figures 2.7 – 2.9. The transient peaks in Figure 2.7
are narrower in the time dimension and wider along the frequency dimension than
the transient peaks in Figure 2.9. This is in agreement with earlier discussions of the
canonical relationship between frequency resolution and temporal resolution. Though
Fine Spectrogram (256 Point Sliding Window)
Fine Spectrogram (256 Point Sliding Window)
0.25
0.25
0.2 0.2
0.15 0.15
Frequency
Frequency
0.1 0.1
35
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
Figure 2.5: Transient test signal, spectrogram, running window of 256 points. (Compare with Figure 2.6.) The information
in the time-frequency plane is blocky and smeared in both time and frequency. The exact onset of the transient signals is
impossible to determine, though it can be estimated with some accuracy.
Coarse Spectrogram (512 Point Sliding Window)
Coarse Spectrogram (512 Point Sliding Window)
0.25
0.25
0.2 0.2
0
0.15 0.15
Frequency
Frequency
0.1 0.1
36
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
Figure 2.6: Transient test signal, spectrogram, running window of 512 points. (Compare with Figure 2.5.) This window
length has much better frequency resolution, but the temporal resolution is very coarse. The transient signals are thereby
smeared across a wider time range, though they are much better represented in terms of frequency content.
Continuous Spectrogram, N/2
Continuous Spectrogram, N/2
0.25
0.25
0
0.2 0.2
0.15 0.15
Frequency
Frequency
0.1 0.1
37
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
Figure 2.7: Transient test signal, Continuous Spectrogram, window of N/2. The continous spectrogram implementation
gives a frequency estimate for each point in the original record. This method still falls prey to the time/frequency resolution
tradeoff of the typical spectrogram (Figures 2.5 and 2.6). The windowing involved with the STFT by definition smears the
arrival of transient signals by a timespan proportional to the windowing function. Therefore, the transient signals appear in
the time-frequency plane before the actual instant of their onset.
Continuous Spectrogram, N/4
Continuous Spectrogram, N/4
0.25
0.25
0
0.2 0.2
0.15 0.15
Frequency
Frequency
0.1 0.1
38
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
Figure 2.8: Transient test signal, Continuous Spectrogram, window of N/4, as in Figure 2.7. For a shorter window, the
temporal resolution increases, though information regarding the transient signals is still presenting before their arrival in the
original signal.
Continuous Spectrogram, N/8
Continuous Spectrogram, N/8
0.25
0.25
0
0.2 0.2
0.15 0.15
Frequency
Frequency
0.1 0.1
39
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
Figure 2.9: Transient test signal, Continuous Spectrogram, window of N/8, as in Figure 2.7. For each successive decrease
in window length, the temporal accuracy has increased at the expense of the accuracy of the frequency estimate in the time-
frequency plane.
40
there is a gain in temporal accuracy for each decrease in window length (N/2 → N/4
→ N/8), the frequency resolution at each level decreases by half as well.
The CWT, shown in Figure 2.10, clearly shows the evolution of the properties of the
signal. This wavelet implementation has more accurate temporal resolution than the
spectrogram methods of Section 2.12.1, and the transient signals are matched in the
time-frequency plane. However, the amplitude of the high-frequency component is
lower than the intermediate component, and there is a diffuse smearing of the time
frequency representation along the frequency axis – the scalogram is biased towards
lower frequencies, and at higher scales the approximation between scale and frequency
becomes less accurate.
Wavelet methods have much better temporal accuracy/support than the spec-
trogram methods, which can be seen in the onset of energy being limited to the
temporal region where the transient signals appear. Even in the logarithmic plot,
which magnifies the expression of energy, the transient signals are correctly located
in the time-frequency plane. The attenuation of high frequencies becomes less ap-
parent in the log plot, as both transient signals stand out clearly from the near-zero
of the field. Some slight smearing is more noticeable for lower frequencies, which
is an artifact of the inexact conversion from “scale” to “frequency.” Of note is the
sudden cutoff for the transient signals, which is not a windowing artifact; this sudden
cutoff is an accurate depiction of the synthetic signal, as it was constructed from five
distinct pieces. The first, third, and fifth segment of the record consist only of the
baseline signal, and it is the excellent temporal resolution of the wavelet that leads
to the crispness of the beginning and ending of these transient signals in the second
and fourth segments.
The diffuse nature of this representation suggests (pseudo-)energy in the time-
frequency plane at frequencies where there is no energy in the original signal. This is
a serious shortcoming for this method despite the crisp temporal resolution.
Continuous Wavelet Transform
Continuous Wavelet Transform
0.25
0.25
0.2 0.2
0.15 0.15
0.1 0.1
41
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
Figure 2.10: Transient test signal, Continuous Wavelet Transform (Morlet Wavelet). The amplitude of this distribution is
unevenly distributed around the known frequency in the signal, biased towards higher frequencies. The wavelet scalogram, in
general, provides excellent temporal resolution. Unfortunately, interpreting the frequency content at each instant is not trivial,
as the basis vectors for a wavelet transformation contain more than one frequency; the representation is limited to a frequency
approximation based on the family of wavelet. The crisp cutoff in the log plot (at 1024 and 2048) is not a windowing artifact,
but reflects the end of the transient signals. Each segment is 512 points, and there is no transient component in the first, third,
or fifth interval.
42
2.12.3 Gabor Transform
The Gabor Transform (Equation 2.36) is a linear transformation that has several
properties similar to the STFT and the CWT. As the Gabor Transform is evenly
sampled in frequency, the Gabor Transform avoids the “spreading” associated with
the CWT at high frequencies. The filterbank-like nature of the Gabor Transform pro-
vides superior temporal localization to the STFT. Figure 2.11 shows good resolution
in the time-frequency plane, though the blocky nature of the Gabor Transform (sim-
ilar to that of a coarsely sampled STFT) reduces the possible maximum resolution.
The Wigner-Ville Distribution for the synthetic signal is shown in Figure 2.12. Inter-
ference from the cross-terms makes this distribution almost unusable for determining
the instantaneous behavior of the signal. Cross-terms in the kernel of the WVD
create interference terms between each component of the signal. These cross terms
appear in the time-frequency plane directly between the auto terms – between each
pair of auto terms is a region of highly oscillatory interference. As the WVD satis-
fies the marginal conditions (Section 2.3), these cross term interference patterns are
completely removed by summing across the time (horizontal) or frequency (vertical)
axis, even though the interference itself is not supported in time or frequency.
0.25 0.25
0.2 0.2
0.15 0.15
Frequency
Frequency
0.1 0.1
43
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
Figure 2.11: Transient test signal, Gabor Transform. The Gabor Transform is equivalent to an ideally-sampled STFT, and
can also be described as a filterbank with linear sampling across the frequency axis. The uniform sampling improves on the
diffuse nature of the CWT and preserves resolution at higher frequencies.
Wigner−Ville Distribution Wigner−Ville Distribution
0.25 0.25
0.2 0.2 0
0.15 0.15
Frequency
Frequency
0.1 0.1
44
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
Figure 2.12: Transient test signal, Wigner-Ville Distribution. In addition to oscillations near the peaks of the signal, there
is cross-term interference located midway between each auto component of the signal – shifted in both time and frequency.
The intermediate “pulse” at time ∼ 1024 and ∼ 0.095Hz is entirely from interference between the two transient peaks. This
interference term is unsupported in either time or frequency: it suggests energy at a frequency where there is known to be no
energy, as shown in the Fourier Transform; and it suggests a temporal onset at an instant when there is no transient signal.
There are other interference terms at ∼ 0.0725Hz between ∼ 256 and ∼ 1536, and at ∼ 0.04067Hz from ∼ 768 to ∼ 2048,
representing interference between the transient terms and the baseline term.
45
baseline amplitude was doubled. Again, these terms are entirely interference, so this
distribution does not accurately represent the character of the signal.
The Reduced Interference Distribution (RID) reduces the expression of the interfer-
ence terms found in the WVD. Reducing the expression of the cross-term interference
is one of the primary goals in using quadratic TFR methods for signal processing.
The RID has a smoothing function in both time (g) and frequency (h). I present
several combinations of time and frequency window lengths in Figures 2.14 – 2.18. A
longer window along the time axis smears information temporally, and a longer win-
dow along the frequency axis increases the bandwidth of the frequency information.
In general form, the RID is similar to the SPWVD, though by adjusting the win-
dowing functions, the RID has better temporal accuracy while avoiding the strongest
expressions of cross-term interference.
In general, the smoothing functions of the RID need to be selected based on each
individual signal being analyzed. There is no combination of window lengths that
performs ideally for all signals, and the TFR method used must be adjusted based
on the frequency content, the character of the signal, etc. However, RID window
lengths of N/4 to N/10 will generally give good results for the signals commonly
encountered in seismology and structural engineering, and in later sections I will be
using windowing functions in this range.
Smoothed Pseudo Wigner−Ville Distribution Smoothed Pseudo Wigner−Ville Distribution
0.25 0.25
0.2 0.2
0.15 0.15
Frequency
Frequency
0.1 0.1
46
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
(a) Smoothed Pseudo Wigner-Ville Distribution (b) Smoothed Pseudo Wigner-Ville Distribution, log10(abs())
Figure 2.13: Transient test signal, Smoothed Pseudo Wigner-Ville Distribution. By eliminating the cross-term interference,
the SPWVD becomes a much more accurate representation of energy in the time-frequency plane. However, this formulation of
the SPWVD brings the WVD closer to a spectrogram representation (c.f. Figure 2.9). The Reduced Interference Distribution
is a related method for removing cross-term interference and in general has superior results (c.f. Figure 2.18).
Reduced Interference Distribution, h=2, g=2 Reduced Interference Distribution, h=2, g=2
0.25 0.25
0.2 0.2
0.15 0.15
Frequency
Frequency
0.1 0.1
47
0.05 0 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
Figure 2.14: Transient test signal, Reduced Interference Distribution, frequency window N/2, time window N/2. The RID
has, in general, better temporal resolution than the SPWVD (Figure 2.13). Some residual interference appears between cross
terms, but the general expression in the time-frequency plane is nearly ideal for this sample signal.
Reduced Interference Distribution, h=4, g=4 Reduced Interference Distribution, h=4, g=4
0.25 0.25
0.2 0.2
0.15 0.15
Frequency
Frequency
48
0.1 0.1
0.05 0 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
Figure 2.15: Transient test signal, Reduced Interference Distribution, frequency window N/4, time window N/4. A shorter
time and frequency windowing function gives better results than Figure 2.14.
Reduced Interference Distribution, h=4, g=10 Reduced Interference Distribution, h=4, g=10
0.25 0.25
0
0.2 0.2
0.15 0.15
Frequency
Frequency
0.1 0.1
49
0.05 0 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
Figure 2.16: Transient test signal, Reduced Interference Distribution, frequency window N/4, time window N/10. For
completeness, I show several similar plots with different time and frequency window lengths. Lengths between N/4 and N/10,
for both time and frequency, generally are suitable for seismic signals.
Reduced Interference Distribution, h=10, g=4 Reduced Interference Distribution, h=10, g=4
0.25 0.25
0.2 0.2
0.15 0.15
Frequency
Frequency
0.1 0.1
50
0
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
Figure 2.17: Transient test signal, Reduced Interference Distribution, frequency window N/10, time window N/4. For
completeness, I show several similar plots with different time and frequency window lengths. Lengths between N/4 and N/10,
for both time and frequency, generally are suitable for seismic signals.
Reduced Interference Distribution, h=10, g=10 Reduced Interference Distribution, h=10, g=10
0.25 0.25
0
0.2 0.2
0.15 0.15
Frequency
Frequency
0.1 0.1
51
0
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
Figure 2.18: Transient test signal, Reduced Interference Distribution, frequency window N/10, time window N/10. This is
the shortest set of time and frequency windowing functions presented in this section. In general, window lengths need to be
adjusted based on the character of the signal, though smoothing functions in the range of N/4 to N/10 are generally useful
starting points.
52
2.12.7 Choi-Williams Distribution
The Choi-Williams Distribution (Figure 2.19), is similar to the RID, and is able to
remove the cross-term interference nearly completely. The temporal behavior matches
the RID, though the crisper temporal onset of the RID makes it better suited for
general signal analysis.
The Born-Jordan Distribution (Figure 2.20) is another smoothed WVD, which re-
moves many cross-terms while keeping the desired temporal and frequency resolution.
0.25 0.25
0.2 0.2
0.15 0.15
Frequency
Frequency
0.1 0.1
53
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
Figure 2.19: Transient test signal, Choi-Williams Distribution. Temporal support is lost, as the transient signals appear before
their instant of onset. Cross term interference, however, is cleanly removed when compared with the Reduced Interference
Distribution (Section 2.12.6).
Born−Jordan Distribution Born−Jordan Distribution
0.25 0.25
0
0.2 0.2
0.15 0.15
Frequency
Frequency
0.1 0.1
54
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
Figure 2.20: Transient test signal, Born-Jordan Distribution. This representation is very similar to the Reduced Interference
Distribution presented in Section 2.12.6. Interference between auto terms is preserved when there are multiple frequency
components at each time instant, which is usually the case in signals of interest to seismology.
Zhao−Atlas−Marks Distribution Zhao−Atlas−Marks Distribution
0.25 0.25
0.2 0.2
0.15 0.15
Frequency
Frequency
55
0.1 0.1
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
Figure 2.21: Transient test signal, Zhao-Atlas-Marks Distribution. This distribution contains more oscillatory terms sur-
rounding the auto terms than the Born-Jordan or Reduced Interference Distributions.
56
impinges on the expression of the auto terms, as they cannot distinguish between
interference and the actual signal in these regions.
The scalogram overlay (scalogram as 2-D filter, Figure 2.22) is biased towards
lower frequencies, and as such does not correctly identify the relative amplitudes
of the two transient signals. The spectrogram overlay (spectrogram as 2-D filter)
is presented for several window sizes (Figures 2.23 – 2.25), and the results are fairly
accurate in the linear plot. The log plot, however, reveals considerable deviation from
the ideal frequency representation for this signal. Interference in the time-frequency
plane is merely reduced, not removed, and the onset of the transient components
still arrive before their true onset in the time domain. These representations fail to
accurately capture the evolution of the signal.
0.25
0.25
0.2 0
0.2
0.15 0.15
Frequency
Frequency
0.1 0.1
57
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
(a) Masked Wigner-Ville Distribution, CWT (b) Masked Wigner-Ville Distribution, CWT, log10(abs())
Figure 2.22: Transient test signal, Masked Wigner-Ville Distribution, Continuous Wavelet Transform (Morlet Wavelet). This
does not remove the interference terms, as the WVD contains cross term oscillations that impinge on the region of auto
terms. In general, other methods (such as the RID, CW, ZAM) reduce cross-term interference first, within the kernel of the
transformation, leading to a better representation for a signal with many closely spaced frequencies.
Masked WVD −− Spectrogram N/2 Masked WVD −− Spectrogram N/2
0.25 0.25
0.2 0.2 0
0.15 0.15
Frequency
Frequency
0.1 0.1
58
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
Figure 2.23: Transient test signal, Masked Wigner-Ville Distribution, Spectrogram N/2. This distribution has many of the
characteristics desired from an ideal time-frequency representation. Using the spectrogram as a 2-D filter suppresses interference
in the time-frequency plane in regions where there is little energy in the signal. The spectrogram is small but nonzero in the
intermediate regions. Since the cross term interference of the WVD can be of large amplitude (in some cases larger than the
auto terms), this method does not properly address the source of the interference, as seen in the plots of log amplitude. This
representation preserves negative values.
Masked WVD −− Spectrogram N/4 Masked WVD −− Spectrogram N/4
0.25 0.25
0.2 0.2 0
0.15 0.15
Frequency
Frequency
59
0.1 0.1
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
Figure 2.24: Transient test signal, Masked Wigner-Ville Distribution, Spectrogram N/4, as in Figure 2.23.
Masked WVD −− Spectrogram N/8 Masked WVD −− Spectrogram N/8
0.25 0.25
0.2 0.2
0
0.15 0.15
Frequency
Frequency
60
0.1 0.1
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
Figure 2.25: Transient test signal, Masked Wigner-Ville Distribution, Spectrogram N/8, as in Figure 2.23.
61
2.12.11 Simple Masking Representation
The continuous wavelet transform (Sections 2.5 and 2.12.2) brings excellent temporal
resolution across nearly the entire frequency range, and the spectrogram (Sections 2.4
and 2.12.1) method gives good temporal and frequency resolution. Combining them
into a simple masking TFR provides an easy way to visually compare the results of
other TFR methods. As it is composed of fundamental techniques that are more
familiar to the engineering community, this example may help to illustrate the types
of compromises necessary in creating a joint time-frequency representation.
This simple TFR has significant drawbacks: most obvious is the attenuation of
higher frequency components. The scale/frequency transformation is an approxima-
tion, as the frequency content of a wavelet is somewhat diffuse and “leaks” across
frequency ranges. Thus, the amplitudes of the two transient signals in the wavelet
plot are different, while the area under the curve (energy integral) is approximately
the same. (The higher frequency component is broader along the frequency axis,
with a smaller peak; the lower frequency component is narrower along the frequency
axis, with a higher peak.) When this is overlain with the spectrogram, the high-
frequency transient signal is decreased relative to the medium-frequency transient –
their amplitudes ought to be roughly the same given the character of the signal (c.f.
the spectrogram plots of Section 2.12.1 and the Fourier Transform of the signal in
Figure 2.4).
Figures 2.26 – 2.28 show Simple Masking TFR plots. The spectrograms are cal-
culated using three different window lengths (N/2, N/4, and N/8), to again demon-
strate the time-frequency resolution tradeoff inherent to the spectrogram method.
The temporal accuracy of the CWT helps to attenuate the Spectrogram energy that
was, non-physically, expressed before the temporal onset of the transient signals. The
Spectrogram, in turn, helps to align the CWT more closely on the actual expressed
frequencies of the signal. However, the shortcomings of both these methods are not
removable by such simple methods. The CWT attenuation of high-frequency energy
is a scaling that extends to the simple TFR, and the temporal lack of support for
62
the Spectrogram is strong enough that the CWT cannot completely suppress the ex-
pression of non-physical energy. The Gabor transform (Sections 2.6 and 2.12.3), as
an ideally sampled spectrogram/scalogram representation, is more accurate than the
simple masking representations.
Simple TFR, CWT Masked with Spectrogram (N/2)
Simple TFR, CWT Masked with Spectrogram (N/2)
0.25
0.25
0.2 0.2
0.15 0.15
Frequency
Frequency
0.1 0.1
63
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
(a) Simple TFR, Scalogram and Spectrogram (N/2) (b) Simple TFR, Scalogram and Spectrogram (N/2), log10(abs())
Figure 2.26: Transient test signal, Simple TFR, Scalogram and Spectrogram (N/2). Combining spectrogram and wavelet
transform methods gives an estimate for the energy in the time-frequency plane. Each method has shortcomings and advantages,
though these simple TFR methods are in general less useful than other established methods, and are presented here as an
illustrative example only. In the linear plot, the temporal resolution seems suitable, though the amplitudes of the transient
components are not equal. The logarithmic scaled plot shows more clearly the smearing along the time-axis.
Simple TFR, CWT Masked with Spectrogram (N/4)
Simple TFR, CWT Masked with Spectrogram (N/4)
0.25
0.25
0.2 0.2
0.15 0.15
Frequency
Frequency
0.1 0.1
64
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
(a) Simple TFR, Scalogram and Spectrogram (N/4) (b) Simple TFR, Scalogram and Spectrogram (N/4), log10(abs())
Figure 2.27: Transient test signal, Simple TFR, Scalogram and Spectrogram (N/4). As in Figure 2.26 for a finer spectrogram.
A shorter window gives better temporal resolution, though there is still energy before the onset of the transient signals. Again,
the wavelet method attenuates energy at high frequencies, so the transient terms do not have the same amplitudes.
Simple TFR, CWT Masked with Spectrogram (N/8)
Simple TFR, CWT Masked with Spectrogram (N/8)
0.25
0.25
0.2 0.2
0.15 0.15
Frequency
Frequency
0.1 0.1
65
0.05 0.05
512 1024 1536 2048 2560 512 1024 1536 2048 2560
Points Points
(a) Simple TFR, Scalogram and Spectrogram (N/8) (b) Simple TFR, Scalogram and Spectrogram (N/8), log10(abs())
Figure 2.28: Transient test signal, Simple TFR, Scalogram and Spectrogram (N/8). As in Figure 2.26 for a finer spectrogram.
At this level of spectrogram, the temporal resolution is more accurate, but still subject to uncertainty considerations. The
Gabor transform, as an ideally sampled STFT (or alternately, as a linearly sampled wavelet filterbank) is in general the best
linear time-frequency transform method, providing the best combination of time and frequency accuracy.
66
2.13 Further Examples
Synthetic signals, with known frequency content, help to illustrate the representations
created by the WVD. A chirp with linearly increasing frequency (Figure 2.29) is a
useful example of a signal with changing frequency content. The WVD of this chirp,
in Figure 2.30, identifies the instantaneous frequency of the original signal, though
some interference is present. The interference is more clearly shown in the logarithmic
plot of Figure 2.31.
The WVD has cross-term interference directly between locations of energy in the
time-frequency plane. This interference is oscillatory in nature, and the Reduced
67
Figure 2.30: Wigner-Ville Distribution of Chirp Function. This is close to the actual
instantaneous frequency content and represents intuitively the correct evolution of the
frequency. Problems include a small amount of noise coming in underneath the signal
(clearer in log plot) and also negative values along the line.
68
Figure 2.32: Double Chirp Function. This is a test function for identifying instanta-
neous frequency content of signals with changing frequencies. Sample Matlab script
to generate this signal is available in Appendix C.
linear chirps separated in the time-frequency plane. The WVD is smoother than
the spectrogram (Figure 2.33), with the shortcoming of interference terms that are
not supported in time or frequency. The RID, however, has excellent resolution in
time and frequency (Figures 2.34 and 2.35), without the large interference terms.
The ability to resolve closely spaced signal components makes the RID suitable for
many signal processing applications, including tracking evolving frequency content of
a nonlinear system.
70
Figure 2.33: Spectrogram and WVD of double chirp function. Neither representa-
tion is accurate – the spectrogram is blocky, and the WVD has interference terms
between the two signals. (Note the oscillatory nature of the interference terms.)
71
Figure 2.34: RID and WVD of double chirp function. The RID has a much clearer
representation of the frequency content – very close to the known values at each
instant. The interference is magnitudes smaller (clearer to see in the 3-D plots of
Figure 2.35). Vertical striations remaining in the RID are small in amplitude and can
be further reduced by refining the smoothing function.
72
Figure 2.35: Double chirp function, 2-D and 3-D Plots of RID and WVD. These
show the improvements of using the RID over the WVD. The 3-D plot presents more
clearly the difference in relative amplitude of the cross term interference.
73
2.13.3 Millikan Library: 2005 Parkfield Earthquake
cm/s2
−5 0
cm/s2
−5
1.5
1.5
1.4
1.4
1.3 0
1.3
1.2 1.2
1.1 1.1
Frequency [Hz]
74
Frequency [Hz]
1 1
0.9 0.9
0.8 0.8
0.7 0.7
50 100 150 200 250 50 100 150 200 250
Time [s] Time [s]
(a) Continuous Spectrogram, N/2 (b) Continuous Spectrogram, window of N/2, log10(abs())
Figure 2.36: Continuous Spectrogram, 2005 Parkfield Earthquake, Millikan Library record. Window of length N/2. For a
longer window, information is smeared across the entire time axis. Energy shows up well before the onset of shaking at t = 50.
Continuous Spectrogram, N/4
Continuous Spectrogram, N/4
5
0 5
cm/s2
−5 0
cm/s2
−5
1.5
1.5
1.4
1.4
1.3 0
1.3
1.2 1.2
1.1 1.1
Frequency [Hz]
Frequency [Hz]
75
1 1
0.9 0.9
0.8 0.8
0.7 0.7
50 100 150 200 250 50 100 150 200 250
Time [s] Time [s]
(a) Continuous Spectrogram, window of N/4 (b) Continuous Spectrogram, window of N/4, log10(abs())
Figure 2.37: Continuous Spectrogram, 2005 Parkfield Earthquake, Millikan Library record. As in Figure 2.36, with a window
of length N/4. There is still energy appearing before the onset of shaking, but the shorter temporal window improves the
representation.
Continuous Spectrogram, N/8
Continuous Spectrogram, N/8
5
0 5
cm/s2
−5 0
cm/s2
−5
1.5
1.5
1.4
1.4
1.3 0
1.3
1.2 1.2
1.1 1.1
Frequency [Hz]
Frequency [Hz]
76
1 1
0.9 0.9
0.8 0.8
0.7 0.7
50 100 150 200 250 50 100 150 200 250
Time [s] Time [s]
(a) Continuous Spectrogram, window of N/8 (b) Continuous Spectrogram, window of N/8, log10(abs())
Figure 2.38: Continuous Spectrogram, 2005 Parkfield Earthquake, Millikan Library record. As in Figure 2.36, window of
length N/8. The onset of shaking at t = 50 is better indicated in the log plot – for a short window, the temporal resolution is
increased at the expense of frequency resolution.
Continuous Wavelet Transform
Continuous Wavelet Transform
5
0 5
−5 0
−5
1.5
1.5
1.4
1.4
0
1.3
1.3
1.2 1.2
1.1 1.1
77
1 1
0.9 0.9
0.8 0.8
0.7 0.7
50 100 150 200 250 50 100 150 200 250
Points Points
Figure 2.39: 2005 Parkfield Earthquake, Millikan Library record. Continuous Wavelet Transform. Wavelet methods are
more diffuse than the spectrogram methods, and it is very difficult to interpret the evolving frequency content of this signal.
However, the onset of shaking at t = 50 is crisply indicated in the log plot, and it would be possible to estimate the peak
frequency at each instant using this method.
Wigner−Ville Distribution Wigner−Ville Distribution
5 5
0 0
cm/s2
cm/s2
−5 −5
1.5 1.5
1.4 1.4
1.3 1.3
0
1.2 1.2
1.1 1.1
Frequency [Hz]
Frequency [Hz]
0
78
1 1
0.9 0.9
0.8 0.8
0.7 0.7
50 100 150 200 250 50 100 150 200 250
Time [s] Time [s]
Figure 2.40: 2005 Parkfield Earthquake, Millikan Library record. Wigner-Ville Distribution. Large-amplitude cross-terms
complicate the interpretation of this distribution. The SPWVD (Figure 2.41) and RID (Figure 2.42) are refinements to the
WVD.
Smoothed Pseudo Wigner−Ville Distribution Smoothed Pseudo Wigner−Ville Distribution
5 5
0 0
cm/s2
cm/s2
−5 −5
1.5 1.5
1.4 1.4
1.3 1.3 0
1.2 1.2
1.1 1.1
Frequency [Hz]
Frequency [Hz]
79
1 1
0.9 0.9
0
0.8 0.8
0.7 0.7
50 100 150 200 250 50 100 150 200 250
Time [s] Time [s]
(a) Smoothed Pseudo Wigner-Ville Distribution (b) Smoothed Pseudo Wigner-Ville Distribution, log10(abs())
Figure 2.41: 2005 Parkfield Earthquake, Millikan Library record. Smoothed Pseudo Wigner-Ville Distribution. The SPWVD
is similar in character to the spectrogram of Figure 2.38. Smoothing the WVD to obtain the SPWVD reduces the expression
of interference terms, but this smoothing also decreases the resolution in the time-frequency plane.
Reduced Interference Distribution, h=4, g=10 Reduced Interference Distribution, h=4, g=10
5 5
0 0
cm/s2
cm/s2
−5 −5
1.5 1.5
1.4 1.4
1.3 1.3
0
1.2 1.2
1.1 1.1
Frequency [Hz]
Frequency [Hz]
1 1
80
0.9 0.9
0.8 0.8
0.7 0.7
50 100 150 200 250 50 100 150 200 250
Time [s] Time [s]
Figure 2.42: 2005 Parkfield Earthquake, Millikan Library record. Reduced Interference Distribution. The RID is not entirely
positive and still suffers from cross-term interference. But (pseudo-)energy is nearly perfectly located in the time-frequency
plane, and the RID is in general a good approximation to the instantaneous frequency content of a signal. The onset of shaking
is clearly indicated at t = 50.
81
In Section 2.3, I presented a brief description of the “marginal conditions,” a
mathematical characterization for the extent of energy in the time-frequency plane.
The Wigner-Ville Distribution, though it contains oscillatory interference and large
negative values, satisfies the marginal conditions as seen graphically in Figures 2.43
and 2.44. Summing across the frequency axis (the frequency integral of the time-
frequency plane) gives the squared envelope of the original time series, and summing
across the time axis (the time integral of the time-frequency plane) gives the squared
Fourier Transform of the original signal. The oscillatory terms impart no net negative
values when integrated along the time or frequency axis. The marginal conditions
are satisfied in part because the interference terms in the WVD are highly oscillatory
– large negative values tend to appear closely spaced with large positive values, such
that if they are summed horizontally or vertically they add to zero. This justifies
smoothing of the WVD as in the SPWVD and RID, where some limited smoothing
of the WVD reduces the expression of interference terms.
Figure 2.43 shows the WVD for the Parkfield event, along with the plots graphi-
cally illustrating the marginal conditions. At the top of the figure is the original time
series, and to the right is the Fourier Transform of the record. Along the left side
of the plot are the squared Fourier Transform and the time integral of the WVD. At
the bottom of the plot are the frequency integral of the WVD and the squared time
series. The integrals and the squared Fourier Transform / Amplitude are presented
separately in Figure 2.44, plotted on the same axis to confirm that the WVD satisfies
the marginal conditions for this record.
82
cm/s2
0
−5
1.5
[Fourier Transform] 2
Fourier Transform
Sum of WVD
1.4
1.3
1.2
Frequency [Hz]
1.1
0
1
0.9
0.8
0.7
Sum of WVD
[Fourier Transform]2
Time Integral of WVD
Figure 2.44: WVD Marginal Conditions, as in Figure 2.43, 2005 Parkfield Earth-
quake at Millikan Library: The time and frequency integrals, for the selected signal,
show more detail as to the match. Matching the time and energy marginals is not an
approximation, but a consequence of the definition of the Wigner-Ville Distribution.
84
2.14 Summary
Time-Frequency Representations such as the Reduced Interference Distribution (RID)
allow for a detailed analysis of evolving signals, with excellent temporal and frequency
resolution. The RID is able to create a point-for-point estimate of the frequency
content in a signal, while reducing the strong interference terms which complicate the
interpretation of the Wigner-Ville Distribution (WVD).
There are various properties that a Time-Frequency Representation (TFR) can
satisfy, some of which (e.g., the marginal conditions) are useful for creating repre-
sentations that match the physical behavior of a system. The TFR methods in this
chapter represent a brief summary of existing techniques. There is no TFR which is
perfectly suited for all tasks; a TFR technique must be selected for each application
based on the type of signal and the information desired about the signal. The modern
TFR methods presented in this thesis are a set of tools that can at times extract more
useful information from a signal than other more widely used techniques such as the
spectrogram and wavelet methods.
85
Chapter 3
Time-Frequency Representations:
Response of Simulated Nonlinear
Systems
The first simulated system is a distributed element model (“Iwan Model,” (Iwan,
1966; Iwan, 1967)), which is a well-studied hysteretic model for a nonlinear-elastic
single degree-of-freedom (SDOF) system with evolving stiffness. Stick-slip elements
fail at specified load levels, causing a loss of stiffness and, therefore, a decreased
natural frequency.
Figure 3.2 presents the response of a distributed element model to a selected input
motion (gaussian noise, first of constant amplitude, then scaled by a cosine envelope).
The hysteretic behavior evolves with time: as the amplitude of motion crosses the
element slip levels, the corresponding spring is lost, weakening the system. For this
system there is hysteretic damping at high amplitudes; once the input function ends,
the system motion is damped out, and eventually the system response returns to
linear oscillations. In Figure 3.3, the same time history (input and response) is shown,
with the estimated secant stiffness for each oscillation of the system. By estimating
the peak-to-peak slope in the force/displacement plane, it is possible to extract an
estimated stiffness for each cycle and correlate the observed natural frequency with
the amplitude of the system response.
Figures 3.4, 3.5, and 3.6 show spectrogram energy estimates for the distributed
element model for both linear and logarithmic scales. The spectrogram windows,
of N/2, N/4, and N/8, respectively, give different levels of temporal and frequency
resolution. For the longest time window (N/2), the evolving stiffness/frequency is
smeared along the time axis. Variations in frequency content are elongated, but this
representation matches the actual stiffness of the system. At the next windowing
level, N/4, the evolving stiffness and instantaneous frequency are qualitatively closer,
echoing the true behavior of the system. The finest window, N/8, has rapid os-
87
Figure 3.1: Distributed Element Model Schematic, after Iwan (1966) and Iwan (1967). Stick-slip elements (“Coulomb Sliders”)
are joined to linear springs – the stick-slip elements fail with increasing amplitudes of motion, at which point the attached
spring contributes no stiffness to the system. This model is a useful analogy for complicated systems with evolving stiffnesses.
88
Hysteretic Behavior
5
Force
−5
−4 −2 0 2 4
Displacement
Figure 3.2: Distributed Element Model, response to enveloped gaussian noise. Input
signal consists of a constant amplitude runup, cosine enveloped gaussian noise, and
zero for the rest of the signal. The hysteretic behavior is defined by a backbone
curve as in 3.1. For shaking in the linear range, there is no damping, and the system
will oscillate perpetually between the minimum and maximum linear portion of the
backbone curve. (For higher levels of input motions, hysteretic damping acts to
dissipate energy until the system reaches the linear range.)
89
Relative stiffness (compared to linear region, secant stiffness at each cycle, from force/displacement hysteresis loop)
1.4
Relative Stiffness
1.2
0.8
0.6
0.4
20 40 60 80 100 120 140 160 180 200
Time [s], Vertical lines at each displacement reversal
Figure 3.3: Distributed Element Model, response to enveloped gaussian noise. For
each loop in the force-displacement plane, it is possible to estimate a peak-to-peak
secant stiffness in the system. This is an approximation for the stiffness that shows
the decrease in stiffness with increasing amplitude and in turn implies a changing
instantaneous natural frequency. The “one-sided” estimate is the final stiffness es-
timated using the distance from the initial conditions. The “two-sided” estimate
uses the slope connecting the positive excursion and the following negative excursion.
These estimates differ in the case when the system does not recover to the pre-event
displacement level due to permanent plastic deformations.
90
cillations in instantaneous frequency – this correctly represents the changing secant
stiffness as the system returns to the linear regime in every cycle.
The continuous wavelet transform (Figure 3.7) shows similar behavior to the N/8
spectrogram. The changing stiffness over each cycle causes rapid changes in the
estimate of instantaneous frequency. Some smearing of energy along the frequency
axis follows from the scale-frequency approximation and the diffuse representation of
high frequency information in the wavelet transformation.
Figure 3.8 shows the same system and the associated Reduced Interference Distri-
bution. The instantaneous frequency estimate again tracks the changing stiffness in
the system. This representation of a SDOF system illustrates the type of connection
between changing stiffness and changing dynamic properties that can be created for
more complicated systems.
Spectrogram, N/2 Spectrogram, N/2, log10()
5 5
0 0
Amplitude
Amplitude
−5 −5
20 40 60 80 100 120 140 160 180 200 20 40 60 80 100 120 140 160 180 200
0.25 0.25
0.2 0.2
91
Frequency [Hz]
Frequency [Hz]
0.15 0.15
0.1 0.1
0.05 0.05
20 40 60 80 100 120 140 160 180 200 20 40 60 80 100 120 140 160 180 200
Time [s] Time [s]
Figure 3.4: Distributed Element Model, Time Series and Spectrogram, window of N/2. Dash-dotted line is linear natural
frequency, and the dark dotted line is the instantaneous frequency estimate for the system. Stiffness of the system, in blue
circles, is estimated from the maximum excursion at each hysteretic loop. A Spectrogram window of N/2 gives coarser temporal
resolution, smearing information along the time axis.
Spectrogram, N/4 Spectrogram, N/4, log10()
5 5
0 0
Amplitude
Amplitude
−5 −5
20 40 60 80 100 120 140 160 180 200 20 40 60 80 100 120 140 160 180 200
0.25 0.25
0.2 0.2
92
Frequency [Hz]
Frequency [Hz]
0.15 0.15
0.1 0.1
0.05 0.05
20 40 60 80 100 120 140 160 180 200 20 40 60 80 100 120 140 160 180 200
Time [s] Time [s]
Figure 3.5: Distributed Element Model, Time Series and Spectrogram, window of N/4. Dash-dotted line is linear natural
frequency, and the dark dotted line is the instantaneous frequency estimate for the system. Stiffness of the system, in blue
circles, is estimated from the maximum excursion at each hysteretic loop. A Spectrogram window of N/4 gives better temporal
information than in Figure 3.4.
Spectrogram, N/8 Spectrogram, N/8, log10()
5 5
0 0
Amplitude
Amplitude
−5 −5
20 40 60 80 100 120 140 160 180 200 20 40 60 80 100 120 140 160 180 200
0.25 0.25
0.2 0.2
93
Frequency [Hz]
Frequency [Hz]
0.15 0.15
0.1 0.1
0.05 0.05
20 40 60 80 100 120 140 160 180 200 20 40 60 80 100 120 140 160 180 200
Time [s] Time [s]
Figure 3.6: Distributed Element Model, Time Series and Spectrogram, window of N/8. Dash-dotted line is linear natural
frequency, and the dark dotted line is the instantaneous frequency estimate for the system. Stiffness of the system, in blue
circles, is estimated from the maximum excursion at each hysteretic loop. A Spectrogram window of N/8 gives fine temporal
resolution, though the resolution in the frequency axis is coarser.
Continuous Wavelet Transform Continuous Wavelet Transform, log10()
5 5
0 0
Amplitude
Amplitude
−5 −5
20 40 60 80 100 120 140 160 180 200 20 40 60 80 100 120 140 160 180 200
0.25 0.25
0.2 0.2
94
Frequency [Hz]
Frequency [Hz]
0.15 0.15
0.1 0.1
0.05 0.05
20 40 60 80 100 120 140 160 180 200 20 40 60 80 100 120 140 160 180 200
Time [s] Time [s]
Figure 3.7: Distributed Element Model, Time Series and Continuous Wavelet Transform (CWT). Dash-dotted line is linear
natural frequency, and the dark dotted line is the instantaneous frequency estimate for the system. Stiffness of the system,
in blue circles, is estimated from the maximum excursion at each hysteretic loop. The wavelet transformation gives a similar
result as the spectrogram for N/8 (3.6), correctly identifying changes in stiffness.
Reduced Interference Distribution Reduced Interference Distribution, log10(abs())
5 5
0 0
Amplitude
Amplitude
−5 −5
20 40 60 80 100 120 140 160 180 200 20 40 60 80 100 120 140 160 180 200
0.25 0.25
0.2 0.2
95
Frequency [Hz]
Frequency [Hz]
0.15 0.15
0.1 0.1
0.05 0.05
20 40 60 80 100 120 140 160 180 200 20 40 60 80 100 120 140 160 180 200
Time [s] Time [s]
(a) Reduced Interference Distribution (b) Reduced Interference Distribution, log10(abs()) of amplitude
Figure 3.8: Distributed Element Model, Time Series and Reduced Interference Distribution. Dash-dotted line is linear natural
frequency, and the dark dotted line is the instantaneous frequency estimate for the system. The instantaneous frequency
estimate closely matches the changing stiffness in the system, though the interference terms complicate identification of the
peak frequency.
96
3.2 Nonlinear Finite Element Building
For the signals of interest to Structural Health Monitoring, the evolution of natural
frequencies becomes important for damage detection and system identification. Small
to moderate earthquakes, for example, will temporarily decrease the apparent natu-
ral frequencies of structures – nonlinearity in the force-displacement relationship will
cause an apparent loss of stiffness with greater excitation levels. This effect can be
seen in buildings under different loading conditions, such as strong winds or forced
vibration testing (e.g., Caltech’s Millikan Library in Chapter B, Clinton (2004), Clin-
ton et al. (2006), and others). Buildings can also be damaged during earthquake
loading, leading to permanent changes in dynamic characteristics.
20
0
−20
0.45
0.4
0.35
Frequency [Hz]
0.3
0.25
0.2
0.15
0.1
0.1
Fracture (N/A) Plastic
0
1 Floors 1−5
Floors 6−10
0
Floors 11−15
−1 Floors 16−20
10
Drift
0
−10
50 100 150 200 250 300
Time (s)
Figure 3.9: RID of synthetic data for the “strong weld” FE model. (Model is not
truly linear, as there is some plastic behavior of the joints. This model was created
without weld fracture, and with material properties that match linear behavior as
well as possible. The small amplitude plots of Section 3.2.1 are closer to true “linear”
building behavior.) Plots as in Figure 3.10.
99
Nonlinear Model:
Roof Displacement, RID, Plastic Rotation, Number of Fractures, and % Drift at selected stories
Displacement
10
0
−10
0.45
0.4
0.35
Frequency [Hz]
0.3
0.25
0.2
0.15
0.1
0.4
Plastic
0.2
40 Floors 1−5
Fracture
Floors 6−10
20 Floors 11−15
0 Floors 16−20
10
Drift
0
−10
50 100 150 200 250 300
Time (s)
Figure 3.10: RID of synthetic data for the nonlinear FE model. The top plot is the
roof displacement, the large center plot is the RID of the data, and the three bottom
plots present the plastic rotation, number of fractured welds, and drift for selected
stories. These damage measures are presented for selected groups of floors, as noted
in the legend.
Nonlinear Model:
Linear Elastic Model: Roof Displacement, RID, Plastic Rotation, Number of Fractures, and % Drift at selected stories
Roof Displacement, RID, Plastic Rotation (small), Number of Fractures (N/A), and % Drift at selected stories
10
20 0
0 −10
Displacement
−20
0.45
Displacement
0.45
0.4
0.4
0.35
0.35
0.3
0.3
0.25
Frequency [Hz]
0.25
Frequency [Hz]
0.2
0.2
0.15
0.15
100
0.1 0.1
0.1 0.4
Floors 1−5 Floors 6−10 Floors 11−15 Floors 16−20
0.2
Plastic
0
1 40 Floors 1−5
Floors 1−5 Floors 6−10 Floors 11−15 Floors 16−20 Floors 6−10
0 20 Floors 11−15
Fracture
−1 0 Floors 16−20
Drift
Drift
0
−10 −10
50 55 60 65 70 75 80 85 90 95 100 50 55 60 65 70 75 80 85 90 95 100
Time (s) Time (s)
Figure 3.11: RID of synthetic data for the “strong-weld” and nonlinear FE model. Plots as in Figure 3.10. The time period
of the initial strong motion pulse is presented to more clearly show the evolution of frequencies on a time scale of seconds. The
change in frequency is much less noticeable for the strong weld model than in the full nonlinear case.
101
of the building based on the change in stiffness in the lower floors; the kink in the
lower floors represents a region of lower stiffness, and as such there are internal waves
generated at the interface between the upper floors and the lower floors. The displaced
shape of the building at each instant shows that the roof displacements are out of
phase with the lower floors of the building during many portions of this record,
particularly between t = 70s and t = 150s. From a parametric standpoint it is
not immediately obvious, a priori, that this mode-splitting would occur. This is
an example of a system where exploratory time-frequency methods can be used to
investigate the behavior of the system and aid in the design of a suitable parametric
model of the evolving dynamic properties.
102
3.2.1 Response to Varying Levels of Input Motion
In order to investigate the connection between shaking levels and damage patterns,
the same input motions as in the previous section were linearly scaled so as to change
the damage levels experienced by the building. Damage patterns at each level of
shaking (Figures 3.12 – 3.18) show the evolution of damage with stronger shaking
and increased damage measures, with the level in percentage representing the scaling
of the ground motions inputs used for the model (from 1% to 150% of the original
shaking). The splitting of the dominant frequency into two frequencies is apparent
at higher scaling levels, though at the largest level of shaking (150% in Figure 3.18),
the permanent offset in the building complicates the interpretation of the results.
For a signal with large offsets, the RID becomes more strongly oscillatory along the
time axis (closely spaced vertical “stripes” in Figure 3.18(a) are a feature of this
interference), and the normalized data is also difficult to interpret. Large amplitudes
at low frequencies, the analog of a Fourier DC-offset, reduce the relative amplitudes
of the RID between 0.1Hz and 0.3Hz, though it is still possible to track the behavior
as similar to that in Figures 3.16 and 3.17. At the strongest levels of shaking, 150%
of the original ground motion amplitudes, the building is within a few percent of total
collapse according to the finite element code. Gravity would likely work to topple the
building in this damage state.
Figures 3.19 and 3.20 summarize the time-frequency behavior of the building under
increasing amplitudes of motion. As the damage grows more severe, the final natural
frequency continues to decrease, and there is a permanent change in frequency at the
end of the event. The time-frequency representation of the roof records provides more
information than simply looking at the Fourier Transform, though there is an obvious
difference in frequency content between the different levels of excitation (the difference
in frequency content can be used as a preliminary investigation tool). The RID adds
further information about the onset of strong shaking, the timespan over which the
frequency decreases and recovers, and identifies splitting in the time-frequency plane
for the unscaled (100%) motions.
Roof Motions, Input Scaled by 0.01 Roof Motions, Input Scaled by 0.01
0.2 0.2
0 0
−0.2 −0.2
Displacement
Displacement
0.45 0.45
0.4 0.4
0.35 0.35
0
0.3 0.3
0
0.25 0.25
Frequency [Hz]
Frequency [Hz]
0.2 0.2
0.15 0.15
0.1 0.1
103
−5 −5
x 10 x 10
3 3
2 2
1 1
Plastic
Plastic
0 0
1 Floors 1−5 1 Floors 1−5
Floors 6−10 Floors 6−10
0 0
Floors 11−15 Floors 11−15
Fracture
Fracture
0.02 0.02
0 0
Drift
Drift
−0.02 −0.02
50 60 70 80 90 100 110 120 130 140 150 50 60 70 80 90 100 110 120 130 140 150
Time [s] Time [s]
Figure 3.12: The input motions of Figure 3.10 were scaled to show the evolution of damage patterns with increasing am-
plitudes. In this case, the amplitudes of 1% resulted in nearly linear behavior, with no fractured welds and and negligible
plasticity.
Roof Motions, Input Scaled by 0.25 Roof Motions, Input Scaled by 0.25
5 5
0 0
−5 −5
Displacement
Displacement
0.45 0.45
0.4 0.4
0.35 0.35
0
0.3 0.3
0.25 0 0.25
Frequency [Hz]
Frequency [Hz]
0.2 0.2
0.15 0.15
104
0.1 −3
0.1 −3
x 10 x 10
4 4
2 2
Plastic
Plastic
0 0
1 Floors 1−5 1 Floors 1−5
Floors 6−10 Floors 6−10
0 0
Floors 11−15 Floors 11−15
Fracture
Fracture
0.5 0.5
0 0
Drift
Drift
−0.5 −0.5
50 60 70 80 90 100 110 120 130 140 150 50 60 70 80 90 100 110 120 130 140 150
Time [s] Time [s]
Figure 3.13: As in Figure 3.12, for scaling of 25% of the original amplitudes. At this amplitude, the plastic behavior begins
to have an effect on the system, but there is still no weld fracture, and the drifts and plasticity are small.
Roof Motions, Input Scaled by 0.50 Roof Motions, Input Scaled by 0.50
10 10
0 0
−10 −10
Displacement
Displacement
0.45 0.45
0.4 0.4
0.35 0.35 0
0.3 0.3
0.25 0.25
Frequency [Hz]
Frequency [Hz]
0.2 0 0.2
0.15 0.15
0.1 0.1
105
0.1 0.1
0.05 0.05
Plastic
Plastic
0 0
15 Floors 1−5 15 Floors 1−5
10 Floors 6−10 10 Floors 6−10
5 Floors 11−15 5 Floors 11−15
Fracture
Fracture
Drift
Drift
−2 −2
50 60 70 80 90 100 110 120 130 140 150 50 60 70 80 90 100 110 120 130 140 150
Time [s] Time [s]
Figure 3.14: As in Figure 3.12, for scaling of 50% of the original amplitudes. Weld fractures on the first 5 floors, some on
floors 6-10, causing a change in the dynamic properties of the roof response. Natural frequency has been permanently changed
by almost 10%, corresponding to a loss of ∼ 19% of the total stiffness.
Roof Motions, Input Scaled by 0.75 Roof Motions, Input Scaled by 0.75
10 10
0 0
−10 −10
−20 −20
Displacement
Displacement
0.45 0.45
0.4 0.4
0.35 0.35 0
0.3 0.3
0.25 0.25
0
Frequency [Hz]
Frequency [Hz]
0.2 0.2
0.15 0.15
0.1 0.1
0.3 0.3
106
0.2 0.2
0.1 0.1
Plastic
Plastic
0 0
Floors 1−5 Floors 1−5
20 Floors 6−10 20 Floors 6−10
10 Floors 11−15 10 Floors 11−15
Fracture
0 Floors 16−20 Fracture 0 Floors 16−20
2 2
0 0
Drift
Drift
−2 −2
−4 −4
50 60 70 80 90 100 110 120 130 140 150 50 60 70 80 90 100 110 120 130 140 150
Time [s] Time [s]
Figure 3.15: As in Figure 3.12, for scaling of 75% of the original amplitudes. Further weld fractures on the first 5 floors,
some on floors 6-10, causing a change in the dynamic properties of the roof response. Natural frequency has been permanently
changed by 26%, corresponding to a loss of 45% of the total stiffness. At this level of excitation, there is a significant change
to the peak of the Fourier Transform (along the Y-axis) compared with the Fourier Transform for lower levels of shaking
(Figures 3.12 – 3.14).
Roof Motions, Input Scaled by 1.00 Roof Motions, Input Scaled by 1.00
20 20
10 10
0 0
−10 −10
Displacement
Displacement
0.45 0.45
0.4 0.4
0
0.35 0.35
0.3 0.3
0
0.25 0.25
Frequency [Hz]
Frequency [Hz]
0.2 0.2
0.15 0.15
0.1 0.1
107
0.4 0.4
0.2 0.2
Plastic
Plastic
Fracture
Fracture
5 5
Drift
Drift
0 0
−5 −5
50 60 70 80 90 100 110 120 130 140 150 50 60 70 80 90 100 110 120 130 140 150
Time [s] Time [s]
(a) Scaled by 100% (Original Ground Motions) (b) Scaled by 100%, Normalized
Figure 3.16: As in Figure 3.12, for the original amplitudes for the synthetic ground motions. Weld fractures on the first 5
floors, some on floors 6-10 and higher, causing a change in the dynamic properties of the roof response. Natural frequency has
been permanently changed by almost 40%, corresponding to a loss of 62% of the total stiffness.
Roof Motions, Input Scaled by 1.25 Roof Motions, Input Scaled by 1.25
20 20
10 10
0 0
−10 −10
Displacement
Displacement
0.45 0.45
0.4 0.4 0
0.35 0.35
0.3 0.3
0
0.25 0.25
Frequency [Hz]
Frequency [Hz]
0.2 0.2
0.15 0.15
0.1 0.1
108
0.4 0.4
0.2 0.2
Plastic
Plastic
0 0
Floors 1−5 Floors 1−5
20 Floors 6−10 20 Floors 6−10
Floors 11−15 Floors 11−15
Fracture
0 Floors 16−20 Fracture 0 Floors 16−20
8 8
6 6
4 4
Drift
Drift
2 2
0 0
50 60 70 80 90 100 110 120 130 140 150 50 60 70 80 90 100 110 120 130 140 150
Time [s] Time [s]
Figure 3.17: As in Figure 3.12, for scaling of 125% of the original amplitudes. Complete weld fractures on the first 5 floors,
with nearly complete weld fractures on floors 6-10. Natural frequency has been permanently changed by approximately 43%,
corresponding to a loss of 67% of the total stiffness. At these amplitudes, the permanent offset in the roof (drift due to plastic
deformation) results in DC-like oscillations at low frequencies, though the evolution of the frequency content is still clear in
the normalized plot.
Roof Motions, Input Scaled by 1.50 Roof Motions, Input Scaled by 1.50
40 40
20 20
0 0
Displacement
Displacement
0.45 0.45
0
0.4 0.4
0.35 0.35
0.3 0.3
0
0.25 0.25
Frequency [Hz]
Frequency [Hz]
0.2 0.2
0.15 0.15
0.1 0.1
0.8 0.8
109
0.6 0.6
0.4 0.4
Plastic
0.2 Plastic 0.2
0 0
Floors 1−5 Floors 1−5
20 Floors 6−10 20 Floors 6−10
Floors 11−15 Floors 11−15
Fracture
Fracture
10 10
5 5
Drift
Drift
0 0
50 60 70 80 90 100 110 120 130 140 150 50 60 70 80 90 100 110 120 130 140 150
Time [s] Time [s]
Figure 3.18: As in Figure 3.12, for scaling of 150% of the original amplitudes. Complete weld fractures on the first 5 floors,
nearly complete weld fracture on floors 6-10. Natural frequency has been permanently changed by 43%, corresponding to a
loss of 67% of the total stiffness. At these amplitudes, the permanent offset in the roof (permanent drift, plastic deformation)
again results in DC-like oscillations at low frequencies. Further work needs to be done in order to create a time-frequency
representation that can properly handle large permanent deformations.
1% (Linear Behavior) 75% (Intermediate Damage) 100%
20
0.2 10 10
0 0
0
cm/s2
cm/s2
cm/s2
!10
!0.2 !10
!20
0.45 0.45 0.45
Frequency (Hz)
Frequency (Hz)
Frequency (Hz)
0.25 0.25 0.25
110
(a) Scaled by 1% (Linear) (b) Scaled by 75%, (Intermediate) (c) Scaled by 100%
Figure 3.19: Summary of evolving damage patterns for increasing levels of input motions and damage. (Compare with
Figure 3.20.) Damage patterns as seen in Figures 3.12, 3.15, and 3.16. In addition to a change in frequency content (the
Fourier Transform along the y-axis), the distribution of energy in the time-frequency plane is significantly different under these
excitation levels. The splitting of the frequency content in (c) represents the different behavior of the roof displacements after
the first few floors are severely damaged.
1% (Linear Behavior) 75% (Intermediate Damage) 100%
20
0.2 10 10
0 0
0
cm/s2
cm/s2
cm/s2
!10
!0.2 !10
!20
0.45 0.45 0.45
0
0.4 0.4 0.4
Frequency (Hz)
Frequency (Hz)
Frequency (Hz)
(a) Scaled by 1% (Linear) (b) Scaled by 75%, (Intermediate) (c) Scaled by 100%
Figure 3.20: Summary of evolving damage patterns for increasing levels of input motions and damage, log10(abs()). (Compare
with Figure 3.19.) Damage patterns as seen in Figures 3.12, 3.15, and 3.16. The logarithmic scaling reveals a more complete
representation of (pseudo-)energy in the time-frequency plane. At each level of excitation the roof record has a different
distribution of energy in the time-frequency plane, from the nearly linear behavior in (a) to the nonlinear behavior in (c).
112
3.3 Conclusions and Discussion
Physical properties of a changing system can be usefully explored through time-
frequency analysis. Losses of stiffness are associated with decreases in natural fre-
quency, and it is possible to identify nonlinear elastic behavior, plastic behavior, and
damage (e.g., weld fractures) by tracking the estimated frequency at each instant.
Distinguishing between nonlinear elastic behavior, plastic behavior, and damage is
not straightforward. If a system’s frequency decreases during strong motions and
then recovers to pre-event levels, it is likely amplitude-related nonlinear elastic be-
havior or plastic behavior with a recovery to the original stiffness. Permanent changes
in the natural frequency are more often correlated with damage.
The simple SDOF of Section 3.1 provides a useful example of using time-frequency
representations to extract information about a system with changing stiffness. A
20-story finite-element building model has many more degrees of freedom than the
simple system, but in Section 3.2 I show that it is possible to extract a similar rela-
tionship between the natural frequencies and the instantaneous stiffness. Increasing
damage and plastic behavior decreases the observed natural frequency, though for
the strongest levels of shaking the information in the time-frequency plane yields
less useful information. This is not a serious shortcoming of using these methods
for damage detection, as at these larger amplitudes it is already immediately ob-
vious from the original time-domain information that large permanent deformations
have occurred and, therefore, detailed time-frequency investigations are not necessary.
The instantaneous behavior of the system allows for an investigation of the fragility
of the building, such as correlating the instant of fracture of an element with the
amplitudes of input motions. In complicated systems such as the 20-story building,
time-frequency representations can also be used to explore building behavior as a first
step to developing a parametric model of the system.
113
Chapter 4
Time-Frequency Representations:
Response of Instrumented
Structures
For instrumented structures, obtaining perfect information about the physical prop-
erties at each instant is impossible. The synthetic systems of Chapter 3 are useful
examples of evolving systems, as the physical properties are perfectly known at each
time instant. In the analysis of the behavior of real systems, physical properties are
typically inferred from the building records. The natural frequencies, mode shapes,
and damping parameters can be estimated from instrument records, either for an en-
tire event of interest or evaluating the record in segments to track the evolution of the
properties. An appropriate parameterization for a system is not always straightfor-
ward, particularly when analyzing a structure with changes in stiffness from damage.
Exploratory methods such as time-frequency analysis are a valuable tool for investi-
gating the behavior of a system and can be used to track evolving frequency content
in an arbitrary signal. These time-frequency methods can then be expanded to pro-
vide estimates of dynamic properties of the system and provide useful insight into the
physical nature of a system – which can in turn lead to a more accurate parameteri-
zation for further system identification.
There are several instrumented civil structures that have undergone strong shak-
ing, some to failure. Applying time-frequency analysis methods to these structures
reveals the instantaneous frequency content in the system, which can then be used to
114
create an estimated damage pattern. The numerical studies of this chapter along with
the synthetic results from Chapter 3 motivate a discussion of how to associate dam-
age with the time-frequency behavior of the system. In general, damaged structures
have a distinctive “fingerprint”’in the time-frequency plane, though further empirical
study is required.
Sample records from several structures are presented with time-frequency anal-
ysis and discussion. I begin with an investigation of Millikan Library (California
Institute of Technology campus, Pasadena, CA), which has experienced several ma-
jor earthquakes since construction in the 1960s. Many of these earthquake events
have had a permanent effect on the natural frequency of the building, as seen in
Figure 1.1. Studies of Millikan’s response to different levels of excitation (e.g., wind,
earthquakes, forced vibration testing) have revealed interesting nonlinear behavior as
well as sensitivity to weather conditions such as rainfall and temperature. With many
potential factors influencing the observed dynamic properties of the building, time-
frequency analysis methods provide a valuable exploratory tool. An introduction to
frequency/amplitude studies is discussed for Millikan Library during the Northridge
Earthquake in order to examine amplitude related nonlinearity during this event.
I also present a brief description of a 52-story office building in Los Angeles, CA,
and the observed behavior during the 1994 Northridge Earthquake. This 52-story
building is a tall, flexible building that remained nearly linear during the Northridge
event and experienced no permanent change in dynamic properties.
The last section is an investigation into the Imperial County Services Building (El
Centro, CA). This instrumented building was damaged (partial collapse) in the 1979
Imperial Valley Earthquake and was later demolished due to the earthquake damage.
In Clinton et al. (2006), the authors present a description of the observed wander of
the natural frequencies of Millikan Library. The natural frequencies are sensitive to
the level of excitation (as seen in, e.g., Kuroiwa (1967) and Appendix B), with higher
amplitudes of motions corresponding to lower observed frequencies. The nonlinear
response between amplitude and observed frequency extends to the response of Mil-
likan Library under wind and earthquake loading, where larger amplitude motions
correspond to increased building response and a lower natural frequency. Millikan
Library is also sensitive to other weather conditions such as rainfall and extreme
temperatures, which complicates inferring changes in structural properties based on
changes in apparent natural frequency.
In the historical behavior plot of Figure 1.1, the data is collected at sparse intervals
throughout decades of investigation. The continuous instrumentation of station MIK
(a triaxial SCEDC station on the 9th floor, see Appendices A and B) has allowed for a
more detailed investigation of the wander in the natural frequencies, as in Figure 4.1.
In this figure, the frequency is estimated from a spectrogram of the continuous data
set, and there is considerable deviation from the mean natural frequency. Apart from
weekly/monthly trends in the frequency, there are temporary changes that occur over
a very short timespan, often associated with earthquakes, forced vibration tests, and
weather patterns, as shown in the bottom plot. An increasing trend over the last
few months of the plot (strongest in the EW mode) corresponds with construction of
117
non-structural partition walls and redistribution of books to repurpose several floors
as office space.
Figure 4.1: Deviation from the mean natural frequency for the 3 fundamental fre-
quencies at Millikan Library station MIK, May 2001 – Nov 2003. Fundamental fre-
quencies and the deviation from average are calculated at each hour. The hourly peak
is shown in thin green lines, the thick black line tracks the daily average. The thick
red horizontal line is the average frequency for the entire time range. Daily rainfall
(black), max wind gust (green), and max temperature (red) from JPL are plotted
along the bottom of the figure. Adapted from Clinton et al. (2006).
Wind storms and rainfall both influence the behavior of Millikan Library, as can
be more clearly seen during a sixty day period during 2003 (Figure 4.2). This time
period had both strong winds and heavy rains. In the top plot is a spectrogram of
the EW component of station MIK, trimmed to show the first EW mode. (The EW
component is the most sensitive to weather effects.) The middle plot is a normalized
spectrogram, scaled by the peak amplitudes, which allows for a clearer investigation
into the evolving natural frequency. The bottom plot shows rainfall and windspeed
for this interval, as measured at the JPL weather station (http://weatherstation.
jpl.nasa.gov/).
118
In the top plot of Figure 4.2, there are several interesting features of the spectro-
gram. Moving from left to right, the strong winds of January 5th and 6th result in
high amplitudes of shaking on those days. There is also a daily cycle in amplitude,
which is tied to the air conditioning equipment; the air conditioner is turned off each
night from midnight to 4:00 a.m., and is responsible for a change in the ambient noise
in the building. A pattern of increasing and decreasing daily amplitude, in intervals
of . . . 5-2-5-2. . . , represents weekday/weekend differences in building use, though this
is somewhat masked by the weather behavior during this timespan. Continuing left
to right, there are noticeable glitches in the data, narrow spikes of high amplitude
shaking, which represent forced vibration tests (February 10th), and a small earth-
quake (22 February 2003, Big Bear M5.4, ∆ = 119km). In the normalized plot it is
easier to identify the changes in frequency correlating with the increased amplitudes
during the windstorm. The overall correlation between rainfall and natural frequency
from February 11th through the 13th also stands out more clearly in the normalized
plot.
Figures 4.3 and 4.4, adapted from Clinton et al. (2006), quantify the observed
changes during the rainfall and windstorm of Figure 4.2. For rainfall, Figure 4.3,
heavy rain increased the natural frequency of the EW and Torsional modes by more
than 3% within a day of the storm, leaving the NS mode largely unaffected. The
second EW mode (not shown) is similarly affected (Clinton et al., 2006). At the end
of the rain event, the system frequency recovers to pre-rain levels over a timespan
of 7-10 days. Though this shift in natural frequency is not fully understood, one
suggested mechanism involves the relation of the foundation to the surrounding soil
at different levels of saturation. Applying Biot’s theory of wave propagation to a shear
beam embedded in a poroelastic half space gives results that qualitatively match the
behavior of the Millikan Library during heavy rainfall (Todorovska and Al Rjoub,
2006a; Todorovska and Al Rjoub, 2006b).
Wind, in Figure 4.4, has the opposite effect on the natural frequency, decreasing
the frequency by 3% during the period of highest windspeeds. The building recovers
from the wind effects immediately, which is consistent with the known nonlinear
119
1.3
Frequency [Hz]
1.2
1.1
MIK: EW Spectrogram
1
1.3
Frequency [Hz]
1.2
1.1
Windspeed [m/s]
Rainfall [mm]
100
20
50
0 0
0 10 20 30 40 50
Time [Days beginning 01−January−2003 00:00GMT]
Figure 4.2: Spectrogram of 60 days of Millikan Library, from January 1st through
February 28th 2003, including strong winds and heavy rainfall. EW Component of
station MIK. In the top plot is the spectrogram, showing the wander of the frequency
content near the first mode of the building. The second plot is a scaled spectrogram,
which more clearly shows the change in observed peak frequency, particularly during
the windstorms of January 5th and 6th, and the rainfall of February 11th through
13th. (Forced vibration testing, February 10th; Big Bear earthquake, M5.4, ∆ = 119,
22 February 2003)
120
Figure 4.3: Deviation from the mean for the natural frequencies of Millikan Library
during February 2003, which includes a major rainstorm. Fundamental frequencies
and the deviation from average are calculated at each hour. The hourly peak is
shown in thin blue lines, the thick dashed black line tracks the daily average. The
dashed red horizontal line is the mean frequency for this time interval. Weather data
is shown at the bottom of the figure, the black bar data is the cumulative hourly
rainfall (re-zeros at midnight). The red line is the maximum hourly temperature,
and the green is the wind gust. At the top of the figure are the natural frequency
amplitudes for the hourly FFT peak. The rainfall coincides with a very sharp rise in
natural frequencies in the EW and Torsional modes, followed by a slow return towards
pre-rainfall levels. Dashed vertical lines represent the start of each new day (12:00
a.m. PST). Frequency spikes are due to instrument glitches (6th, 21st, and 27th of
February), forced vibration testing (February 10th) and the Big Bear earthquake of
22 February 2003. No major increase in excitation amplitude occurs during rainfall
events not associated with high winds. Adapted from Clinton et al. (2006).
121
behavior of the building. Again, the EW and Torsional modes are most strongly
affected. As wind tends to decrease the observed natural frequency and rain tends to
increase the observed frequency, the behavior of the building during a typical storm
(heavy rain and strong winds at the same time) would be difficult to model. The
days selected for this analysis represent windy days with no rainfall and a rainstorm
with unusually low windspeeds, allowing for a separation of wind and rain effects.
Figure 4.4: Deviation from the mean for the natural frequencies of Millikan Library
in January 2003, as in Figure 4.3, for strong winds. Fundamental frequencies and
the deviation from average are calculated at each hour. The hourly peak is shown in
thin blue lines, the thick dashed black line tracks the daily average. The dashed red
horizontal line is the mean frequency for this time interval. Weather data is plotted
at the bottom of the figure. The natural frequencies of the Library dramatically
decrease for the duration of the most intense windstorm, most notably in the EW
direction but also significantly for the NS and torsional modes. Natural frequencies
recover to pre-wind levels quickly. Adapted from Clinton et al. (2006).
Figure 4.5: Deviation from the mean for the natural frequencies of Millikan Library
in August and September of 2002, as in Figure 4.3, for high temperatures. Funda-
mental frequencies and the deviation from average are calculated at each hour. The
hourly peak is shown in thin blue lines, the thick dashed black line tracks the daily
average. The dashed red horizontal line is the mean frequency for this time interval.
Weather data is plotted at the bottom of the figure. These temperatures represent
typical daily variations during the summer months. All 3 natural frequencies increase
during the hottest days. Adapted from Clinton et al. (2006).
123
4.1.2 Millikan Library: 1971 San Fernando Earthquake
In Figure 1.1, I presented the historical behavior of Millikan Library since construc-
tion, and the most striking feature in this plot was the change in natural frequency
during the San Fernando event (9 February 1971, Mw = 6.6, ∆ ∼ 31km), particularly
in the EW component. During this event, there was a permanent loss of stiffness in
the building, which has been suggested to be the result of non-structural damage to
the superstructure and cracking along the foundation/basement levels. The exterior
walls on the East and West side of the building are shear walls that resist NS motions,
and the central elevator core provides further NS stiffness as well as resisting EW
motions. (Please see Appendix B for Millikan structural diagrams.) Non-structural
window panels on the North and South faces of the building provided additional stiff-
ness in the EW direction at the time of construction (though they were not designed
to carry lateral loads). After the San Fernando event there was a small amount of
observed damage to these panels (Jennings, 1971). There is also evidence that some
foundation elements were permanently damaged in the event, such as the concrete
mat at the ground floor and the interface between underground steam tunnels and
the Millikan Library basement.
2
0 0
cm/s
cm/s2
−200 −200
4 4
3.5 3.5
3 3
2.5 2.5
2 2
Frequency (Hz)
Frequency (Hz)
126
1.5 1.5
0
1 1
0.5 0.5
10 20 30 40 50 60 70 80 90 10 20 30 40 50 60 70 80 90
Time (s) Time (s)
Figure 4.6: 1971 San Fernando Earthquake, Millikan Library Response. Channel 136, Roof NS. Pre-event natural frequency
of 1.9Hz is temporarily decreased to less than 1.6Hz, recovering to ∼ 1.7Hz by the end of the record. (Behavior is clearer in
the logarithmic plot of Figure 4.7.)
1971 San Fernando Earthquake, Millikan Library − Channel 136 − Roof, NS Direction, RID 1971 San Fernando Earthquake, Millikan Library − Channel 136 − Roof, NS Direction, Spectrogram
200 200
2
0 0
cm/s
cm/s2
−200 −200
4 4
3.5 3.5
3 3
2.5 2.5
2 2
Frequency (Hz)
Frequency (Hz)
0
127
1.5 1.5
1 1
0.5 0.5
10 20 30 40 50 60 70 80 90 10 20 30 40 50 60 70 80 90
Time (s) Time (s)
Figure 4.7: 1971 San Fernando Earthquake, Millikan Library Response. Channel 136, Roof NS. As in Figure 4.6, for
logarithmic scaling, to more clearly show the location of (pseudo-)energy in the time-frequency plane.
1971 San Fernando Earthquake, Millikan Library − Channel 137 − Roof, EW Direction, RID 1971 San Fernando Earthquake, Millikan Library − Channel 137 − Roof, EW Direction, Spectrogram
200 200
2
0 0
cm/s
cm/s2
−200 −200
4 4
3.5 3.5
3 3
2.5 2.5
2 2
0
Frequency (Hz)
Frequency (Hz)
128
1.5 1.5
1 1
0.5 0.5
10 20 30 40 50 60 70 80 90 10 20 30 40 50 60 70 80 90
Time (s) Time (s)
Figure 4.8: 1971 San Fernando Earthquake, Millikan Library Response. Channel 137, Roof EW. The EW component was
more strongly affected than the NS component: pre-event frequency of near 1.5Hz was temporarily decreased to ∼ 1Hz,
recovering to 1.2Hz by the end of the record. Compare with the logarithmic plots of Figure 4.9.
1971 San Fernando Earthquake, Millikan Library − Channel 137 − Roof, EW Direction, RID 1971 San Fernando Earthquake, Millikan Library − Channel 137 − Roof, EW Direction, Spectrogram
200 200
2
0 0
cm/s
cm/s2
−200 −200
4 4
3.5 3.5
3 3
2.5 2.5
2 2
Frequency (Hz)
Frequency (Hz)
0
129
1.5 1.5
1 1
0.5 0.5
10 20 30 40 50 60 70 80 90 10 20 30 40 50 60 70 80 90
Time (s) Time (s)
Figure 4.9: 1971 San Fernando Earthquake, Millikan Library Response. Channel 137, Roof EW. As in Figure 4.8, for
logarithmic scaling. This allows for a better estimate of instantaneous frequency content during the record.
1971 San Fernando Earthquake, Millikan Library − Channel 138 − Roof, Z Direction, RID 1971 San Fernando Earthquake, Millikan Library − Channel 138 − Roof, Z Direction, Spectrogram
2
50 50
0 0
cm/s
cm/s2
−50 −50
4 4
3.5 3.5
3 3
2.5 2.5
2 2
130
Frequency (Hz)
Frequency (Hz)
1.5 0 1.5
1 1
0.5 0.5
10 20 30 40 50 60 70 80 90 10 20 30 40 50 60 70 80 90
Time (s) Time (s)
Figure 4.10: 1971 San Fernando Earthquake, Millikan Library Response. Channel 138, Roof Z.
1971 San Fernando Earthquake, Millikan Library − Channel 138 − Roof, Z Direction, RID 1971 San Fernando Earthquake, Millikan Library − Channel 138 − Roof, Z Direction, Spectrogram
2
50 50
0 0
cm/s
cm/s2
−50 −50
4 4
3.5 3.5
3 3
2.5 2.5
0
0
2 2
131
Frequency (Hz)
Frequency (Hz)
1.5 1.5
1 1
0.5 0.5
10 20 30 40 50 60 70 80 90 10 20 30 40 50 60 70 80 90
Time (s) Time (s)
Figure 4.11: 1971 San Fernando Earthquake, Millikan Library Response. Channel 138, Roof Z. Logarithmic scaling.
1971 San Fernando Earthquake, Millikan Library − Channel 133 − Basement, NS Direction, RID 1971 San Fernando Earthquake, Millikan Library − Channel 133 − Basement, NS Direction, Spectrogram
2
100 100
0
cm/s
0
cm/s2
−100 −100
4 4
3.5 3.5
3 3
2.5 2.5
2 2
132
Frequency (Hz)
Frequency (Hz)
0
1.5 1.5
1 1
0.5 0.5
10 20 30 40 50 60 70 80 90 10 20 30 40 50 60 70 80 90
Time (s) Time (s)
Figure 4.12: 1971 San Fernando Earthquake, Millikan Library Response. Channel 133, Basement NS.
1971 San Fernando Earthquake, Millikan Library − Channel 133 − Basement, NS Direction, RID 1971 San Fernando Earthquake, Millikan Library − Channel 133 − Basement, NS Direction, Spectrogram
2
100 100
0 0
cm/s
cm/s2
−100 −100
4 4
3.5 3.5
3 3
2.5 2.5
2 2
0
133
Frequency (Hz)
Frequency (Hz)
1.5 1.5
1 1
0.5 0.5
10 20 30 40 50 60 70 80 90 10 20 30 40 50 60 70 80 90
Time (s) Time (s)
Figure 4.13: 1971 San Fernando Earthquake, Millikan Library Response. Channel 133, Basement NS. Logarithmic scaling.
1971 San Fernando Earthquake, Millikan Library − Channel 134 − Basement, EW Direction, RID 1971 San Fernando Earthquake, Millikan Library − Channel 134 − Basement, EW Direction, Spectrogram
2
100 100
0 0
cm/s
cm/s2
−100 −100
4 4
3.5 3.5
3 3
2.5 2.5
2 2
134
Frequency (Hz)
Frequency (Hz)
0
1.5 1.5
1 1
0.5 0.5
10 20 30 40 50 60 70 80 90 10 20 30 40 50 60 70 80 90
Time (s) Time (s)
Figure 4.14: 1971 San Fernando Earthquake, Millikan Library Response. Channel 134, Basement EW.
1971 San Fernando Earthquake, Millikan Library − Channel 134 − Basement, EW Direction, RID 1971 San Fernando Earthquake, Millikan Library − Channel 134 − Basement, EW Direction, Spectrogram
2
100 100
0 0
cm/s
cm/s2
−100 −100
4 4
3.5 3.5
3 3
2.5 2.5
0 0
2 2
135
Frequency (Hz)
Frequency (Hz)
1.5 1.5
1 1
0.5 0.5
10 20 30 40 50 60 70 80 90 10 20 30 40 50 60 70 80 90
Time (s) Time (s)
Figure 4.15: 1971 San Fernando Earthquake, Millikan Library Response. Channel 134, Basement EW. Logarithmic scaling.
1971 San Fernando Earthquake, Millikan Library − Channel 135 − Basement, Z Direction, RID 1971 San Fernando Earthquake, Millikan Library − Channel 135 − Basement, Z Direction, Spectrogram
2
50 50
0 0
cm/s
cm/s2
−50 −50
4 4
3.5 3.5
3 3
2.5 2.5
2 2
136
Frequency (Hz)
Frequency (Hz)
1.5 1.5
1 1
0.5 0.5
10 20 30 40 50 60 70 80 90 10 20 30 40 50 60 70 80 90
Time (s) Time (s)
Figure 4.16: 1971 San Fernando Earthquake, Millikan Library Response. Channel 135, Basement Z.
1971 San Fernando Earthquake, Millikan Library − Channel 135 − Basement, Z Direction, RID 1971 San Fernando Earthquake, Millikan Library − Channel 135 − Basement, Z Direction, Spectrogram
2
50 50
0 0
cm/s
cm/s2
−50 −50
4 4
3.5 3.5
3 3
2.5 2.5
0
0
2 2
137
Frequency (Hz)
Frequency (Hz)
1.5 1.5
1 1
0.5 0.5
10 20 30 40 50 60 70 80 90 10 20 30 40 50 60 70 80 90
Time (s) Time (s)
Figure 4.17: 1971 San Fernando Earthquake, Millikan Library Response. Channel 135, Basement Z. Logarithmic scaling.
138
4.1.3 Millikan Library: 1987 Whittier Narrows Earthquake
Compared with the San Fernando Earthquake, the 1987 Whittier Narrows event (01
October 1987, Mw = 6.1, ∆ ∼ 18km) had slightly smaller roof accelerations in the EW
direction, but a much larger response in the NS component. The Whittier event also
had a smaller amount of permanent change to the dynamic properties, likely because
the San Fernando event had already damaged the most vulnerable components of the
building. The most fragile elements were damaged in the first significant earthquake,
permanently decreasing the system frequency due to loss of stiffness. Subsequent
shaking of the same amplitude will affect the structure less severely.
Figure 4.18 presents the RID and scaled RID for the EW component of the Whit-
tier Narrows event. The natural frequency drops from a pre-event natural frequency
of ∼ 1.21Hz to well below 1Hz during the strongest shaking, but recovers to above
1Hz by the end of the record. Based on subsequent tests of the building, the EW nat-
ural frequency recovered to 1.18Hz in the following months, for a permanent decrease
of only 2.5%.
139
(a) RID of Millikan Library 9th floor records for 1987 Whittier Narrows event.
(b) RID of Millikan Library, scaled by peak amplitude at each time instant.
Figure 4.18: Millikan Library, EW Component of 9th Floor Response during 1987
Whittier Narrows event.
140
4.1.4 Millikan Library: 1991 Sierra Madre Earthquake
The 1991 Sierra Madre earthquake (28 June 1991, Mw = 5.6, ∆ ∼ 19km) was the
next strong shaking that Millikan experienced. The amplitudes were smaller than the
Whittier Narrows event, and there was only a small change in natural frequency. This
is in agreement with the explanation that the most fragile elements were damaged in
the San Fernando event, and subsequent shaking has had a smaller damaging effect.
Figures 4.19 – 4.29 show the RID for this event. The horizontal basement records
(Figures 4.27 and 4.29) show a significant amount of energy at 0.5Hz. This energy
is also represented in the 6th floor and roof records (Figures 4.19 – 4.26). When
analyzing roof records, the observed roof response is a combination of the building
response, the system response (including rocking and soil structure interaction), and
the input motion (which has its own spectrum). A potential weakness of applying
time-frequency analysis techniques to roof records is that it is not immediately pos-
sible to distinguish between the different factors that contribute to changes in the
apparent natural frequencies of the system.
Channel A − Roof, NS Direction Channel A − Roof, NS Direction
200 200
2
0 0
cm/s
cm/s2
−200 −200
4 4
3.5 3.5
3 3
2.5 2.5
2 2
141
Frequency (Hz)
Frequency (Hz)
1.5 1.5
0
1 0 1
0.5 0.5
2 4 6 8 10 12 14 16 18 20 22 2 4 6 8 10 12 14 16 18 20 22
Time (s) Time (s)
Figure 4.19: Millikan Response to 1991 Sierra Madre Earthquake, Channel A – Roof, NS Direction
Channel B − Roof, Z Direction Channel B − Roof, Z Direction
200 200
2
100 100
0 0
cm/s
cm/s2
−100 −100
4 4
3.5 3.5
3 3
2.5 2.5
2 0 2
142
Frequency (Hz)
Frequency (Hz)
1.5 1.5
1 1
0.5 0.5
2 4 6 8 10 12 14 16 18 20 22 2 4 6 8 10 12 14 16 18 20 22
Time (s) Time (s)
Figure 4.20: Millikan Response to 1991 Sierra Madre Earthquake, Channel B – Roof, Z Direction
Channel C − Roof, EW Direction Channel C − Roof, EW Direction
200 200
2
0 0
cm/s
cm/s2
−200 −200
4 4
3.5 3.5
3 3
2.5 2.5
2 2
143
Frequency (Hz)
Frequency (Hz)
0
1.5 1.5
1 0 1
0.5 0.5
2 4 6 8 10 12 14 16 18 20 22 2 4 6 8 10 12 14 16 18 20 22
Time (s) Time (s)
Figure 4.21: Millikan Response to 1991 Sierra Madre Earthquake, Channel C – Roof, EW Direction
Channel D − 6th Floor, NS Direction Channel D − 6th Floor, NS Direction
200 200
2
100 100
0 0
cm/s
cm/s2
−100 −100
4 4
3.5 3.5
3 3
2.5 2.5
2 2
0
144
Frequency (Hz)
Frequency (Hz)
1.5 1.5
1 1
0.5 0.5
2 4 6 8 10 12 14 16 18 20 22 2 4 6 8 10 12 14 16 18 20 22
Time (s) Time (s)
Figure 4.22: Millikan Response to 1991 Sierra Madre Earthquake, Channel D – 6th Floor, NS Direction
Channel E − 6th Floor, Z Direction Channel E − 6th Floor, Z Direction
100 100
2
0 0
cm/s
cm/s2
−100 −100
4 4
3.5 3.5
3 3
2.5 2.5
2 2
145
Frequency (Hz)
Frequency (Hz)
1.5 1.5
1 1
0.5 0.5
2 4 6 8 10 12 14 16 18 20 22 2 4 6 8 10 12 14 16 18 20 22
Time (s) Time (s)
Figure 4.23: Millikan Response to 1991 Sierra Madre Earthquake, Channel E – 6th Floor, Z Direction
Channel F − 6th Floor, EW Direction Channel F − 6th Floor, EW Direction
100 100
2
0 0
cm/s
cm/s2
−100 −100
4 4
3.5 3.5
3 3
2.5 2.5
2 2
146
Frequency (Hz)
Frequency (Hz)
0
0
1.5 1.5
1 1
0.5 0.5
2 4 6 8 10 12 14 16 18 20 22 2 4 6 8 10 12 14 16 18 20 22
Time (s) Time (s)
Figure 4.24: Millikan Response to 1991 Sierra Madre Earthquake, Channel F – 6th Floor, EW Direction
Channel G − Roof, NS Direction Channel G − Roof, NS Direction
200 200
2
0 0
cm/s
cm/s2
−200 −200
4 4
3.5 3.5
3 3
2.5 2.5
2 2
147
Frequency (Hz)
Frequency (Hz)
1.5 1.5
1 0 1
0.5 0.5
2 4 6 8 10 12 14 16 18 20 22 2 4 6 8 10 12 14 16 18 20 22
Time (s) Time (s)
Figure 4.25: Millikan Response to 1991 Sierra Madre Earthquake, Channel G – Roof, NS Direction
Channel H − Roof, EW Direction Channel H − Roof, EW Direction
200 200
2
0 0
cm/s
cm/s2
−200 −200
4 4
3.5 3.5
3 3
2.5 2.5
2 2
148
Frequency (Hz)
Frequency (Hz)
1.5 1.5 0
1 0 1
0.5 0.5
2 4 6 8 10 12 14 16 18 20 22 2 4 6 8 10 12 14 16 18 20 22
Time (s) Time (s)
Figure 4.26: Millikan Response to 1991 Sierra Madre Earthquake, Channel H – Roof, EW Direction
Channel I − Basement, NS Direction Channel I − Basement, NS Direction
100 100
2
0 0
cm/s
cm/s2
−100 −100
4 4
3.5 3.5
3 3
2.5 2.5
2 2
149
Frequency (Hz)
Frequency (Hz)
0
1.5 0 1.5
1 1
0.5 0.5
2 4 6 8 10 12 14 16 18 20 22 2 4 6 8 10 12 14 16 18 20 22
Time (s) Time (s)
Figure 4.27: Millikan Response to 1991 Sierra Madre Earthquake, Channel I - Basement, NS Direction
Channel J − Basement, Z Direction Channel J − Basement, Z Direction
2
50 50
0 0
cm/s
cm/s2
−50 −50
−100 −100
4 4
3.5 3.5
3 3
2.5 2.5
2 2
0
150
Frequency (Hz)
Frequency (Hz)
1.5 1.5
1 1
0.5 0.5
2 4 6 8 10 12 14 16 18 20 22 2 4 6 8 10 12 14 16 18 20 22
Time (s) Time (s)
Figure 4.28: Millikan Response to 1991 Sierra Madre Earthquake, Channel J - Basement, Z Direction
Channel K − Basement, EW Direction Channel K − Basement, EW Direction
100 100
2
0 0
cm/s
cm/s2
−100 −100
4 4
3.5 3.5
3 3
2.5 2.5
2 2
151
Frequency (Hz)
Frequency (Hz)
0
1.5 1.5
1 1
0.5 0.5
2 4 6 8 10 12 14 16 18 20 22 2 4 6 8 10 12 14 16 18 20 22
Time (s) Time (s)
Figure 4.29: Millikan Response to 1991 Sierra Madre Earthquake, Channel K - Basement, EW Direction
152
4.1.5 Millikan Library:
1994 Northridge Earthquake and Aftershock
The 1994 Northridge Earthquake (17 January 1994, Mw = 6.7, ∆ ∼ 34km) had very
large accelerations, and permanently changed the EW and NS natural frequencies
by a small amount. It is interesting to compare the time-frequency plots for the
Northridge event (Figures 4.30 – 4.40) with a Mw = 5.6 aftershock that was recorded
11 hours later (Figures 4.41 – 4.51), particularly when investigating the change in
observed natural frequency.
At the end of the Northridge event, the NS roof records (Figures 4.30 and 4.36)
show an estimated natural frequency of ∼ 1.4Hz, compared with the pre-event fre-
quency of 1.7Hz. The NS natural frequency briefly dropped to 1.3Hz during the
strongest shaking. The exact amount of time that it takes for Millikan to recover
after an event is unknown, but a triggered aftershock of the Northridge event (Fig-
ures 4.41 and 4.47) shows that the pre-event natural frequency, 11 hours after the
Northridge event, is above 1.5Hz. (It is difficult to identify the exact initial natural
frequency of the event due to interference in the RID, but it is unambiguously greater
than 1.5Hz at the beginning of the record.)
In the EW roof components for the Northridge event (Figures 4.32 and 4.37), the
pre-event frequency of 1.18Hz decreases to 0.9Hz during the strongest shaking and
recovers to ∼ 1Hz by the end of the record. The EW aftershock records (Figures 4.43
and 4.48) do not show the same recovering behavior as the NS components; the
resonant frequency during the event is ∼ 1Hz, and it is not obvious that the natural
frequency was larger than 1Hz at the beginning of the record.
Channel A − Roof, NS Direction Channel A − Roof, NS Direction
200 200
2
0 0
cm/s
cm/s2
−200 −200
4 4
3.5 3.5
3 3
2.5 2.5
2 2
153
Frequency (Hz)
Frequency (Hz)
1.5 1.5
0
0
1 1
0.5 0.5
5 10 15 20 25 30 35 5 10 15 20 25 30 35
Time (s) Time (s)
2
0 0
cm/s
cm/s2
−100 −100
4 4
3.5 3.5
3 3
2.5 2.5
2 0 2
0
154
Frequency (Hz)
Frequency (Hz)
1.5 1.5
1 1
0.5 0.5
5 10 15 20 25 30 35 5 10 15 20 25 30 35
Time (s) Time (s)
100 100
2
0 0
cm/s
cm/s2
−100 −100
4 4
3.5 3.5
3 3
2.5 2.5
2 2
155
Frequency (Hz)
Frequency (Hz)
1.5 1.5
0
0
1 1
0.5 0.5
5 10 15 20 25 30 35 5 10 15 20 25 30 35
Time (s) Time (s)
2
100 100
0 0
cm/s
cm/s2
−100 −100
−200 −200
4 4
3.5 3.5
3 3
2.5 2.5
2 2
156
Frequency (Hz)
Frequency (Hz)
1.5 1.5
0 0
1 1
0.5 0.5
5 10 15 20 25 30 35 5 10 15 20 25 30 35
Time (s) Time (s)
Figure 4.33: Millikan Response to Northridge Earthquake, Channel D – 6th Floor, NS Direction
Channel E − 6th Floor, Z Direction Channel E − 6th Floor, Z Direction
50 50
2
0 0
cm/s
cm/s2
−50 −50
4 4
3.5 3.5
3 3
2.5 2.5
0
2 2
0
157
Frequency (Hz)
Frequency (Hz)
1.5 1.5
1 1
0.5 0.5
5 10 15 20 25 30 35 5 10 15 20 25 30 35
Time (s) Time (s)
Figure 4.34: Millikan Response to Northridge Earthquake, Channel E – 6th Floor, Z Direction
Channel F − 6th Floor, EW Direction Channel F − 6th Floor, EW Direction
2
50 50
0 0
cm/s
cm/s2
−50 −50
4 4
3.5 3.5
3 3
2.5 2.5
2 2
158
Frequency (Hz)
Frequency (Hz)
0
1.5 0 1.5
1 1
0.5 0.5
5 10 15 20 25 30 35 5 10 15 20 25 30 35
Time (s) Time (s)
Figure 4.35: Millikan Response to Northridge Earthquake, Channel F – 6th Floor, EW Direction
Channel G − Roof, NS Direction Channel G − Roof, NS Direction
2
200 200
0 0
cm/s
cm/s2
−200 −200
4 4
3.5 3.5
3 3
2.5 2.5
2 2
159
Frequency (Hz)
Frequency (Hz)
1.5 1.5
0
1 1
0.5 0.5
5 10 15 20 25 30 35 5 10 15 20 25 30 35
Time (s) Time (s)
100 100
2
0 0
cm/s
cm/s2
−100 −100
4 4
3.5 3.5
3 3
2.5 2.5
2 2
160
Frequency (Hz)
Frequency (Hz)
1.5 1.5 0
1 1
0.5 0.5
5 10 15 20 25 30 35 5 10 15 20 25 30 35
Time (s) Time (s)
2
0 0
cm/s
cm/s2
−100 −100
4 4
3.5 3.5
3 3
2.5 2.5
2 2
161
Frequency (Hz)
Frequency (Hz)
0
1.5 1.5
1 1
0.5 0.5
5 10 15 20 25 30 35 5 10 15 20 25 30 35
Time (s) Time (s)
2
50 50
0 0
cm/s
cm/s2
−50 −50
4 4
3.5 3.5
3 3
2.5 2.5
0
2 2
162
Frequency (Hz)
Frequency (Hz)
0
1.5 1.5
1 1
0.5 0.5
5 10 15 20 25 30 35 5 10 15 20 25 30 35
Time (s) Time (s)
2
50 50
0 0
cm/s
cm/s2
−50 −50
4 4
3.5 3.5
3 3
2.5 2.5
2 2
0
163
Frequency (Hz)
Frequency (Hz)
0
1.5 1.5
1 1
0.5 0.5
5 10 15 20 25 30 35 5 10 15 20 25 30 35
Time (s) Time (s)
2
0 0
cm/s
cm/s2
−50 −50
4 4
3.5 3.5
3 3
2.5 2.5
2 2
164
Frequency (Hz)
Frequency (Hz)
0
1.5 1.5
0
1 1
0.5 0.5
2 4 6 8 10 12 14 2 4 6 8 10 12 14
Time (s) Time (s)
2
20 20
0 0
cm/s
cm/s2
−20 −20
4 4
3.5 3.5
3 3
2.5 2.5
0
2 2
165
Frequency (Hz)
Frequency (Hz)
0
1.5 1.5
1 1
0.5 0.5
2 4 6 8 10 12 14 2 4 6 8 10 12 14
Time (s) Time (s)
2
0 0
cm/s
cm/s2
−20 −20
−40 −40
4 4
3.5 3.5
3 3
2.5 2.5
0
2 2
166
Frequency (Hz)
Frequency (Hz)
0
1.5 1.5
1 1
0.5 0.5
2 4 6 8 10 12 14 2 4 6 8 10 12 14
Time (s) Time (s)
2
20 20
0 0
cm/s
cm/s2
−20 −20
4 4
3.5 3.5
3 3
2.5 2.5
2 2
167
Frequency (Hz)
Frequency (Hz)
0
1.5 1.5
1 0 1
0.5 0.5
2 4 6 8 10 12 14 2 4 6 8 10 12 14
Time (s) Time (s)
Figure 4.44: Millikan Response to Northridge Aftershock, Channel D – 6th Floor, NS Direction
Channel E − 6th Floor, Z Direction Channel E − 6th Floor, Z Direction
2
10 10
0 0
cm/s
cm/s2
−10 −10
−20 −20
4 4
3.5 3.5
3 3
2.5 2.5
0
2 2
168
Frequency (Hz)
Frequency (Hz)
1.5 1.5
1 1
0.5 0.5
2 4 6 8 10 12 14 2 4 6 8 10 12 14
Time (s) Time (s)
Figure 4.45: Millikan Response to Northridge Aftershock, Channel E – 6th Floor, Z Direction
Channel F − 6th Floor, EW Direction Channel F − 6th Floor, EW Direction
20 20
2
0 0
cm/s
cm/s2
−20 −20
4 4
3.5 3.5
3 3
2.5 2.5
0
2 2
169
Frequency (Hz)
Frequency (Hz)
1.5 0 1.5
1 1
0.5 0.5
2 4 6 8 10 12 14 2 4 6 8 10 12 14
Time (s) Time (s)
Figure 4.46: Millikan Response to Northridge Aftershock, Channel F – 6th Floor, EW Direction
Channel G − Roof, NS Direction Channel G − Roof, NS Direction
50 50
2
0 0
cm/s
cm/s2
−50 −50
4 4
3.5 3.5
3 3
2.5 2.5
2 2
170
Frequency (Hz)
Frequency (Hz)
0
0
1.5 1.5
1 1
0.5 0.5
2 4 6 8 10 12 14 2 4 6 8 10 12 14
Time (s) Time (s)
2
0 0
cm/s
cm/s2
−20 −20
−40 −40
4 4
3.5 3.5
3 3
2.5 2.5
2 2
0
171
Frequency (Hz)
Frequency (Hz)
0
1.5 1.5
1 1
0.5 0.5
2 4 6 8 10 12 14 2 4 6 8 10 12 14
Time (s) Time (s)
20 20
2
0 0
cm/s
cm/s2
−20 −20
4 4
3.5 3.5
3 3
2.5 2.5
2 0 2
172
Frequency (Hz)
Frequency (Hz)
1.5 1.5
1 1
0.5 0.5
2 4 6 8 10 12 14 2 4 6 8 10 12 14
Time (s) Time (s)
20 20
2
0 0
cm/s
cm/s2
−20 −20
4 4
3.5 3.5
3 3
2.5 2.5
2 2
173
Frequency (Hz)
Frequency (Hz)
1.5 1.5
1 1
0.5 0.5
2 4 6 8 10 12 14 2 4 6 8 10 12 14
Time (s) Time (s)
2
0 0
cm/s
cm/s2
−20 −20
4 4
3.5 3.5
3 3
2.5 2.5
0
2 2
174
Frequency (Hz)
Frequency (Hz)
1.5 0 1.5
1 1
0.5 0.5
2 4 6 8 10 12 14 2 4 6 8 10 12 14
Time (s) Time (s)
The Millikan Library has been instrumented since its construction in the 1960s, cre-
ating a valuable database of earthquake records. From these historical records, it
is possible to extract information regarding the nonlinear relationship between the
amplitude of excitation and the natural frequency of the building response. For exam-
ple, in (Clinton et al., 2006), historical behavior of Millikan was compiled into a figure
similar to Figure 4.52. The apparent natural frequency for each event was estimated
from the free response of the building, and a best fit trendline is plotted through the
data. There is a trend towards lower frequencies with higher rooftop accelerations,
implying amplitude-related nonlinearity of the building system (including the foun-
dation and soil-structure interaction). This is consistent with the known behavior of
Millikan during forced vibration testing (Appendix B), where an increase in applied
force results in a decreased natural frequency.
As the Northridge event had large amplitudes with a small permanent offset, I
selected it as a sample event to investigate the relation between instantaneous ampli-
tude and stiffness. The results for Millikan during the Northridge event (Figures 4.55
and 4.56) show a weak trend towards lower frequencies with increasing amplitudes,
though there is a large amount of scatter in the estimated response for each excursion
during the event. This scatter implies that the nonlinearity of Millikan Library under
176
1.5
1.4
Frequency [Hz]
1.3
1.2
1.1
0.9
NS Data
EW Data
0.8
0 1 2 3
10 10 10 10
Rooftop Acceleration [cm/s2]
cm/s
cm/s
−20 −20
2 2
1.8 1.8
1.6 1.6
1.4 1.4
1.2 1.2 0
1 0 1
Frequency (Hz)
Frequency (Hz)
0.8 0.8
178
0.6 0.6
0.4 0.4
0.2 0.2
5 10 15 20 25 30 35 5 10 15 20 25 30 35
Time (s) Time (s)
Figure 4.53: Northridge Earthquake, velocity records, Channel A – Roof, NS Direction. The spectrogram is a smoothed
WVD, which smears information in the time-frequency plane. The RID has better resolution, but interference terms complicate
analysis. Velocity records are typically used for better resolution when calculating instantaneous frequency.
Channel C − Roof, EW Direction Channel C − Roof, EW Direction
10 10
0 0
cm/s
cm/s
−10 −10
−20 −20
2 2
1.8 1.8
1.6 1.6
1.4 1.4
0
1.2 1.2
0
1 1
179
Frequency (Hz)
Frequency (Hz)
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
5 10 15 20 25 30 35 5 10 15 20 25 30 35
Time (s) Time (s)
Figure 4.54: Northridge Earthquake, velocity records, Channel C – Roof, EW Direction. As in Figure 4.53.
Channel A − Roof, NS Direction Channel A − Roof, NS Direction
20 20
0 0
−20 −20
Amplitude [cm/s]
Amplitude [cm/s]
5 10 15 20 25 30 35 5 10 15 20 25 30 35
5 10 15 20 25 30 35 5 10 15 20 25 30 35
50 50
Peak to Peak Peak to Peak
Single Side Single Side
0 0
Frequency [Hz]
Frequency [Hz]
1.2 1.2
5 10 15 20 25 30 35 5 10 15 20 25 30 35
180
1.5 1.5
1.4 1.4
1.3 1.3
Frequency [Hz]
Frequency [Hz]
Peak−to−Peak Peak−to−Peak
One Sided One Sided
1.2 1.2
0 1 0 1
10 10 10 10
Amplitude [cm/s] Amplitude [cm/s]
Figure 4.55: Frequency vs. amplitude relation, Millikan response to Northridge Earthquake, NS Component. The top plot is
the velocity of the record, and the second plot is the amplitude envelope. In the third plot I present peak-to-peak estimates of
amplitude and the single-side amplitude estimate, and the fourth plot is the frequency as estimated from the time-frequency
plane. The last plot is amplitude vs. natural frequency, which reveals a linear trend.
Channel C − Roof, EW Direction Channel C − Roof, EW Direction
10 10
0 0
−10 −10
−20 −20
Amplitude [cm/s]
Amplitude [cm/s]
5 10 15 20 25 30 35 5 10 15 20 25 30 35
5 10 15 20 25 30 35 5 10 15 20 25 30 35
0 0
Frequency [Hz]
Frequency [Hz]
0.8 0.8
5 10 15 20 25 30 35 5 10 15 20 25 30 35
181
1.1 1.1
1 1
0.9 0.9
Frequency [Hz]
Frequency [Hz]
Peak−to−Peak Peak−to−Peak
One Sided One Sided
0.8 0.8
0 1 0 1
10 10 10 10
Amplitude [cm/s] Amplitude [cm/s]
Figure 4.56: Frequency vs. amplitude, Millikan response to Northridge Earthquake, EW Component. The top plot is the
velocity of the record, and the second plot is the amplitude envelope. In the third plot I present peak-to-peak estimates of
amplitude and the single-side amplitude estimate, and the fourth plot is the frequency as estimated from the time-frequency
plane. The last plot is amplitude vs. natural frequency, which is far more diffuse than in the NS component.
Channel A − Roof, NS Direction Channel A − Roof, NS Direction
20 20
0 0
−20 −20
Amplitude [cm/s]
Amplitude [cm/s]
Amplitude Envelope Amplitude Envelope
20 20
15 15
10 10
Amplitude [cm/s]
Amplitude [cm/s]
5 5
Instantaneous Frequency (Estimated from RID) Instantaneous Frequency (Estimated from Spectrogram)
1.5 1.5
1.45 1.45
182
1.4 1.4
1.35 1.35
1.3 1.3
Frequency [Hz]
Frequency [Hz]
1.25 1.25
1.2 1.2
10 15 20 25 30 35 10 15 20 25 30 35
Time [s] Time [s]
Figure 4.57: Frequency vs. amplitude, Millikan response to Northridge Earthquake, NS Component. These plots are selected
portions of the data from Figure 4.55, focused on a shorter timespan. The spectrogram has poorer resolution, and the RID has
interference terms which complicate the interpretation.
Channel C − Roof, EW Direction Channel C − Roof, EW Direction
10 10
0 0
−10 −10
−20 −20
Amplitude [cm/s]
Amplitude [cm/s]
Amplitude Envelope Amplitude Envelope
15 15
10 10
Amplitude [cm/s]
Amplitude [cm/s]
5 5
Instantaneous Frequency (Estimated from RID) Instantaneous Frequency (Estimated from Spectrogram)
1.1 1.1
183
1.05 1.05
1 1
0.95 0.95
0.9 0.9
Frequency [Hz]
Frequency [Hz]
0.85 0.85
0.8 0.8
10 15 20 25 30 35 10 15 20 25 30 35
Time [s] Time [s]
Figure 4.58: Frequency vs. amplitude, Millikan response to Northridge Earthquake, EW Component. These plots are selected
portions of the data from Figure 4.56, focused on a shorter timespan. (As in Figure 4.57.)
184
Figure 4.59: 52-Story Office Building, Los Angeles, CA – photograph and instru-
ment layout. CSMIP Station 24602.
For a tall, flexible building under moderate excitation, the dynamic behavior
should be approximately linear with respect to amplitude. In the time-frequency plane
(Figures 4.60 – 4.63), the fundamental modes of the building do not have a significant
change over the duration of the event. The building is fairly symmetric, and the
natural frequencies in the EW and NS directions are nearly identical; the NS natural
frequencies are within ∼ 2% of the EW natural frequencies for the fundamental
modes and the first two overtones. (By the third overtone, the EW and NS natural
frequencies begin to differ by larger amounts.) The fundamental modes in the EW
185
and NS direction are near 0.16Hz, but the mode the largest response was the first
overtone near 0.54Hz.
During the strongest motions between 15 and 25 seconds, the building changed
from the pre-event natural frequencies to a slightly lower apparent frequency (clearest
in the first overtones near 0.55Hz). This change (of less than 1%) is much smaller
than the observed changes in the other systems in this thesis, as the building behaved
nearly linearly for the Northridge Earthquake.
1994 Northridge Earthquake, 52 Story Office Building, Los Angeles − Roof, EW Direction, RID 1994 Northridge Earthquake, 52 Story Office Building, Los Angeles − Roof, EW Direction, Spectrogram
2
100 100
0 0
cm/s
cm/s2
−100 −100
−200 −200
2 2
1.8 1.8
1.6 1.6
1.4 1.4
1.2 1.2
1 1
Frequency (Hz)
Frequency (Hz)
186
0.8 0.8
0
0.6 0.6
0.4 0.4
0.2 0.2
20 40 60 80 100 120 140 160 180 20 40 60 80 100 120 140 160 180
Time (s) Time (s)
Figure 4.60: 1994 Northridge Earthquake – 52-Story Office Building, Channel 19, Roof EW. This building behaves nearly
linearly for this event, though there is a small change in frequency between 15 and 25 seconds (clearest in the first overtone).
1994 Northridge Earthquake, 52 Story Office Building, Los Angeles − Roof, EW Direction, RID 1994 Northridge Earthquake, 52 Story Office Building, Los Angeles − Roof, EW Direction, Spectrogram
2
100 100
0 0
cm/s
cm/s2
−100 −100
−200 −200
2 2
1.8 1.8
1.6 1.6
1.4 1.4
1.2 1.2
0
1 1
Frequency (Hz)
Frequency (Hz)
0
187
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
20 40 60 80 100 120 140 160 180 20 40 60 80 100 120 140 160 180
Time (s) Time (s)
Figure 4.61: 1994 Northridge Earthquake – 52-Story Office Building, Channel 19, Roof EW. As in Figure 4.60, with logarithmic
scaling to more clearly show the information in the time-frequency plane and compare the RID with the Spectrogram.
1994 Northridge Earthquake, 52 Story Office Building, Los Angeles − Roof, NS Direction, RID 1994 Northridge Earthquake, 52 Story Office Building, Los Angeles − Roof, NS Direction, Spectrogram
200 200
2
0 0
cm/s
cm/s2
−200 −200
2 2
1.8 1.8
1.6 1.6
1.4 1.4
1.2 1.2
1 1
Frequency (Hz)
Frequency (Hz)
188
0
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
20 40 60 80 100 120 140 160 180 20 40 60 80 100 120 140 160 180
Time (s) Time (s)
Figure 4.62: 1994 Northridge Earthquake – 52-Story Office Building, Channel 20, Roof NS. The NS natural frequencies are
very similar to the EW frequencies (Figure 4.60).
1994 Northridge Earthquake, 52 Story Office Building, Los Angeles − Roof, NS Direction, RID 1994 Northridge Earthquake, 52 Story Office Building, Los Angeles − Roof, NS Direction, Spectrogram
200 200
2
0 0
cm/s
cm/s2
−200 −200
2 2
1.8 1.8
1.6 1.6
1.4 1.4
1.2 1.2
0
1 1
Frequency (Hz)
Frequency (Hz)
0
189
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
20 40 60 80 100 120 140 160 180 20 40 60 80 100 120 140 160 180
Time (s) Time (s)
Figure 4.63: 1994 Northridge Earthquake – 52-Story Office Building, Channel 20, Roof NS. As in Figure 4.62, with logarithmic
scaling to more clearly show the information in the time-frequency plane and compare the RID with the Spectrogram.
190
4.3 Imperial County Services Building, El Centro,
CA – (1979 Imperial Valley Earthquake)
Time-frequency representations allow for a detailed investigation into the evolving
frequency content of records from structures. As in Chapter 3, one possible use for
TFR methods such as the Wigner-Ville Distribution is to correlate recorded motions
with changes to the physical properties of an instrumented building.
One instrumented building that has experienced strong earthquake motions is the
Imperial County Services Building in El Centro (Figure 4.64); this six-story reinforced
concrete building was famously damaged during the 1979 Imperial Valley earthquake,
and later was demolished due to partial collapse during the event. At the time of
the event, the building was instrumented with 13 accelerometers on 4 levels in the
building and 3 channels at a free-field site (http://www.consrv.ca.gov/cgs/smip/).
In Figures 4.65(a) and 4.65(b), the building suffers a significant change in natural
frequency between 6 and 12 seconds. An ambient vibration test of the building
several months before the event identified the NS fundamental frequency at 2.24Hz,
and the observed natural frequency drops to 1.5Hz at 7 seconds into the record (a
decrease of 33%). From t = 7s to t = 11s, the natural frequency drops to less than
0.9Hz (a further loss of 27% for a total loss of 60% of the original frequency). The
frequency recovered slightly to 1.1Hz by the end of the record, for a final drop of
p
∼50%. As natural frequency, ω, is proportional to k/m (k = stiffness and m =
mass), a decrease in 50% in natural frequency corresponds to a loss of 75% of the
global system stiffness. Most of the damage, and therefore most of the loss of stiffness,
was due to damage in the first story, including partial collapse of a column.
Past investigations into this building have suggested that damage was initiated
∼ 7 seconds into the event, and columns collapsed ∼ 11 seconds into this record
(Rojahn and Mork, 1982; Todorovska and Trifunac, 2006a; Todorovska and Trifunac,
2006b). This agrees with the information as presented in the time-frequency plane
(Figure 4.65(b)), with an evolving natural frequency that agrees with the general
known damage pattern of the event.
191
Figure 4.64: Imperial County Services Building, El Centro, CA – photograph and instrument layout. CSMIP Station 01260.
192
(a) RID of Imperial County Services Building roof records for 1979 Imperial
Valley event.
(b) RID of Imperial County Services Building records, scaled by peak amplitude
at each time instant. This is an estimate for instantaneous frequency, and
can be tied to the behavior of the structure, such as partial collapse between
t = 7s and t = 11s.
Figure 4.65: Imperial County Services Building (El Centro, CA). NS Component
of roof response during 1979 Imperial Valley Event.
193
4.4 Conclusions and Discussion
Modern time-frequency representations allow for a detailed investigation into the be-
havior of an evolving signal. When those signals are from instrumented structures,
instantaneous frequency content can be correlated with physical changes in a struc-
ture. For damaged buildings, time-frequency analysis methods such as the Reduced
Interference Distribution provide increased resolution in the time-frequency plane,
which allows for a detailed, instant-for-instant analysis of the building behavior.
The 52-story office building (Los Angeles, CA), by contrast with Millikan Library,
remains nearly linear during the moderate excitations of the 1994 Northridge Earth-
quake. With tall flexible buildings, time-frequency analysis methods can help to
constrain the dynamic properties.
Instrumented buildings under strong earthquake excitation are rare, but one well-
studied example is the Imperial County Services Building, El Centro, CA. This build-
ing was severely damaged in the 1979 Imperial Valley Earthquake and was later de-
molished. As the building was instrumented during the event, it is possible to use
time-frequency analysis methods to investigate the changing dynamic properties of
194
the building response and compare the observed behavior in the time-frequency plane
with previous investigations into the onset and extent of damage.
In all these cases, time-frequency analysis methods allow for an exploration of
building behavior, adding another tool to the evaluation and identification of build-
ing properties. The Reduced Interference Distribution adds resolution in the time-
frequency plane when compared with spectrogram and wavelet methods and provides
a more exploratory tool as compared to parametric system identification techniques.
195
Chapter 5
Conclusions
Time-frequency analysis methods allow for the investigation of signals with evolv-
ing frequency content. The Wigner-Ville Distribution (including related techniques
such as the Reduced Interference Distribution) allows for higher resolution in the
time-frequency plane compared with Fourier and wavelet methods, within the con-
straints of the uncertainty principle. This increase in resolution makes advanced
time-frequency analysis methods a useful tool for investigating the dynamic proper-
ties of instrumented systems.
Various sample signals in Chapter 2 illustrate the advantages and shortcomings
of the time-frequency representations described in this thesis. The smearing of the
spectrogram and the diffuse nature of the wavelet transforms complicate the interpre-
tation of these time-frequency representations. The Reduced Interference Distribu-
tion overcomes many of these complications, though it introduces shortcoming such
as interference and negative values. There is no single ideal time-frequency repre-
sentation for every application; the operator needs to select a representation that
most clearly reveals the information of interest for a given signal. In the case of
instrumented structures under earthquake excitation, the Reduced Interference Dis-
tribution is nearly optimal for revealing information about the onset and extent of
changes in observed apparent frequency. The changes in frequency content of the
signal can then be correlated with physical changes to the system.
In Chapter 3 I presented synthetic systems under various types of excitation
and used time-frequency representations to explore the dynamic behavior of the
196
system response and correlate changes in the system with changes in the observed
frequency content of the signal. These synthetic systems demonstrate the value of
time-frequency analysis techniques for analyzing signals that correspond to changing
physical systems. Transform methods provide an essential framework for investigating
the properties of these systems, particularly when the stiffness of the system changes.
Chapter 4 provides a more detailed investigation into instrumented systems of
interest, including Caltech’s Millikan Library. Millikan is a useful testbed for investi-
gation of nonlinear systems, as it has been instrumented for more than 40 years and
has well-established nonlinear behavior under small to moderate earthquake excita-
tion as well as sensitivity to weather conditions such as wind, rain, and temperature.
Appendix B presented a study of the nonlinear behavior of Millikan Library during
a forced vibration experiment, where increasing the applied force through a different
shaker configuration led to a decreased natural frequency (though the mechanism
for exciting the building with a roof shaker and with earthquake excitations at the
ground level are different, the general trend is consistent). This nonlinearity agrees
with historical records showing a trend towards decreased natural frequency with
increased amplitude of motion, both in earthquake records and in connection with
strong winds. Further investigating the nonlinearity using the Reduced Interference
Distribution allows for an instant-for-instant description of observed frequency with
respect to changing amplitudes of excitation. A summary of digitized records, using
modern time-frequency techniques, is presented in order to help build a catalog of
building records under earthquake excitation. The collection of Millikan and other
building records is part of an attempt to empirically derive a distinctive “fingerprint”
of damage in the time-frequency plane. The character of the time-frequency repre-
sentation will in turn aid in developing parametric models of the system.
Changes in the physical properties of a system (from nonlinearity, loss of stiffness,
or damage) are represented in the time-frequency plane as changes in frequency con-
tent. The Wigner-Ville Distribution provides a mathematical framework to describe
evolving frequencies in a signal, where the signal can be a multi-component synthetic
(as in Chapter 2), a simulation of a nonlinear system (Chapter 3), or records from an
197
instrumented building (Chapter 4). Time-frequency analysis techniques are a valuable
exploratory tool, providing guidance for the development of parametric models. Fur-
ther work in connecting parametric and time-frequency analysis techniques can help
build an empirical and theoretical connection between changing frequency content of
a signal and changing stiffness of a system.
198
199
Appendix A
1
(4g)( ) = 2.38 × 10−7 g = 2.34 × 10−4 cm/s2 . (A.1)
224
200
As typical ambient vibrations in an instrumented building can range from ±0.1cm/s2
to ±0.01cm/s2 , there is an effective signal-to-noise ratio of ∼ 50 to 500 for the analysis
of ambient building behavior.
The value of 24-bit digitization is demonstrated in Figure A.1 for sample data from
Caltech’s Millikan Library (station MIK). MIK is digitized at 24bits, and the data is
shown for the equivalent 20-bit and 19-bit digitization. At 19-bits (32x coarser than
the 24-bit data), the signal is essentially zero-crossings, as the amplitude is roughly
equal to the digitization cutoffs. 20-bit data (16x coarser than the 24-bit data, 2x
finer than the 19-bit data) is less blocky, though the discrete cutoffs still harm the
character of the signal.
There are five instrumented buildings on the Caltech campus, with various instrument
types. The five buildings are the Robert A. Millikan Memorial Library (“Millikan Li-
brary”), the Broad Center for the Biological Sciences (“Broad Center”), a telescope
pit beneath the Robinson Laboratory (“Robinson Pit”), the Caltech Athenæum, and
the basement of the United States Geological Survey offices at 525 S. Wilson Avenue.
Of these, only the Broad Center and Millikan Library have instruments through-
out the structure. The other stations are located beneath the buildings. The Mil-
likan Library and Broad Center continuously record data to the Southern California
Earthquake Data Center (SCEDC), and waveforms can be freely downloaded from
http://www.data.scec.org/STP/stp.html (data from other stations can be down-
loaded as described in the following sections).
Caltech’s Millikan Library (Figure A.2), one of the world’s most heavily researched
and instrumented buildings, has been studied extensively since its construction in the
1960s. It is a very stiff, 9-story, reinforced concrete building with one basement level.
201
Counts
0 0
−96 −96
1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 1020 1040 1060 1080 1100 1120 1140 1160 1180 1200
Station MIK, 20−bit data (16x coarser than 24−bit data) Station MIK, 20−bit data (16x coarser than 24−bit data)
96 96
80 80
64 64
48 48
32 32
Counts
Counts
16 16
0 0
−16 −16
−32 −32
−48 −48
−64 −64
−80 −80
−96 −96
1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 1020 1040 1060 1080 1100 1120 1140 1160 1180 1200
Station MIK, 19−bit data (32x coarser than 24−bit data) Station MIK, 19−bit data (32x coarser than 24−bit data)
96 96
64 64
32 32
Counts
Counts
0 0
−32 −32
−64 −64
−96 −96
1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 1020 1040 1060 1080 1100 1120 1140 1160 1180 1200
Points [20sps data] Points [20sps data]
(a) Digitization example, 100s (2000 points) of (b) Digitization example, 10s (200 points) of
MIK data MIK data
Figure A.1: Digitization example. Data from MIK (02-Nov-2005), showing 24-
bit data and equivalent data as would have been obtained from a 20-bit and 19-bit
digitization system for the same instruments. The 20-bit data is somewhat coarser
than the 24-bit data, as the bit level resolution is 16x less than the original data. The
19-bit data is 32x coarser than the 24-bit data and is essentially providing only the
zero-crossing information. (Gain for station MIK is 4280counts/(cm/s2 ). ±96counts
= ±0.022cm/s2 .)
202
A synchronized vibration generator (“shaker”) is installed on the roof for forced vibra-
tion tests. On the 9th floor is SCEDC station MIK: a triaxial Episensor accelerometer
and 24-bit Quanterra datalogger. Three channels of data from station MIK are con-
tinuously telemetered to the SCEDC archives for permanent storage of 20sps data
and 80sps data for triggered events.
The USGS/Caltech Dense Instrumentation Network was installed in 1998. This
dense network consists of 36 instruments: three horizontal accelerometers on each
floor (from basement to roof) and three vertical accelerometers in the basement.
The dense array is digitized through two parallel systems: a 19-bit Mt. Whitney
system run by USGS and a 16-bit Digitexx Datalogger (a proprietary digitization
and telemetry system).
The Broad Center (Figure A.3), completed in 2002, is a three story steel building
with two deep basement levels. Unbonded braced frames in the North-South and
East-West direction form an interior core. SCEDC station CBC was installed during
construction and has been recording data since 2003. Station CBC consists of three
triaxial Episensors and a 24-bit Quanterra datalogger. The instruments are located
at the North-West corner of the internal braced frame core at the ground floor, the
North-West corner of the core at the roof level, and the South-East corner of the core
at the roof level. (The South-East roof sensor records the horizontal components and
discards the vertical information, as the digitizer can only record eight channels of
data.) These 8 channels are continuously telemetered to the SCEDC archives and
stored at 20sps (100sps for triggered events).
In the basement of the Robinson Building is an unused telescope pit. This pit houses
SCEDC station CRP, which is currently used as a testing station, so the instrumen-
tation varies.
203
Figure A.2: Caltech Millikan Library. View from the South-East. The dark grey
walls in the foreground are structural shear walls, and the white panels on the south
(left) side of the building are non-structural marble cladding.
204
Figure A.3: Caltech Broad Center. View from the South-West. The South wing of
the building has an irregular floorplan with increasing height; the main building is to
the North.
SCEDC Station CAC, located in the basement of the Athenæum faculty club on the
Caltech campus, is digitized through a 19-bit K-2 datalogger. Only triggered data
from this station is stored in the SCEDC archives.
The USGS Pasadena office is located at 525 S. Wilson Ave, across the street from
Caltech (200m due West of Millikan Library and 300m SSW of the Broad Center). In
the basement is SCEDC station GSA (24-bit Quanterra digitizer and triaxial Episen-
sor accelerometer), which is commonly used as a reference station for data from other
Caltech stations.
The Caltech Online Monitoring and Evaluation Testbeds (COMET) site, http://
comet.caltech.edu, was developed as an educational tool as well as a testbed for
the development of structural health monitoring techniques.
The COMET site acts as an interface between the Southern California Earthquake
Data Center (SCEDC) database and a web-browser, using an implementation of the
205
Earthworm Processing System. Real-time data from two campus buildings, Millikan
Library (SCEDC station MIK) and the Broad Center (SCEDC station CBC), is
archived from the SCEDC at 80sps and 100sps, respectively, and stored on a local
hard drive. COMET users can log onto the site and download selected data, display
these data, and perform modal analysis using Caltech’s MODE-ID modal analysis
software (Beck, 1996). The COMET interface was implemented by Dr. H. F. Lam,
City University of Hong Kong, during a postdoctoral scholarship at Caltech.
The Doris and Louis Factor Building (“Factor Building”) is a 15-story moment-
resisting steel frame building on the UCLA Campus. The Factor building has 72
channels of continuous accelerometer data, four horizontal instruments on each floor,
plus two horizontal instruments and two vertical instruments on each of the two
basement levels. Instruments are uniaxial FBA-11s, digitized with 24-bit Quanterras.
Figure A.4 shows the distribution of instruments and their orientations. There is
also a 100m borehole located approximately 25m east of the Factor Building, with
triaxial Episensors at the top and the bottom of the borehole, telemetered through a
Quanterra digitizer.
This dense instrumentation network makes the Factor Building one of the most
heavily-instrumented buildings in the world. Data is stored in a buffer at 500sps
and downsampled to 100sps for permanent storage. The 100sps Factor Building
data is freely available from the IRIS Data Management Center, http://www.iris.
edu/. For further information on the Factor Building, please see http://factor.
gps.caltech.edu and Kohler et al. (2006).
A.3 Summary
(a) UCLA Factor Building, view from the North-East (b) Layout of instruments from the UCLA Factor Building
Figure A.4: UCLA Factor Building. (Instrument layout from http://factor.gps.caltech.edu.) The Factor Building is a
moment frame building with an irregular floorplan that increases in area at the 10th floor. This causes an internal reflection
of traveling waves.
207
pendent on the advances in seismic instrumentation and digitization; the 24-bit sys-
tems in common use allow for investigations of ambient vibrations as well as moderate
to strong earthquake motions. Several buildings on the Caltech and UCLA campuses
have been instrumented, and selected data is further analyzed in this thesis.
208
209
Appendix B
Adapted from Bradford, Clinton, Favela, and Heaton 2004, Earthquake Engineering
Research Library report on research supported by the California Institute of Tech-
nology.
ACKNOWLEDGEMENTS
The authors would like to acknowledge Arnie Acosta for data triggering and retrieval
and thank him for his support. We thank Caltech’s Structural Monitoring Group for
their input during this project, and we also thank the Southern California Earthquake
Center and the Portable Broadband Instrumentation Center at the University of
California Santa Barbara for the loan of the portable instrument. We acknowledge
the SCEDC for the MIK data.
210
ABSTRACT
This report documents an investigation into the dynamic properties of Millikan Li-
brary under forced excitation. On July 10, 2002, we performed frequency sweeps from
1 Hz to 9.7 Hz in both the East-West (E-W) and North-South (N-S) directions using
a roof level vibration generator. Natural frequencies were identified at 1.14 Hz (E-W
fundamental mode), 1.67 Hz (N-S fundamental mode), 2.38 Hz (Torsional fundamen-
tal mode), 4.93 Hz (1st E-W overtone), 6.57 Hz (1st Torsional overtone), 7.22 Hz (1st
N-S overtone), and at 7.83 Hz (2nd E-W overtone). The damping was estimated at
2.28% for the fundamental E-W mode and 2.39% for the N-S fundamental mode. On
August 28, 2002, a modal analysis of each natural frequency was performed using
the dense instrumentation network located in the building. For both the E-W and
N-S fundamental modes, we observe a nearly linear increase in displacement with
height, except at the ground floor, which appears to act as a hinge. We observed
little basement movement for the E-W mode, while in the N-S mode 30% of the roof
displacement was due to basement rocking and translation. Both the E-W and N-S
fundamental modes are best modeled by the first mode of a theoretical bending beam.
The higher modes are more complex and not well represented by a simple structural
system.
211
B.1 Introduction
This report documents the forced vibration testing of the Robert A. Millikan Memo-
rial Library (“Millikan Library”) located on the California Institute of Technology
campus. It also provides a historical backdrop to put our results in perspective.
During and immediately after the construction of the library in the late 1960s,
numerous dynamic analyses were performed (Kuroiwa, 1967; Foutch et al., 1975; Tri-
funac, 1972; Teledyne-Geotech-West, 1972). In these analyses, fundamental modes
and damping parameters were identified for the library, and higher order modes were
suggested, but not investigated. It has been established that the fundamental fre-
quencies of the library vary during strong motion (Luco et al., 1986; Clinton et al.,
2006). Some drift in the long-term behavior of the building has also been observed in
compiled reports of modal analysis from the CE180 class offered every year at Caltech
(Clinton, 2004).
The temporal evolution of the building’s dynamic behavior, as well as the much
improved density and quality of instrumentation, led to an interest in a complete dy-
namic investigation into the properties of the system. Our experiments were designed
to provide an updated account of the fundamental modes and to identify and explore
the higher order modes and modeshapes. A better understanding of the dynamic
behavior of the Millikan Library will aid in the increased research being performed on
the building and provide a better understanding of the data currently being recorded
by the instruments in the library.
An initial test was performed on July 10, 2002, for which we performed a full
frequency sweep of the building (from 9.7Hz, the limit of the shaker, to 1Hz), in both
the E-W and N-S directions. Frequencies of interest were explored in more detail,
with a finer frequency spacing and different weight configurations, during a second
test on August 28, 2002.
212
B.2 Millikan Library
The Millikan Library (Figure B.1) is a nine-story reinforced concrete building, ap-
proximately 44m tall and 21m by 23m in plan. Figure B.2 shows plan views of the
foundation and a typical floor as well as cross-section views of the foundation and a
N-S cross-section.
The building has concrete moment frames in both the E-W and N-S directions. In
addition, there are shear walls on the East and West sides of the building that provide
most of the stiffness in the N-S direction. Shear walls in the central core provide added
stiffness in both directions. More detailed descriptions of the structural system may
be found in Kuroiwa (1967), Foutch (1976), Luco et al. (1986), Favela (2004), and
Clinton (2004).
Figure B.1: Robert A. Millikan Memorial Library: View from the northeast. The
dark colored walls in the foreground comprise the East shear wall.
213
(c) Floor Plan and Instrumentation of Millikan Library, dense instrument array shown in
red, station MIK (on 9th floor) shown in black. Shaker position (roof level) also shown.
The Millikan Library has been extensively monitored and instrumented since its com-
pletion in 1966 (Kuroiwa, 1967; Trifunac, 1972; Foutch, 1976; Luco et al., 1986;
Chopra, 1995). Clinton (2004) has summarized some of the previous data on Mil-
likan Library’s behavior under forced and ambient vibrations in Appendix B.6. The
evolution of the building behavior, including some dramatic shifts in the fundamental
modes, is documented in Clinton (2004) and is reproduced here in Table B.1 and
Figure B.3. A drop of 21% and 12% for the E-W and N-S fundamental modes since
construction is noted. The primary cause for these shifts appears to be a permanent
loss of structural stiffness which occurs during strong ground motions, most noticeably
the San Fernando (1979) and Whittier Narrows (1987) events. Small fluctuations in
natural frequencies have also been noted which can depend on weather conditions at
the time of testing (Clinton et al., 2006). In particular, the E-W and torsional funda-
mental frequencies have increased by ∼3% in the days following large rainfalls. The
frequencies observed during ambient studies also differ from the frequencies observed
during forced vibration tests (Clinton, 2004).
A1 = 235.73 N · sec2
Fi = Ai f 2 sin (2πf t) A2 = 1518.67 N · sec2 (B.1)
A3 = 3575.89 N · sec2 .
215
Table B.1: History of Millikan Library strong motion behavior. Fundamental Modes.
“%diff1” is the difference between the recorded frequency and that obtained in the
first forced vibration tests (from (Kuroiwa, 1967)). “%diff2” is the difference between
the recorded frequency and that obtained in the most recent forced vibration test prior
to the event. (Adapted from Clinton (2004).)
216
2
34
LC M5.3 !=57
WN M6.1 !=19
SM M5.8 !=18
BB M5.4 !=119
NR M6.7 !=34
BH M4.2 !=26
Frequency, Hz
SS M6.5!=323
351
1.4
Figure B.3: Graphical depiction of Table B.1, the historical behavior of Millikan
Library. Dashed lines represent the E-W natural frequencies, and the dashed-dotted
lines represent the N-S natural frequencies. Shaded area is the likely range of natural
frequencies taking into consideration errors in measurement due to various factors -
weight configuration in the shaker, weather conditions at the time of the test, and
experimental error. Crosses indicate actual time forced test was made. Circles in-
dicate natural frequency estimates from the strong motion record during earthquake
events, and numbers in italics are peak acceleration recorded for the event (cm/s2 ).
[Earthquake Abbreviations: LC: Lytle Creek, SF: San Fernando, WN: Whittier Nar-
rows, SM: Santa Monica, NR: Northridge, BH: Beverly Hills, BB: Big Bear, SS: San
Simeon.] (Adapted from Clinton (2004).)
217
Frequency, f , is in Hz; Ai (a shaker constant) is in N·sec2 ; and the resulting force,
Fi , is in units of N. Table B.2 lists the values of Ai and the limiting frequency for
each weight configuration. For our test we used three shaker levels: A3 , full weights
with the buckets loaded at 100% of capacity; A2 , an intermediate configuration with
two large weights in each of the large weight sections of each bucket, corresponding
to 42.5% of the mass of the full buckets; and A1 , empty buckets, which corresponds
to a shake factor of 6.6% of the full weight configuration.
We can strongly excite the torsional modes through E-W shaking, as the shaker
is located ∼6.1 meters to the South of the building’s N-S line of symmetry (Fig-
ure B.2(c)). The shaker is located ∼0.3 meters to the East of the building’s E-W
centerline, and therefore we do not expect shaking in the N-S direction to effectively
excite the building in torsion.
0 235.73 [9.7] 429.31 [7.2] 622.88 [6.0] 816.45 [5.2] 1010.03 [4.7]
1 877.20 [5.0] 1070.77 [4.6] 1264.35 [4.2] 1457.92 [3.9] 1651.49 [3.7]
2 1518.67 [3.8] 1712.24 [3.6] 1905.81 [3.4] 2099.39 [3.3] 2292.96 [3.1]
3 2160.13 [3.2] 2353.71 [3.1] 2547.28 [3.0] 2740.85 [2.8] 2934.43 [2.8]
4 2801.60 [2.8] 2995.17 [2.7] 3188.75 [2.6] 3382.32 [2.6] 3575.89 [2.5]
Table B.2: Shaker constant, Ai (N ·sec2 ), and limiting frequencies (Hz) for different
configurations of lead weights in the shaker. Bold type indicates the configurations
used in these experiments.
In 1996 the United States Geological Survey (USGS) and the Caltech Department
of Civil Engineering installed a permanent dense array of uni-axial strong-motion
instruments (1g and 2g Kinemetrics FBA-11s) in the Millikan Library, with 36 chan-
nels recording on two 19-bit 18-channel Mt. Whitney digitizers. The instruments are
distributed throughout the building, with three horizontal accelerometers located on
each floor and three vertical instruments in the basement. This dense array is recorded
by the Mt. Whitney digitizer system, providing local hard-drive storage of triggered
events. In 2001 a 3-component Episensor together with a 24-bit Q980 data logger was
installed on the 9th floor. Data from this sensor is continuously telemetered to the
Southern California Seismic Network (SCSN) as station MIK. Figure B.2(c) provides
a schematic of the instrument locations.
A frequency sweep of Millikan Library was performed on July 10, 2002. This test was
designed to identify the natural frequencies and their damping; modeshapes would be
determined with later detailed testing. The building response was recorded using the
SCSN station MIK on the 9th floor and a Mark Products L4C3D seismometer with a
16 bit Reftek recorder on the roof (provided by the Southern California Earthquake
Center (SCEC) Portable Broadband Instrument Center located at the University of
219
California, Santa Barbara). We also used a Ranger seismometer with an oscilloscope
at the roof level to provide an estimate of roof level response during our experiment.
We began with a N-S frequency sweep and the shaker set with empty buckets,
starting near the frequency limit of the shaker at 9.7Hz. We held the frequency con-
stant for approximately 60 seconds to allow the building response to approach steady
state and then lowered the shaker frequency in either .05Hz or .1Hz increments, again
pausing for 60 seconds at each frequency. Once we reached 3.8Hz, we turned off the
shaker and loaded it with two large weights in each of the large weight compartments
in each of the buckets (the intermediate 42.5% loading configuration). We then con-
tinued the frequency sweep from 3.7Hz to 1.5Hz. This procedure was repeated for the
E-W direction, driving the empty shaker from 9.7Hz to 3.8Hz, then sweeping from
3.7Hz to 1Hz using the intermediate configuration.
Figure B.5 shows normalized peak displacement curves for the frequency sweeps.
For each frequency, a representative section from the steady state portion of the data
was selected, bandpass filtered (0.2Hz above and below each frequency, using a 2-pass
3-pole butterworth filter), and fit to a sine wave to estimate the exact frequency, am-
plitude, and phase. These sinusoidal amplitudes were then normalized by the applied
shaker force for the particular frequency and weight combination (Equation B.1).
Furthermore, damping was determined by applying the half-power (band-width)
method (Meirovitch, 1986). We estimated the peak displacement frequency from a
cubic interpolation of the normalized data, as our data sampling is somewhat sparse
for a frequency/amplitude curve. As the higher modes have too much lower mode
participation to determine the half-power points, damping was only determined for
the fundamental modes. Damping is estimated to be 2.28%, 2.39%, and 1.43% for
the E-W, N-S, and Torsional fundamental modes, respectively.
We performed a forced excitation test of Millikan Library on August 18, 2002, record-
ing data using the dense instrumentation network operated by the USGS and station
220
Figure B.5: Lin-Log normalized peak displacement curves for the frequency sweep
performed on July 10, 2002. Amplitudes for E-W and N-S shaking are normalized by
the force factor corresponding to the weight configuration and frequency, as calculated
in Equation B.1. Roof response is shown in solid blue lines, and station MIK (9th
floor) is shown in dashed red lines.
221
MIK. We compare the behavior of the library with the behavior of uniform shear and
bending beams (see Appendix B.5), but it is important to note that these are simple
structural approximations which do not include the behavior of the foundation or the
true structural system of the library.
We began the experiment with the shaker buckets fully loaded and set to excite the
E-W direction. We excited the building at frequencies near the fundamental E-W and
torsional modes, in frequency increments of .03-.05 Hz (again holding for 60s at each
frequency to allow the building to approach steady-state response). With full buckets
we then set the shaker to excite in the N-S direction to examine the fundamental
N-S mode. We then repeated the excitation of the fundamental E-W, N-S, and
torsional modes with the intermediate 42.5% weight configuration, to examine the
shift in natural frequencies with different levels of exciting force. With empty shaker
222
buckets, we excited the first and second E-W overtones, the first torsional overtone,
and the first N-S overtone. Table B.3 summarizes our testing procedure and results.
There are two parallel arrays of instruments in the N-S direction: one set located
on the east side of the library and the other on the west side of the library, as shown
in Figure B.2. In the E-W orientation there is one array, located on the west side of
the building. The two N-S arrays are positioned towards the East and West edges
of the building, far from the E-W centerline, while the E-W array on the west side
of the building is located only 1m from the N-S centerline. Therefore, we expect to
observe torsional response as strong, out of phase motion from the N-S arrays, with
relatively small motion observed from the E-W array.
For each frequency, we selected a representative section from the the steady-state
portion of the data, bandpassed the data (1/2 octave above and below each frequency
using a 2-pass 2-pole butterworth filter), and integrated twice to obtain displacement
values. We created resonance curves by fitting the displacement data to a sine wave
to estimate frequency and amplitude and then normalizing the response based on the
applied force for each frequency and weight combination (Section B.2.2). The mode
shape snapshots in Figures B.6 to B.11 depict the behavior of the building at the
point of maximum roof displacement for each frequency. Using the geometry of the
basement and the position of the vertical basement sensors, we were able to estimate
the rigid body rocking of the building and use it to correct our mode shapes. The
horizontal basement sensors were used to correct for rigid building translation. Our
mode shape figures present the raw results from all three instrument arrays, and the
corrected results with basement translation and rocking removed.
For the torsional modes, Figures B.8 and B.11, we present a snapshot of the
displacement records and also provide a snapshot in terms of rotation angle, θ, at
each floor. The rotation angle at each floor was calculated by subtracting the western
N-S array displacement values (D2 ) from the eastern N-S array displacement values
(D1 ) and dividing by the E-W length (LE-W ) between the arrays (Equation B.2). For
our rotation angle figures, we present the rotation angle at each floor, the basement
223
rotation angle, and the rotation angle corrected for basement rotation.
D1 − D2
θ ≈ tan θ = . (B.2)
LE-W
The instruments in the Eastern N-S array malfunctioned on floors 2 and 8, and
as a result we show a linearly interpolated value for those floors in our mode shape
diagrams.
Figures B.6(a) and B.6(b) show the resonance curve obtained from forced E-W shak-
ing with full weights and 42.5% weights, respectively. Figures B.6(c) and B.6(d)
show the respective mode shapes observed at the resonant frequencies for the differ-
ent weight configurations. Shapes from all three sets of channels are shown on the
same plot – the E-W response clearly dominates during E-W excitation. The observed
mode shapes for different weight configurations are similar, but due to the nonlinear
force-response behavior of the building, the resonant frequency shifts from 1.11Hz
with full weights to 1.14Hz with 42.5% weights. This shift in resonant frequency with
respect to changing the applied force is small, and though obvious, is at the limit of
the resolution of our survey.
Figures B.6(c) and B.6(d) show the mode shapes for both the raw displacements
and the displacements corrected for translation and tilt. The mode shapes have
a strong linear component and closely resemble the theoretical mode shape for a
bending beam, with the inclusion of the kink at the ground floor. See Appendix B.5
for a brief summary and comparison of bending and shear beam behavior. Tilting
and translation effects in this mode account for 3% of the roof displacement.
224
180 180
EW(W) comp
NS(W) comp
160 NS(E) comp 160
Force Normalized Roof Displacement [cm/N]x10−7
120 120
100 100
80 80
60 60
EW(W) comp
40 40 NS(W) comp
NS(E) comp
20 20
0 0
1.05 1.1 1.15 1.2 1.25 1.05 1.1 1.15 1.2 1.25
Frequency, Hz Frequency, Hz
(a) Resonance curve for E-W Shaking (b) Resonance curve for E-W Shaking
with full weights with intermediate (42.5%) weights
R R
9 9
8 8
7 7
6 6
Floor
Floor
5 5
4 4
3 3
Figure B.6: Resonance curves and mode shapes for the E-W fundamental mode
under two loading conditions. Mode shapes are shown corrected (for rigid body
motion) and uncorrected. The mode shapes and resonance curves are shown for the
east-west array located on the west side of the building, EW(W); the western north-
south array, NS(W); and the eastern north-south array, NS(E). Force is calculated as
in Equation B.1 based on the frequency and loading configuration of the shaker.
225
North-South Fundamental Mode
Figures B.7(a) and B.7(d) contain the resonance curves for N-S shaking with full
weights and 42.5% weights, respectively. The fundamental N-S mode is also non-linear
with respect to applied force, and we observe a resonant frequency shift from 1.64Hz
for full weights to 1.67Hz with 42.5% weights. The mode shapes, Figure B.7(c) and
Figure B.7(d) are near identical, and show a more pronounced hinge behavior than the
first E-W mode. When compared to the theoretical mode shapes of Appendix B.5, the
observed shape most closely resembles theoretical bending beam behavior, differing
near the ground floor due to the pronounced hinge behavior in this mode shape. We
also observe that the two N-S arrays are exhibiting in-phase motion and that the E-W
response to N-S shaking is small, as expected. Foundation compliance becomes much
more important for this mode, as we observe that ∼25% of the roof displacement is
due to tilting of the library, and ∼5% is due to translation of the base of the library.
Similar observations for the rigid-body rotation and translation of the building were
made by Foutch et al. (1975).
The fundamental torsional mode involves the twisting of the building and therefore
has more complicated three-dimensional behavior. Due to the positioning of the in-
struments, a small amplitude response is observed from the accelerometers in the
E-W array, while the two N-S arrays recorded a large amplitude out of phase re-
sponse. Figure B.8(a) shows the resonance curve for the fundamental torsional mode.
Figure B.8(b) gives the displacement records for the torsional mode shapes, and Fig-
ure B.8(c) shows the torsional mode shapes in terms of twist angle, θ (as defined in
Equation B.2), instead of displacement. In Figure B.8(b), the two N-S arrays display
the expected out of phase displacements, although some asymmetry is observed.
226
60 60
EW(W) comp EW(W) comp
NS(W) comp NS(W) comp
NS(E) comp NS(E) comp
−7
40 40
30 30
20 20
10 10
0 0
1.6 1.65 1.7 1.75 1.8 1.6 1.65 1.7 1.75 1.8
Frequency, Hz Frequency, Hz
(a) Resonance curve for N-S Shaking (b) Resonance curve for N-S Shaking
with full weights with intermediate (42.5%) weights
R R
9 9
8 8
7 7
6 6
Floor
Floor
5 5
4 4
3 3
Figure B.7: Resonance curves and mode shapes for the N-S fundamental mode
under two loading conditions. Mode shapes are shown corrected (for rigid body
motion) and uncorrected. The mode shapes and resonance curves are shown for the
east-west array, located on the west side of the building, EW(W); the western north-
south array, NS(W); and the eastern north-south array, NS(E). Force is calculated as
in Equation B.1 based on the frequency and loading configuration of the shaker.
R R
25
EW(W) comp
NS(W) comp 9 9
NS(E) comp
8 8
20
7 7
6 6
15
5 5
Floor
Floor
4 4
10
3 3
2 2
5
Figure B.8: Resonance curves and mode shapes for the Torsional fundamental mode. Force is calculated as in Equation B.1
based on the frequency and loading configuration of the shaker.
228
B.4.3 Higher Order Modes
Prior to the installation of the dense instrument array, the higher order mode shapes
were difficult to observe; determining the modeshapes and frequencies for these higher
order modes was one of the primary goals of our suite of experiments.
The first E-W overtone (second E-W mode) has a broad resonance peak, with a
maximum response at 4.93Hz (Figure B.9(a)). The mode shape, seen in Figure B.9(c),
is typical of the second mode shape of a beam in bending (Appendix B.5). The ratio
of the frequency of the second mode to the first mode is 4.32, much lower than the
theoretical ratio for a bending beam of 6.26. For comparison, the theoretical ratio for
a shear beam is 3.
Also observed during our testing was the second E-W overtone (third E-W mode).
Figure B.9(b) shows a resonance peak with a maximum response at 7.83Hz. The mode
shape for this frequency is presented in Figure B.9(d) and is typical of the second
mode of a theoretical shear beam (Appendix B.5). The ratio of the frequency of the
third mode to the second mode is 1.59, lower than the theoretical ratio for a bending
beam of 2.80 and closer to the theoretical ratio for a shear beam of 1.67. The ratio
of third mode to first mode frequencies for a bending beam is 17.55, the ratio for a
shear beam is 5, and for our observed building behavior the ratio is 6.87.
As can be seen in Figure B.10(a), the first N-S overtone (second N-S mode) also has
a broad resonance peak. The resonance curves for the two N-S arrays did not have
their peaks at the same frequency, so this test did not provide a single resonance
peak. However, based on the frequency sweep of Section B.3, and the shapes of the
two resonance curves, we selected 7.22Hz as the modal frequency. The mode shape at
this frequency, shown in Figure B.10(b), is qualitatively typical of a bending beam’s
second mode, but we see that the two N-S arrays have very different amplitudes
229
2 0.5
1.8 0.45
−7
1.2
0.3
1
0.25
0.8
0.2
0.6
0.15
0.4
0.2 0.1
0 0.05
4.8 4.85 4.9 4.95 5 7.7 7.75 7.8 7.85 7.9
Frequency, Hz Frequency, Hz
(a) Second E-W mode. Resonance curve (b) Third E-W mode. Resonance curve for
for E-W shaking with empty buckets E-W shaking with empty buckets
R R
9 9
8 8
7 7
6 6
Floor
Floor
5 5
4 4
3 3
Figure B.9: Second and third E-W modes (first and second E-W overtones), with
resonance curves and mode shapes. Mode shapes are shown corrected (for rigid body
motion) and uncorrected. The mode shapes and resonance curves are shown for the
east-west array, located on the west side of the building, EW(W); the western north-
south array, NS(W); and the eastern north-south array, NS(E). Force is calculated as
in Equation B.1 based on the frequency and loading configuration of the shaker.
230
with zero crossings at different heights. For the N-S second mode, the eastern and
western arrays should have similar shapes and amplitudes (cf. the first N-S mode
in Figure B.7), as the building is approximately symmetric. This implies that we
did not excite the exact modal frequency or that this mode has a more complicated
three-dimensional response than the first N-S mode. The ratio of the frequency for
the second mode (approximate) to the first mode is 4.32, which is close to the ratio
of frequencies observed in E-W bending and is also lower than the theoretical ratio
for the first two modes of a bending beam.
R
0.9
9
0.8
8
Force Normalized Roof Displacement [cm/N]x10−7
0.7
7
0.6
6
0.5
Floor
0.4 4
EW(W) comp
NS(W) comp
0.3 NS(E) comp 3
Figure B.10: Resonance curves and mode shapes for the second NS mode (first NS
overtone). Mode shapes are shown corrected (for rigid body motion) and uncorrected.
The mode shapes and resonance curves are shown for the east-west array, located on
the west side of the building, EW(W); the western north-south array, NS(W); and
the eastern north-south array, NS(E). Force is calculated as in Equation B.1 based
on the frequency and loading configuration of the shaker.
The first torsional overtone was difficult to excite in the building and difficult to
observe. We excited the torsional mode using E-W excitation and expected small
231
torsional response on the E-W channels and large out of phase response from the two
N-S arrays. However, the observed response was dominated by E-W motion from the
E-W shaking used to excite the system, which drove the building in a mode shape
similar to that of the second E-W mode. We observed out of phase motion in the
two N-S arrays, but the response of the N-S arrays was much smaller than the E-W
response. Figure B.11(a) shows the response curve for this mode, which is dominated
by the E-W motion. Figure B.11(b) shows the mode shapes, and Figure B.11(c)
shows the response in terms of twist angle, θ, as defined in Section B.4.1. As with the
N-S overtone, the resonance curve did not clearly identify a modal frequency, but we
chose 6.57Hz as the frequency of interest based on the shapes of the resonance curves
and the results of the frequency sweep of Section B.3.
R R
0.4
9 9
0.35
8 8
0.3
7 7
6 6
0.25
EW(W) comp 5 5
Floor
NS(W) comp
Floor
NS(E) comp
0.2 4 4
3 3
0.15
2 2
Figure B.11: Resonance curves and mode shapes for the second Torsional mode (first Torsional overtone). Force is calculated
as in Equation B.1 based on the frequency and loading configuration of the shaker.
233
B.4.4 Modeshapes Summary
Table B.4 contains a summary of the ratios of frequencies found for Millikan Library,
along with theoretical results for bending and shear beams. Appendix B.5 presents
a summary of theoretical bending and shear beam behavior. To further analyze the
data, the mode shapes were fit using theoretical bending and shear beam behavior by
a modified least squares method. Figure B.12 shows the results of the least squares
curve-fitting for the E-W and N-S modes. The experimental data and best fit are
shown, along with the theoretical mode shapes which are scaled according to their
participation in the best fit curve. Both the fundamental E-W and N-S modes, Fig-
ures B.12(a) and B.12(d), are dominated by the bending component, as are the second
E-W and N-S modes, Figures B.12(b) and B.12(e). The third E-W mode was not
matched well using the third theoretical bending and shear modes; a fit including the
second theoretical bending and shear modes is presented in Figure B.12(c), implying
that the mode shape is best approximated by the second mode of a theoretical shear
beam.
Table B.4: Ratio of frequencies for bending beam behavior, shear beam behavior,
and the observed behavior of Millikan Library.
234
R R R
9 9 9
8 8 8
7 7 7
6 6 6
Floor
Floor
Floor
5 5 5
Least−Square
Fit
4 4 4 Experimental
Data
2nd Bending
3 3 3
Mode
2nd Shear
2 Mode
2 2
3rd Bending
Least−Square Fit Mode
Least−Square Fit
Experimental Data Experimental Data 3rd Shear
1 1 1 Mode
1st Bending Mode 2nd Bending Mode
Translation
1st Shear Mode 2nd Shear Mode
B B B
0 0.2 0.4 0.6 0.8 1 −1 −0.5 0 0.5 −1 −0.5 0 0.5 1
Normalized Displacement Normalized Displacement Normalized Displacement
(a) Least square fit for funda- (b) Least square fit for second (c) Least square fit for third
mental E-W mode. E-W mode (first E-W over- E-W mode (second E-W
tone). overtone). Translational
term included in curve-
fitting to accommodate
large kink at ground level.
R R
9 9
8 8
7 7
6 6
Floor
Floor
5 5
4 4
3 3
2 2
(d) Least square fit for funda- (e) Least square fit for second
mental N-S mode. N-S mode (first N-S over-
tone).
Figure B.12: Least squares curve fitting for E-W and N-S modes. Linear tilt and
translation removed when calculating the best fit.
235
B.5 Theoretical Beam Behavior
The mode shapes (Xn (z)) and associated natural frequencies (ωi ) for a cantilevered
fixed-free bending beam are found by solving the differential equation:
∂ 4 X(z) ω2m
− β 4 X(z) = 0 β4 = , (B.3)
∂z 4 EI
where m is mass per unit length, E is Young’s Modulus and I is the moment of
inertia.
∂X(z)
X(0) = 0 |x=0 = 0 at the fixed end, (B.4)
∂z
(B.5)
∂ 2 X(z) ∂ 3 X(z)
|x=L = 0 |x=L = 0 at the free end. (B.6)
∂z 2 ∂z 3
which can be solved analytically to give the following values for the first three modes:
9 9 9
8 8 8
7 7 7
6 6 6
Floor
Floor
Floor
5 5 5
4 4 4
3 3 3
2 2 2
1 1 1
B B B
−1 0 1 −1 0 1 −1 0 1
X1 X2 X3
Figure B.13: From left to right, theoretical mode shapes for the fundamental mode
(1st mode) and the first two overtones (2nd and 3rd modes) for a cantilevered bending
beam. Mode shapes Xn are normalized such that the maximum displacement is equal
to 1.
(cos βn L + cosh βn L)
Xn (z) = C1 (sin βn z − sinh βn z) + (cos βn z − cosh βn z) .
(sin βn L − sinh βn L)
(B.12)
9 9 9
8 8 8
7 7 7
6 6 6
Floor
Floor
Floor
5 5 5
4 4 4
3 3 3
2 2 2
1 1 1
B B B
−1 0 1 −1 0 1 −1 0 1
X1 X2 X3
Figure B.14: From left to right, theoretical mode shapes for the fundamental mode
(1st mode) and the first two overtones (2nd and 3rd modes) for a cantilevered shear
beam. Mode shapes Xn are normalized such that the maximum displacement is equal
to 1.
Theoretical shear beam behavior is as follows, with the deformed shape being portions
of a sine curve:
(2n − 1)π
Xn (z) = C1 (sin z) n = 1, 2, 3...
2
Tables B.5 and B.6 provide a summary of various studies into the frequencies and
damping of Millikan Library. Table B.7 contains the references used to compile Ta-
bles B.5 and B.6. These studies include ambient and forced vibration testing as well
as data recorded from earthquake ground motions.
239
Table B.5: Summary of Millikan Library Modal Frequency and Damping Analy-
sis Experiments 1967-1994. fo and f1 are the fundamental frequency and the first
overtone in Hz. ζo and ζ1 are the corresponding damping ratios in %. References are
found in Table B.7. A: Ambient, M: Man Excited, F: Forced Vibration, E: Earthquake
Motions [LC: Lytle Creek Earthquake].
240
Table B.6: Summary of Millikan Library Modal Frequency and Damping Analy-
sis Experiments 1987-2003. fo and f1 are the fundamental frequency and the first
overtone in Hz. ζo and ζ1 are the corresponding damping ratios in %. References are
found in Table B.7. A: Ambient, M: Man Excited, F: Forced Vibration, E: Earthquake
Motions [BH: Beverly Hills Earthquake, BB: Big Bear Earthquake].
241
Table B.7: References which correspond to footnote numbers in Tables B.5 and B.6.
242
243
Appendix C
There are many online resources for time-frequency analysis. This appendix contains
some useful collections of tools and examples as well as MATLAB
R functions. (Links
C.1 Links
WAVELAB
http://www-stat.stanford.edu/∼wavelab/
A collection of Matlab functions that contains a Wigner-Ville implementation,
along with a very complete wavelet analysis toolbox. [Good online documentation
for installing the .m files and setting up your path properly, but please note that some
of the .mex files will not compile using modern versions of Matlab, and some routines
will generate error messages when using obsolete commands/syntax.]
http://www.mathworks.com/matlabcentral/link exchange/
User-contributed toolboxes and .m files. Of particular interest is René Laterveer’s
wvd.m package. I have modified this package and include my version below.
http://tftb.nongnu.org/
244
This toolbox contains the most applicable functions for general time-frequency
analysis. The Reduced Interference Distributions in this thesis were created from a
largely unchanged “tfrridh.m” function, which stands for “Time Frequency Response,
Reduced Interference Distribution, Hanning Smoothing Kernel.” Other smoothing
and TFR options are available and there is very thorough online documentation.
C.2 Downloads
Wigner-Ville Distribution
%
% Adapted from Rene Laterveer, 1999, wvd.m package available from:
% http://www.mathworks.com/matlabcentral/link_exchange/
% MATLAB/Signal_processing/
%
245
z=hilbert(x);
% make even number of points, at given resolution
npts = floor(floor(length(z)/res)/2)*2;
% make window an integer
win=floor(win);
% round window length down to nearest odd integer
oddwin = (floor((win-1)/2)*2)+1;
% half point (for indexing reasons we need it later, we’re
% filling two columns per loop, so we only index through half)
halfwin = (oddwin+1)/2-1;
% create tt and ff
tt=[0:npts-1]*res/sps;
ff=[0:(win-1)]*(sps/2)/(win-1);
% pad with zeros
z = [zeros(1,oddwin-1), z, zeros(1,oddwin-1)];
% initialize (important when creating large arrays)
wv = zeros(win,npts);
R = zeros(1, win);
idx = 1:halfwin;
for n=0:npts/2-1
t = 2*n*res+oddwin;
R(1) = z(t)*conj(z(t)) + i*z(t+res)*conj(z(t+res));
v1 = z(t+idx).*conj(z(t-idx));
v2 = z(t+res+idx).*conj(z(t+res-idx));
R(idx+1) = v1+i*v2;
R(win-idx+1) = conj(v1)+i*conj(v2);
RF = fft(R, win);
wv(:,2*n+1) = real(RF);
wv(:,2*n+2) = imag(RF);
end
return;
Double Chirp
dcexample.m – Sample script that will test that the wvdc and doublechirp com-
mands are working correctly. This shows proper syntax for these commands.
[yy,tt]=doublechirp;
[wv,wvfrequency,wvtime]=wvdc(yy,1,length(yy)/2,1000);
figure;
subplot(2,1,1);
imagesc(wvtime,wvfrequency,wv);axis xy;colormap bone
subplot(2,1,2);
surf(wvtime,wvfrequency,wv);axis tight;shading interp;colormap bone
sac2mat.m – You can use this sac2mat.m file to read .sac files into Matlab. The
syntax is: [data,header,station.name]=sac2mat(’1409094992.CI.MIK.HLE.sac’); use
“help sac2mat” for more information. The other sac2mat versions deal with little/big
endian byte-ordering for windows machines and allow you to select your start and
stop points.
function [data,header,kstnm] = sac2mat(sacfile);
% sac2mat loads a SAC2000 file into Matlab, along with
% header information.
% [data,header,kstnm] = sac2mat(sacfile);
%
% sacfile must be a STRING, eg.
% [fff,hhh,name]=sac2mat(’dir1/file1.sac’)
%
% Output:
% data= the data from the file
% header=the first 105 header variables, from the SAC Home Page
% http://www.llnl.gov/sac/
% http://www.llnl.gov/sac/SAC_Manuals/FileFormatPt1.html
% kstnm= KSTNM, the station name (string)
% (header and kstnm are optional)
%
% Some common/useful headers:
% header(1:3) DELTA,DEPMIN,DEPMAX
% header(57) DEPMEN
% header(6:7) B,E
% header(8) Event origin time (seconds relative to reference time)
% header(32:39) STLA,STLO,STEL,STDP,EVLA,EVLO,EVEL,EVDP
% header(51:54) DIST,AZ,BAZ,GCARC
% header(71:76) NZYEAR,NZDAY,NZHOUR,NZMIN,NZSEC,NZMSEC
% header(80) NPTS
% Case Bradford, 8 June 2004
sacfid = fopen(sacfile,’r’);
if sacfid==-1;error=’Not a valid path.’
247
return;end
% The first 70 header variables are floating point
header1 = fread(sacfid,70,’float32’);
% The next 35 are integers (after that they’re a mix of strings)
header2 = fread(sacfid,35,’int32’);
% Outputs the number of points being read to the screen
header=[header1;header2];
npts = header(80)
% kstnm, station name, is a string
% stored in the 110th header variable
% kstnm = KSTNM (First 3 letters of the station name)
fseek(sacfid,4*110,-1);
kstnm = char(fread(sacfid,3)’);
% Now read the data...
fseek(sacfid,158*4,-1);
data(1:npts) = fread(sacfid,npts,’float32’);
status = fclose(sacfid);
return
sac2mata.m – Adjusts the start/stop points of a file. This version is useful when
reading long files.
“help sac2mata” for more information on syntax.
Bibliography
Boashash, B. (1988, September). Note on the Use of the Wigner Distribution for
Time-Frequency Signal Analysis. IEEE Transactions on Acoustics, Speech, and
Signal Processing 36 (9), 1518–1521.
Bracewell, R. N. (2000). The Fourier Transform and its Applications. McGraw Hill
Higher Education.
Clinton, J. F., S. C. Bradford, T. H. Heaton, and J. Favela (2006, February). The ob-
served wander of the natural frequencies in a structure. Bulletin of the Seismological
Society of America 96 (1), 237–257.
Iwan, W. D. (1967). On a class of models for the yielding behavior of continuous and
composite systems. American Society of Mechanical Engineers, Journal of Applied
Mechanics (89), 612–617.
Luco, J., M. Trifunac, and H. Wong (1987). On the apparent change in dynamic
behavior of a 9- story reinforced-concrete building. BULLETIN OF THE SEISMO-
LOGICAL SOCIETY OF AMERICA 77 (6), 1961–1983.
Yang, J., T. H. Heaton, and J. Hall (2006). Simulated Nonlinear Response of High-
rise Buildings for the 2003 Tokachi-oki Earthquake. In 100th Anniversary Earthquake
Conference, April 18-22, 2006, commemorating the 1906 San Francisco Earthquake,
San Francisco, California.