Chapter 3: Time-Domain Signal Processing

1Overview¶

In bridge SHM, sensors continuously collect measurements that vary over time. These time-varying signals contain information about structural behavior, environmental conditions, and potential damage. Time-domain signal processing analyzes how a signal evolves over time without transforming it to other domains and is often the first step when examining raw sensor data. This chapter covers techniques for preprocessing and interpreting time-series data (e.g., acceleration or strain measurements) directly in the time domain. These time-domain fundamentals form the foundation for the more advanced analysis methods covered in subsequent chapters.

2What is a Time Series?¶

A time series is a sequence of data points indexed in time order, typically collected at equally spaced intervals. For SHM applications, common time series include:

Acceleration data: Measured in m/s² or g-units from accelerometers
Strain data: Measured in microstrain ( $\mu\varepsilon$ ) from strain gauges
Displacement data: Measured in mm or cm from displacement sensors
Temperature data: Environmental measurements affecting structural response

In mathematical terms, a digital (or discrete-time) signal is a sequence of numbers like $\{x[0], x[1], \ldots, x[N-1]\}$ of length $N$ obtained by sampling a continuous-time signal $x(t)$ at equally spaced time instants. If the sampling period is $T_s > 0$ seconds, then the sampling frequency is $F_s = \frac{1}{T_s}$ hertz, and the $n$ -th sample at time instant $t_n = nT_s$ is

x[n]\;\triangleq\;x(nT_s), \quad n = 0,1,\ldots, N-1.

(1)

<Figure size 1400x1200 with 3 Axes> — Signal sampling from continuous to discrete time series.

2.1Pulse amplitude modulation¶

The sampling operation can be represented in continuous time using an impulse (Dirac-comb) model. The sampled digital signal is

x_s(t)\;=\;\sum_{n=0}^{N-1} x[n]\,\delta(t-nT_s),

(2)

This idealizes pulse-amplitude modulation by infinitely short pulses whose areas equal the samples.

<Figure size 1000x600 with 2 Axes> — Figure 2:ontinuous-time signal $x(t)$ and the sampling pulse train with spacing $T_s$ .

import numpy as np, pandas as pd, matplotlib.pyplot as plt, matplotlib.dates as mdates

# Okabe–Ito (colorblind-safe) palette
COL = {
    "black":      "#000000",
    "orange":     "#E69F00",
    "sky":        "#56B4E9",
    "green":      "#009E73",
    "yellow":     "#F0E442",
    "blue":       "#0072B2",
    "vermillion": "#D55E00",
    "purple":     "#CC79A7",
    "grey":       "#7F7F7F",
}

plt.rcParams.update({
    "font.family": "sans-serif",
    "font.sans-serif": ["Arial", "Liberation Sans", "Nimbus Sans", "Nimbus Sans L", "Helvetica", "DejaVu Sans"],
    "mathtext.fontset": "dejavusans",
    "figure.dpi": 150, "savefig.dpi": 300,
    "axes.linewidth": 1.0,
    "axes.titlesize": 18, "axes.labelsize": 16,
    "xtick.labelsize": 13, "ytick.labelsize": 13, "legend.fontsize": 12,
})

def style_axes(ax, grid=True):
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    if grid:
        ax.grid(True, which="major", linewidth=0.6, alpha=0.25)

3Basic principles of sampling¶

Discrete-time analysis uses two frequency units in parallel: analog frequency $f$ (Hz) with $\Omega=2\pi f$ (rad/s), and digital (normalized) frequency $k=f/F_s$ (cycles/sample) with $\Omega_k=2\pi k$ (rad/sample). Sampling the analog sinusoid $x(t)=A\sin(2\pi f t+\theta)$ at rate $F_s$ gives the digital sinusoid

x[n]=A\sin(2\pi k\,n+\theta), \qquad k=\frac{f}{F_s}.

(3)

Because $\sin(2\pi r +\theta)=\sin(\theta)$ for any integer $r=k \times n$ , digital sinusoids whose frequencies differ by an integer multiple of $2\pi$ (equivalently, whose $k$ differ by an integer) are identical. A unique digital frequency can therefore be restricted to the unit band

-\pi\le \Omega_k < \pi \quad\Longleftrightarrow\quad -\tfrac12 \le k < \tfrac12 .

(4)

If the analog signal is bandlimited to $f_{\max}$ , ideal alias-free sampling requires the sampling theorem

F_s > 2 f_{\max}.

(5)

The boundary $F_s/2$ is the Nyquist frequency; sampling above this rate prevents the spectral replicas generated by sampling from overlapping.

3.1Aliasing¶

If (5) is violated, different analog sinusoids collapse to the same digital sinusoid: aliasing.

<Figure size 1650x900 with 2 Axes> — Figure 3:Digital sinusoid identity.

Practical Sampling Guidelines for SHM

In SHM, aliasing can cause a high-frequency mode or event to appear as a lower-frequency phenomenon. Good practice requires:

Choose $f_s$ well above the highest expected significant frequency. Often a factor of 5–10 above the highest mode of interest gives safe margin.
Use anti-alias filters. Most DAQ hardware includes analog filters just before the ADC that cut off frequencies above $f_s/2$ . For instance, if sampling at 100 Hz, a low-pass filter with cutoff 45 Hz (with a certain roll-off rate) will remove most content above 50 Hz.
If downsampling (resampling to a lower $f_s$ ) is needed to reduce data size or focus on low frequencies, always low-pass filter the data first. For example, to decimate from 200 Hz to 50 Hz, apply a low-pass filter at 25 Hz before picking every 4th sample.

4Components of Time-Series Data¶

A time-series can be decomposed into underlying components: trend, seasonality, cyclic variations, and irregular (noise) components. Identifying these components aids preprocessing (e.g., removing trends or seasonal effects) and correct data interpretation. Consider a single foil strain gauge mounted along the girder axis at midspan on the bottom flange of a simply supported steel bridge. Under positive bending from traffic, this location is typically in tension. We express strain in microstrain ( $\mu\varepsilon$ ). For context, a common serviceability target is to keep peak service strains near or below about 0.12% ≈ 1200 $\mu\varepsilon$ ; “yield-level” strain depends on the steel grade: $\varepsilon_y \approx f_y/E$ with $E \approx 200$ GPa, so for $f_y=250\text{–}350$ MPa, $\varepsilon_y \approx 0.125\%\text{–}0.175\%$ (1250–1750 $\mu\varepsilon$ ). Thus, a generic 0.20% (2000 $\mu\varepsilon$ ) line is above yield for many steels and should be treated as an illustrative “do-not-exceed” proxy rather than a universal safety limit.

4.1Additive model¶

The measured strain $X(t)$ is commonly treated as an additive sum of four components,

X(t) \;=\; T(t) \;+\; S(t) \;+\; C(t) \;+\; R(t),

(6)

where $T$ is trend, $S$ is seasonal (regular, fixed-period patterns), $C$ is cyclic (repetitive but not strictly periodic), and $R$ is the residual. This additive model is appropriate because superposition is a good approximation for strain in linear elastic service.

4.1.1Trend, $T(t)$ ¶

The trend is the slow, long-term movement of the baseline strain under comparable environmental and operational conditions. In a bridge application it may creep upward over months if average traffic grows, if stiffness is slowly lost (e.g., fatigue cracking or corrosion reducing section), or if supports settle. The trend evolves on the scale of weeks to years and is usually estimated by a smooth filter (e.g., long moving average or LOESS). Detrending prevents gradual, expected drift from confounding anomaly detection.

4.1.2Seasonality, $S(t)$ ¶

Seasonality repeats with a known period and must not be mistaken for damage when its amplitude is within expected limits. Two drivers dominate for a midspan bottom-flange gauge:

Temperature: apparent/mechanical thermal strain follows $\varepsilon_{\text{th}} \approx \alpha\,\Delta T$ with $\alpha \approx 11\text{–}13$ $\mu\varepsilon$ /°C for steel. A day–night swing of 30 °C therefore produces on the order of $30\times 12 \approx 360$ $\mu\varepsilon$ change, large enough to shift the baseline strain but still normal if combined peaks remain below service limits. Uniform temperature change mainly produces free expansion (little stress in an unconstrained member) but the gauge still reports apparent strain; through-depth thermal gradients add curvature and thus real bending strain at the flange.
Regular usage cycles: diurnal rush-hour loading and weekday/weekend differences repeat predictably (periods of 24 h and 7 d). For anomaly detection, these patterns should be modeled and removed (or conditioned on: |T, S$) so the SHM system does not alarm on normal daily or seasonal swings. For fatigue assessment, the magnitude and count of these seasonal cycles are directly relevant and should be retained for damage accumulation calculations (e.g., rainflow + S-N).

4.1.3Cyclic, $C(t)$ ¶

Cyclic variations are repeatable but not locked to a fixed period. Think of school terms, harvest seasons, tourist waves, or medium-term detours that temporarily reroute trucks. Their amplitude and timing wander; forcing them into the seasonal template will distort $S(t)$ and inflate residuals. Instead, capture them as a separate slow, quasi-periodic component (often via a broader-band smoother than $T$ ).

4.1.4Residual, $R(t)$ ¶

After explaining $T$ , $S$ , and $C$ , the remainder contains high-frequency noise, unmodeled short events (single heavy-truck bursts, braking near the gauge), and truly unusual behavior. This is where sudden, non-recurring changes appear, so robust outlier detection and event logic are applied to $R(t)$ . In contrast to trend shifts, residual spikes are brief and often sparse. Always guard against false alarms from sensor faults or temperature-induced zero-shifts by cross-checking with temperature channels and neighboring sensors.

The following figure shows the decomposition of midspan bottom-flange strain time-series data from a simply supported bridge girder into its underlying components.

<Figure size 2400x2000 with 1 Axes> — Figure 4:Decomposition of measured strain at the midspan bottom flange of a simply supported bridge girder.

4.2Multiplicative models¶

In some time-series, the amplitude of periodic variation grows or shrinks with the overall level. In those cases a multiplicative model gives a cleaner separation of structure and noise than an additive model:

X(t)=T(t)\times S(t)\times C(t)\times R(t)

(7)

When to prefer multiplicative over additive:

Heteroscedasticity is present (variance rises with the mean)
Effects are naturally proportional (e.g., “+15% summer uplift”)
Seasonality is stable in ratios after dividing out the level.

For midspan bottom-flange strain, the additive model is often adequate because elastic strains superpose linearly. The multiplicative view is useful when analyzing relative strain (percent of a baseline) or when operational influences scale the amplitude (e.g., warmer seasons and heavier traffic both raise the baseline and the size of daily swings). The multiplicative model is log-additive:

\log X(t)=\log T(t)+\log S(t)+\log C(t)+\log R(t),

(8)

The same tools used for additive decomposition (moving averages, LOESS, STL, harmonic regression) apply in log space; exponentiate the components afterward (optionally with a small log-normal bias correction) for interpretation and alarm setting.

4.3Implementation of the methods¶

Decomposition proceeds as follows: estimate a smooth trend, remove it, extract the periodic component over the relevant cycle, and isolate the residual. The steps in Algorithm 1 specify when to choose additive vs. multiplicative models and how to validate the residual as noise.

Algorithm 1 (Time-series decomposition (additive vs. multiplicative))

Input: Discrete series $X(t_i)$ ; sampling interval $\Delta t$ ; known/assumed period(s) $\mathcal{P}$ (e.g., daily, weekly, yearly); chosen model type (“additive” or “multiplicative”).

Output: Trend $T(t)$ , seasonal $S(t)$ , (optional) cyclic $C(t)$ , residual $R(t)$ .

Model choice & prechecks.
- Inspect amplitude–level coupling (e.g., plot local sd vs. local mean, try Box–Cox).
- If variance grows with the level → prefer multiplicative; otherwise additive.
- Handle missing data; for multiplicative ensure $X(t)>0$ (use tension channel, $|X|$ , or add a small positive offset).
Trend estimation. Fit a smooth $T(t)$ with a window $\gg$ the dominant period(s) (e.g., LOESS, moving average, spline, or regression with calendar/temperature covariates).
Detrend.
- Additive: $X_{\mathrm{det}}(t)=X(t)-T(t)$ .
- Multiplicative: $X_{\mathrm{det}}(t)=\dfrac{X(t)}{T(t)}$ (or work in logs: $\log X-\log T$ ).
Seasonal extraction. Estimate $S(t)$ over each period in $\mathcal{P}$ (template averaging, harmonic regression, or STL).
- Normalize: additive $S$ has mean $\approx 0$ ; multiplicative $S$ has mean $\approx 1$ .
(Optional) Cyclic component. Model slower, non-strictly periodic variation as $C(t)$ via broader smoothing of $X_{\mathrm{det}}/S$ (multiplicative) or $X_{\mathrm{det}}-S$ (additive).
Residual.
- Additive: $R(t)=X(t)-T(t)-S(t)$ (minus $C$ if used).
- Multiplicative: $R(t)=\dfrac{X(t)}{T(t)\,S(t)}$ (divide also by $C$ if used). In log-space: $\log R=\log X-\log T-\log S\,[ -\log C]$ .
Validation. Check $R(t)$ is “noise-like”: near-zero mean (additive) or near-one mean (multiplicative), reduced variance, weak autocorrelation (e.g., Ljung–Box), and no remaining seasonality. For multiplicative back-transform, apply log-normal bias correction if needed ( $\times \exp(\tfrac{1}{2}\sigma^2_{\log R})$ ).

<Figure size 1650x1575 with 1 Axes> — Figure 6:Additive decomposition of midspan bottom-flange strain (no cyclic): original X(t), detrended X(t)-T(t), and residual R(t).

# Global plotting style (Nature-like) + Okabe–Ito palette
import numpy as np, pandas as pd, matplotlib.pyplot as plt, matplotlib.dates as mdates

COL = {
    "black":      "#000000",
    "orange":     "#E69F00",
    "sky":        "#56B4E9",
    "green":      "#009E73",
    "yellow":     "#F0E442",
    "blue":       "#0072B2",
    "vermillion": "#D55E00",
    "purple":     "#CC79A7",
    "grey":       "#7F7F7F",
}

plt.rcParams.update({
    "font.family": "sans-serif",
    "font.sans-serif": ["Arial", "Liberation Sans", "Nimbus Sans", "Nimbus Sans L", "Helvetica", "DejaVu Sans"],
    "mathtext.fontset": "dejavusans",
    "figure.dpi": 150,
    "savefig.dpi": 300,
    "axes.linewidth": 1.0,
    "axes.titlesize": 18,
    "axes.labelsize": 16,
    "xtick.labelsize": 13,
    "ytick.labelsize": 13,
    "legend.fontsize": 12,
})

def style_time_axes(ax):
    """Nature-like axes, concise time ticks."""
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    ax.grid(True, which="major", linewidth=0.6, alpha=0.25)
    locator = mdates.AutoDateLocator(minticks=4, maxticks=8)
    ax.xaxis.set_major_locator(locator)
    ax.xaxis.set_major_formatter(mdates.ConciseDateFormatter(locator))
    for lab in ax.get_xticklabels():
        lab.set_rotation(0)

5Detrending¶

Detrending removes a known or modeled trend from a time-series so that analysis can focus on the fluctuations of interest. In bridge SHM, trends may arise from environmental effects or long-term structural changes. When tracking slow deterioration the trend itself is of interest; otherwise, separating the slowly-varying component from the dynamic response is the goal.

5.1Subtracting the mean¶

The simplest detrend is removing the constant offset (mean). For many vibration signals this is done by default (AC coupling), since the analysis focuses on oscillations around zero. Subtracting the mean produces a zero-average time-series, which many downstream analyses require.

# Data generation for mean detrending example
import numpy as np, pandas as pd, os

# Create data directory if it doesn't exist
data_dir = "chapters/ch03/ch03_data"
os.makedirs(data_dir, exist_ok=True)

# Define file paths
mean_data_file = os.path.join(data_dir, "mean_detrend_data.npy")
mean_time_file = os.path.join(data_dir, "mean_detrend_time.npy")

# Check if data already exists
if not os.path.exists(mean_data_file):
    # 48 hours at 5-minute sampling
    t = pd.date_range("2024-07-01 00:00", periods=48*12, freq="5min")
    n = len(t)
    rng = np.random.default_rng(1)

    baseline = 420.0  # $\mu\varepsilon$
    thermal = 150.0*np.sin(2*np.pi*(t.hour + t.minute/60)/24)  # daily thermal
    traffic  =  40.0*np.sin(2*np.pi*(t.hour-8)/12)             # rush-hour-ish
    noise    = np.random.normal(0, 6.0, n)

    X = baseline + thermal + traffic + noise  # original strain ($\mu\varepsilon$)
    
    # Save the generated data
    np.save(mean_data_file, X)
    np.save(mean_time_file, t)

5.2Linear detrend¶

When the data has an approximately linear trend, fit a straight line (least squares) to the data and subtract it. This leaves the residual fluctuations around that line. For example, bridge strain data over a day may show a roughly linear increase due to a temperature rise from morning to afternoon; removing that linear component highlights other deviations.

# Data generation for linear detrending example
import numpy as np, pandas as pd, os

# Define file paths
linear_data_file = os.path.join(data_dir, "linear_detrend_data.npy")
linear_time_file = os.path.join(data_dir, "linear_detrend_time.npy")

# Check if data already exists
if not os.path.exists(linear_data_file):
    t = pd.date_range("2024-07-02 00:00", periods=48*12, freq="5min")
    n = len(t)
    rng = np.random.default_rng(2)

    baseline = 400.0
    linear   = np.linspace(0, 60.0, n)                  # ~+60 $\mu\varepsilon$ over 48 h
    thermal  = 120.0*np.sin(2*np.pi*(t.hour+ t.minute/60)/24)
    noise    = rng.normal(0, 7.0, n)
    X = baseline + linear + thermal + noise
    
    # Save the generated data
    np.save(linear_data_file, X)
    np.save(linear_time_file, t)

5.3Polynomial detrend¶

When the trend is nonlinear (e.g., a quadratic temperature curve over a day, or a logarithmic creep curve), fit higher-order polynomials or other nonlinear functions and subtract them to obtain stationary residuals.

# Data generation for polynomial detrending example
import numpy as np, pandas as pd, os

# Define file paths
poly_data_file = os.path.join(data_dir, "poly_detrend_data.npy")
poly_time_file = os.path.join(data_dir, "poly_detrend_time.npy")

# Check if data already exists
if not os.path.exists(poly_data_file):
    t = pd.date_range("2024-07-03 00:00", periods=48*12, freq="5min")
    n = len(t)
    rng = np.random.default_rng(3)

    baseline = 405.0
    # Nonlinear drift (quadratic hump) + diurnal + noise
    u = np.linspace(-1, 1, n)
    nonlinear = -70.0*(u**2) + 40.0                         # inverted parabola
    thermal   = 110.0*np.sin(2*np.pi*(t.hour + t.minute/60)/24)
    noise     = rng.normal(0, 7.0, n)
    X = baseline + nonlinear + thermal + noise
    
    # Save the generated data
    np.save(poly_data_file, X)
    np.save(poly_time_file, t)

5.4Moving average (or moving median) baseline removal¶

A moving average with a large window captures slow changes (the trend); subtracting it from the signal acts as a high-pass filter. A moving median can be used similarly to reduce the influence of outliers.

# Data generation for moving average detrending example
import numpy as np, pandas as pd, os

# Define file paths
ma_data_file = os.path.join(data_dir, "ma_detrend_data.npy")
ma_time_file = os.path.join(data_dir, "ma_detrend_time.npy")

# Check if data already exists
if not os.path.exists(ma_data_file):
    t = pd.date_range("2024-07-04 00:00", periods=48*12, freq="5min")
    n = len(t)
    rng = np.random.default_rng(4)

    baseline = 410.0
    slow_drift = 50.0*np.sin(2*np.pi*np.linspace(0, 1, n))  # very slow swing across 48h
    thermal    = 130.0*np.sin(2*np.pi*(t.hour + t.minute/60)/24)
    noise      = rng.normal(0, 6.0, n)
    X = baseline + slow_drift + thermal + noise
    
    # Save the generated data
    np.save(ma_data_file, X)
    np.save(ma_time_file, t)

5.5Differencing¶

Taking the first difference ( $y[n] = x[n] - x[n-1]$ ) removes a constant trend (it becomes zero) or converts a linear trend to a constant. In SHM, differencing acceleration yields jerk, which is not typical, but for slowly varying strains differencing can highlight incremental changes.

# Data generation for differencing detrending example
import numpy as np, pandas as pd, os

# Define file paths
diff_data_file = os.path.join(data_dir, "diff_detrend_data.npy")
diff_time_file = os.path.join(data_dir, "diff_detrend_time.npy")

# Check if data already exists
if not os.path.exists(diff_data_file):
    t = pd.date_range("2024-07-05 00:00", periods=48*12, freq="5min")
    n = len(t)
    rng = np.random.default_rng(5)

    baseline = 400.0
    linear   = np.linspace(0, 40.0, n)                   # slow rise
    thermal  = 100.0*np.sin(2*np.pi*(t.hour + t.minute/60)/24)
    noise    = rng.normal(0, 6.0, n)
    X = baseline + linear + thermal + noise
    
    # Save the generated data
    np.save(diff_data_file, X)
    np.save(diff_time_file, t)

6Statistical Descriptors in the Time Domain¶

Statistical descriptors quantify the properties of time-domain signals and serve as damage-sensitive features that can indicate changes in structural condition. They capture characteristics of the signal’s distribution, energy content, and extreme values.

6.1Mean (Average)¶

The mean $\mu_x$ of a signal $x[n]$ over a time window is defined as the arithmetic average of all values:

\mu_x = \mathbb{E}\{X\} = \frac{1}{N} \sum_{n=0}^{N-1} x[n]

(9)

where $x[n]$ are the sampled values of the random variable $X$ and $N$ is the number of samples. In acceleration signals the mean is often close to zero after removing any DC bias. For strain signals the mean represents a baseline due to static loads. Removing the mean (zero-centering) prevents a non-zero mean from distorting other statistics such as variance.

6.2Variance and standard deviation¶

The variance $\sigma_x^2$ measures the spread of the data around the mean:

\sigma_x^2 = \mathbb{E}\{(X - \mu_x)^2\} = \frac{1}{N-1} \sum_{n=0}^{N-1} (x[n] - \mu_x)^2

(10)

The standard deviation $\sigma_x$ is simply the square root of variance:

\sigma_x = \sqrt{\sigma_x^2}

(11)

These metrics quantify the energy or variability in the signal.

6.3Root Mean Square (RMS)¶

RMS quantifies the signal’s effective magnitude and is widely used in vibration analysis. It is defined as:

x_{\mathrm{RMS}} = \sqrt{\frac{1}{N} \sum_{n=0}^{N-1} x[n]^2} = \sqrt{\mathbb{E}\{X^2\}}

(12)

The relationship between RMS, mean, and standard deviation is:

x_{\mathrm{RMS}}^2=\mu_x^2+\sigma_x^2.

(13)

6.4Crest Factor¶

The crest factor is the ratio of the peak amplitude (absolute peak) to the RMS value:

\text{CF} = \frac{\max |x[n]|}{x_{\mathrm{RMS}}}.

(14)

This dimensionless number indicates how “peaky” the signal is. A pure sinusoid has $\text{CF} \approx 1.414$ (peak = 1, RMS = $1/\sqrt{2}$ ). Occasional spikes drive the crest factor high because the peak greatly exceeds the RMS. The crest factor combines peak and energy information and is less sensitive to absolute scaling than the peak alone.

6.5Skewness (Third Moment)¶

Skewness measures the asymmetry of the signal value distribution. It is defined as:

S = \mathbb{E}\{(\frac{X - \mu_x}{\sigma_x})^3\} = \frac{1}{N} \sum_{n=0}^{N-1} (\frac{x[n] - \mu_x}{\sigma_x})^3.

(15)

$S = 0$ : Symmetric distribution
$S > 0$ : Right-skewed (tail extends toward positive values)
$S < 0$ : Left-skewed (tail extends toward negative values)

When $S=0$ the distribution is symmetric around the mean. A positive skew indicates a longer right tail; a negative skew indicates a longer left tail.

<Figure size 2250x750 with 3 Axes> — Figure 12:Visual illustration of skewness in probability distributions.

6.6Kurtosis (Fourth Moment)¶

Kurtosis $K$ measures the “tailedness” of the distribution, i.e., how heavy or light the tails are relative to a normal distribution. One definition (for excess kurtosis) is:

K = \mathbb{E}\{(\frac{X - \mu_x}{\sigma_x})^4\}- 3 = \frac{1}{N} \sum_{n=0}^{N-1} (\frac{x[n] - \mu_x}{\sigma_x})^4 - 3,

(16)

Subtracting 3 makes $K=0$ for a Gaussian (the Pearson definition does not subtract 3, giving a normal distribution kurtosis of 3). High kurtosis indicates frequent extreme deviations: a sharp peak and heavy tails. Low kurtosis indicates a flatter distribution with light tails. Because a single large spike can substantially raise kurtosis, tracking it over time can detect the onset of irregular events.

$K = 0$ : Normal distribution
$K > 0$ : Heavy tails (more extreme values)
$K < 0$ : Light tails (fewer extreme values)

6.7Implementation example¶

The following example examines a realistic acceleration measurement from a highway bridge under ambient traffic loading.

import numpy as np, matplotlib.pyplot as plt

# Okabe–Ito palette and Nature-like typography
COL = {"black":"#000000","orange":"#E69F00","sky":"#56B4E9","green":"#009E73",
       "yellow":"#F0E442","blue":"#0072B2","vermillion":"#D55E00","purple":"#CC79A7","grey":"#7F7F7F"}
plt.rcParams.update({"font.family":"sans-serif","font.sans-serif":["Arial","Liberation Sans","Nimbus Sans","Nimbus Sans L","Helvetica","DejaVu Sans"],"mathtext.fontset":"dejavusans","figure.dpi":150,"savefig.dpi":300,
                     "axes.linewidth":1.0,"axes.titlesize":18,"axes.labelsize":16,
                     "xtick.labelsize":13,"ytick.labelsize":13,"legend.fontsize":12})

def ar2_modal(Fs, fn, zeta, N, sigma=1.0, rng=None):
    """Stationary narrowband via AR(2) with poles r·e^{±jω0} (exact pole mapping from ζ, fn)."""
    if rng is None: rng = np.random.default_rng()
    dt = 1/Fs
    w_n = 2*np.pi*fn
    w_d = w_n*np.sqrt(max(0.0, 1.0-zeta*zeta))
    r   = np.exp(-zeta*w_n*dt)
    w0  = w_d*dt
    a1, a2 = 2*r*np.cos(w0), -r*r                  # y[n] = a1 y[n-1] + a2 y[n-2] + b e[n]
    
    # Create strongly non-Gaussian excitation with very heavy tails and asymmetry
    # Moderate, left-skewed excitation using skew-normal (aim skew ≈ -1)
    try:
        from scipy.stats import skewnorm
        e_skew = skewnorm.rvs(a=-6.0, size=N)
    except Exception:
        # Fallback: Gaussian plus small negative lognormal for skew
        e_skew = rng.normal(0.0, 1.0, N) - 0.6*rng.lognormal(0.0, 0.6, N)

    # Small Gaussian baseline
    e_gaussian = rng.normal(0.0, 0.3, N)

    # Occasional small transients with negative bias (kept mild)
    event_mask = rng.random(N) < 0.004  # 0.4% events
    e_shot = event_mask * (rng.exponential(2.0, N)) * rng.choice([1, -1], size=N, p=[0.3, 0.7])

    # Combine to target moderate skew and limited tails
    e = 0.8*e_skew + 0.3*e_gaussian + 0.2*e_shot
    
    y  = np.zeros(N); b = sigma
    for n in range(2, N): y[n] = a1*y[n-1] + a2*y[n-2] + b*e[n]
    return y

def generate_realistic_noise(N, rng, base_std=0.004):
    """Generate strongly non-Gaussian noise with intermittent events typical of real structures."""
    # Base noise: skew-normal with mild negative skew (target γ1 ≈ -1)
    try:
        from scipy.stats import skewnorm
        base = skewnorm.rvs(a=-6.0, size=N)
    except Exception:
        base = rng.normal(0.0, 1.0, N) - 0.6*rng.lognormal(0.0, 0.6, N)
    base_noise = base * (base_std * 6.0)

    # Intermittent micro-events (skewed lognormal magnitudes)
    event_prob = 0.004
    events = rng.random(N) < event_prob
    event_magnitudes = rng.lognormal(0.0, 0.9, size=N) * rng.choice([1, -1], size=N, p=[0.35, 0.65])
    intermittent = events * event_magnitudes * base_std * 6

    # Bursty clusters using short Poisson windows
    burst_onsets = np.where(rng.random(N) < 0.001)[0]
    bursts = np.zeros(N)
    for idx in burst_onsets:
        L = int(rng.integers(10, 40))
        end = min(N, idx + L)
        amp = rng.pareto(2.2) * base_std * 10 * rng.choice([1, -1], p=[0.4, 0.6])
        envelope = np.exp(-np.linspace(0, 3, end - idx))
        bursts[idx:end] += amp * envelope

    # Low-frequency drift with sawtooth-like asymmetry
    drift_freq = 0.01
    t = np.arange(N) / 200.0  # Assuming Fs=200
    saw = (2*(t*drift_freq - np.floor(0.5 + t*drift_freq)))
    drift = -base_std * 0.6 * np.sign(saw) * (np.abs(saw) ** 0.5)  # mild negative bias

    return base_noise + intermittent + bursts + drift

# Data management and synthesis
import os
Fs, T = 200.0, 180.0
N     = int(Fs*T)
t     = np.arange(N)/Fs

# Create data directory if it doesn't exist #
data_dir = "chapters/ch03/ch03_data"
os.makedirs(data_dir, exist_ok=True)

# Define file paths (versioned)
acc_data_file = os.path.join(data_dir, "acceleration_data_v2.npy")
time_data_file = os.path.join(data_dir, "time_vector_v2.npy")

# Check if data already exists
if not os.path.exists(acc_data_file):    
    rng = np.random.default_rng(23)
    
    # Generate modal responses with non-Gaussian excitation
    y1 = ar2_modal(Fs, fn=2.5, zeta=0.025, N=N, sigma=0.15, rng=rng)
    y2 = ar2_modal(Fs, fn=7.5, zeta=0.035, N=N, sigma=0.10, rng=rng)
    
    # Add realistic non-Gaussian noise with intermittent events
    noise = generate_realistic_noise(N, rng, base_std=0.004)
    # Apply light smoothing to simulate sensor response
    noise = np.convolve(noise, np.ones(3)/3, mode="same")
    
    # Combine modal responses and noise
    x = y1 + y2 + noise
    
    # Add more frequent and varied transient events for strong non-Gaussianity
    # Large transient events (door slams, footsteps, impacts)
    transient_times = rng.choice(N, size=int(N*0.001), replace=False)  # Increased to 0.1% of samples
    for idx in transient_times:
        if idx < N-100:  # Ensure we don't go out of bounds
            # Damped exponential transient with random characteristics
            decay_samples = min(rng.integers(20, 100), N-idx)  # Variable decay length
            decay_rate = rng.uniform(0.05, 0.3)  # Variable decay rate
            decay = np.exp(-np.arange(decay_samples) / (Fs * decay_rate))
            # Use log-normal for strongly skewed amplitudes
            amplitude = rng.lognormal(0, 0.8) * 0.03 * rng.choice([1, -1], p=[0.7, 0.3])
            x[idx:idx+decay_samples] += amplitude * decay
    
    # Add micro-transients (small clicks, electrical noise spikes)
    micro_transient_times = rng.choice(N, size=int(N*0.005), replace=False)  # 0.5% of samples
    for idx in micro_transient_times:
        if idx < N-5:
            # Short sharp spikes
            spike_length = rng.integers(1, 6)
            spike_shape = np.exp(-np.arange(spike_length) / 2)
            amplitude = rng.pareto(2.0) * 0.005 * rng.choice([1, -1], p=[0.6, 0.4])
            end_idx = min(idx + spike_length, N)
            x[idx:end_idx] += amplitude * spike_shape[:end_idx-idx]
    
    # Calibrate distribution skewness toward target value (moderate left skew)
    # We use a monotonic quadratic transform x -> z + b*(z^2 - 1) on standardized data
    def _compute_skew(u):
        u = u - np.mean(u)
        m2 = np.mean(u*u)
        if m2 <= 0: return 0.0
        m3 = np.mean(u*u*u)
        return m3 / (m2**1.5)
    
    target_skew = -10.0
    tol = 0.05
    z = (x - np.mean(x)) / (np.std(x) + 1e-12)
    lo, hi = -2.0, 2.0
    best_b, best_err = 0.0, 1e9
    for _ in range(18):
        b = 0.5*(lo+hi)
        zt = z + b*(z*z - 1.0)
        s = _compute_skew(zt)
        err = abs(s - target_skew)
        if err < best_err:
            best_b, best_err = b, err
        if s > target_skew:
            # too positive -> push more negative
            lo = b
        else:
            hi = b
        if err < tol:
            break
    z = z + best_b*(z*z - 1.0)
    # de-standardize using original RMS after normalization below
    x = z
    
    # Target RMS
    target_rms = 0.056
    x = x - np.mean(x)
    x *= target_rms / (np.sqrt(np.mean(x**2)) + 1e-12)
    
    # Save the generated data
    np.save(acc_data_file, x)
    np.save(time_data_file, t)

<Figure size 1650x450 with 1 Axes> — Figure 14:Loading and plotting a sample acceleration data with non-Gaussian characteristics.

Now let’s analyze this acceleration time series as a random variable and examine its probability distribution.

<Figure size 1200x900 with 1 Axes> — Figure 15:Statistical distribution of the reference singal.

import pandas as pd, numpy as np

x_series = pd.Series(x)
summary = pd.DataFrame({
    "μ̂ = mean(x)":           [x_series.mean()],
    "s² = var(x) (N−1)":      [x_series.var(ddof=1)],
    "s = std(x) (N−1)":       [x_series.std(ddof=1)],
    "x_RMS":                  [np.sqrt(np.mean(x**2))],
    "Peak |x|":               [np.max(np.abs(x))],
    "CF = Peak/RMS":          [np.max(np.abs(x))/np.sqrt(np.mean(x**2))],
    "Skewness (sample)":      [x_series.skew()],
    "Kurtosis (excess)":      [x_series.kurtosis()],
    "min(x)":                 [x_series.min()],
    "max(x)":                 [x_series.max()],
    "Range":                  [x_series.max()-x_series.min()],
}).T.rename(columns={0:"Value"}).round(6)

summary

7Correlation Analysis¶

In signal processing, correlation analysis identifies similarities between sensor measurements, detects periodicities, and reduces noise effects.

7.1Similarity of Signals¶

Similarity between two zero-mean time-series can be formalized by the inner product (or covariance):

\langle \mathbf{x}, \mathbf{y} \rangle \;=\; \text{Cov}(\mathbf{x}, \mathbf{y}) = \frac{\sum_{n=0}^{N-1} x[n]\,y[n]}{N},

(17)

The induced norm $\|\mathbf{x}\|=\sqrt{\langle \mathbf{x},\mathbf{x}\rangle} = \sigma_x$ yields a scale-invariant similarity (cosine similarity or Pearson’s correlation coefficient)

\rho_{xy} \;=\; \frac{\langle \mathbf{x}, \mathbf{y} \rangle}{\|\mathbf{x}\|\,\|\mathbf{y}\|} = \frac{\text{Cov}(\mathbf{x}, \mathbf{y})}{\sigma_x . \sigma_y}\in[-1,1],

(18)

which measures how much $x[n]$ aligns with $y[n]$ . Large $|\rho_{xy}|$ means the signals agree strongly in phase and shape at zero lag; $\rho_{xy}=0$ indicates orthogonality (no linear similarity).

<Figure size 1800x1050 with 4 Axes> — Figure 16:Illustrative scatter plots of correlation. (a) Two uncorrelated signals (ρ = 0). (b) Two perfectly correlated signals (ρ = 1).

Additional Learning: Correlation Matrix

For signals measured at $N_s$ different locations, $x_1[n], x_2[n], x_3[n], \ldots$ , we construct a correlation matrix $\mathbf{C}$ where each element $C_{ij} = \langle x_i[n], x_j[n] \rangle$ is the correlation between locations $i$ and $j$ .

In practice, the number of independent random variables driving these signals is often smaller than the number of sensors. Suppose we have a set of zero-mean independent stochastic signals $\mathbf{e}(t)$ with unit variance, so $\mathbb{E}[\mathbf{e}(t) \mathbf{e}^T(t)] = \mathbf{I}$ . Multiplying by a matrix $\mathbf{V}$ generates a linear combination of these independent sources:

\mathbf{y}(t) = \mathbf{V} \mathbf{e}(t)

(19)

This transformation shows how a small set of independent random variables can drive a larger set of signals, and it forms the basis of dimensionality reduction techniques such as Principal Component Analysis (PCA).

7.2Auto/Cross-correlation¶

Cross-correlation measures the similarity between two signals as a function of time lag ( $\tau=hT_s$ ):

r_{xy}[h]\;=\;\frac{\sum_{n=0}^{N-1} x[n]\,y[n+h]}{N},\qquad h\in\mathbb{Z}.

(20)

$r_{xy}[h]$ slides $y[n]$ by $h$ and computes the inner product at each shift; peaks indicate time alignments where the two signals share structure even when not synchronized at $h=0$ .

Figure 17:Cross-correlation operation’s intuition [Source].

A normalized measure that removes scale is:

\rho_{xy}[h]\;=\;\frac{r_{xy}[h]}{\sqrt{r_{xx}[0]\,r_{yy}[0]}}\in[-1,1],

(21)

1
2
3
4
5
6
7
8
9
10
import numpy as np
from scipy import signal

# Method 1: Using np.correlate (most common)
xcorr_np = np.correlate(x, y, mode='full')
xcorr_norm = xcorr_np / np.sqrt(np.max(np.correlate(x,x,mode='full')) * np.max(np.correlate(y,y,mode='full')))

# Method 2: Using scipy.signal.correlate
xcorr_scipy = signal.correlate(x, y, mode='full')
xcorr_scipy_norm = xcorr_scipy / xcorr_scipy.max()

The autocorrelation is the special case $y=x$ . The autocorrelation function is symmetric; its value at zero lag equals the variance ( $\sigma_x^2$ ), and it approaches zero for large lags.

r_{xx}[h]\;=\;\mathbb{E}\!\{X_t . X_{t+\tau}\},\quad r_{xx}[0]=\mathbb{E}\{X^{2}\}\;\;(=\sigma_x \text{ for zero-mean } X),

(22)

The following figure demonstrates autocorrelation by comparing a signal with its time-lagged copies.

<Figure size 2250x1500 with 4 Axes> — Figure 18:Demonstration of autocorrelation showing signal alignment at different lags and the complete autocorrelation function

7.3Period Detection Using Autocorrelation¶

Autocorrelation is useful for detecting the period of a signal when periodicity is not obvious from the time-domain waveform. When a signal contains periodic components, the autocorrelation function shows peaks at lags equal to the fundamental period and its multiples, because shifting a periodic signal by one complete period produces perfect alignment and maximum correlation. For a signal with period $T_0$ , peaks appear at lags $h = kT_0$ for $k = 0, \pm 1, \pm 2, \ldots$ ; the spacing between consecutive peaks directly reveals $T_0$ , even in the presence of noise.

import matplotlib.pyplot as plt

# Configure for clean plots
plt.rcParams.update({
    'font.size': 11,
    'axes.labelsize': 12,
    'axes.titlesize': 13,
    'xtick.labelsize': 10,
    'ytick.labelsize': 10,
    'legend.fontsize': 10,
    'figure.titlesize': 14,
    'lines.linewidth': 2.0,
    'axes.linewidth': 1.2,
    'grid.linewidth': 0.8,
    'font.family': 'sans-serif',
    'font.sans-serif': ['Arial', 'Liberation Sans', 'Nimbus Sans', 'Nimbus Sans L', 'Helvetica', 'DejaVu Sans'],
    'mathtext.fontset': 'dejavusans'
})

<Figure size 2100x1500 with 2 Axes> — Figure 19:Period detection using autocorrelation analysis of a composite sinusoidal signal

8Noise Modeling and Rejection Techniques¶

Noise is inherent in sensor measurements due to electrical and environmental effects, thermal fluctuations, quantization errors, and electromagnetic interference. What we measure in practice is the corrupted signal:

x[n] = s[n] + n_s[n]

(23)

where $s[n]$ is the desired noise-free signal and $n_s[n]$ is the additive noise component. Correlation analysis provides tools to reduce noise effects and recover information about the original signal structure.

8.1Signal-to-Noise Ratio (SNR)¶

Measurement quality is commonly quantified using the Signal-to-Noise Ratio (SNR), defined as the ratio of signal power to noise power:

\text{SNR} = \frac{\sigma_s^2}{\sigma_{n_s}^2}

(24)

The SNR is typically expressed in decibels (dB):

\text{SNR}_{\text{dB}} = 10 \log_{10}\left(\frac{\sigma_s^2}{\sigma_{n_s}^2}\right) \text{ dB}

(25)

Higher SNR values indicate cleaner measurements; lower values indicate significant noise contamination.

8.2White noise¶

A zero-mean Gaussian white-noise process $n_s[n]$ with variance $\sigma_{n_s}^2$ is defined by a delta autocorrelation:

r_{n_sn_s}[h]\;=\;\sigma_{n_s}^{2}\,\delta[h],\qquad \delta[0]=1,\;\delta[h\neq 0]=0,

(26)

If $n_{s1}[n]$ and $n_{s2}[n]$ are independent white noises, then their cross-correlation is identically zero:

r_{n_{s1} n_{s2}}[h]\;=\;0\quad\forall h.

(27)

Likewise, if $x[n]$ is any WSS process independent of $n_s[n]$ , then

r_{xn_s}[h]\;=\;\mathbb{E}\{x[n]\,n_s[n+h]\} \;=\;0\quad\forall h,

(28)

import matplotlib.pyplot as plt
import seaborn as sns

# Set publication-quality style
plt.style.use('seaborn-v0_8-whitegrid')
sns.set_palette("Set2")

# Configure matplotlib for publication quality
plt.rcParams.update({
    'font.size': 11,
    'axes.labelsize': 12,
    'axes.titlesize': 13,
    'xtick.labelsize': 10,
    'ytick.labelsize': 10,
    'legend.fontsize': 10,
    'figure.titlesize': 14,
    'lines.linewidth': 2.0,
    'axes.linewidth': 1.2,
    'grid.linewidth': 0.8,
    'font.family': 'sans-serif',
    'font.sans-serif': ['Arial', 'Liberation Sans', 'Nimbus Sans', 'Nimbus Sans L', 'Helvetica', 'DejaVu Sans'],
    'mathtext.fontset': 'dejavusans'
})

8.3Other Common Noise Models¶

8.3.1Pink Noise (1/f Noise)¶

Unlike white noise, in which each sample is independent, pink noise exhibits temporal correlations that decay gradually with increasing lag following a power law. Distant samples retain a weak but measurable correlation. This slow decay means that power integrated over logarithmic frequency bands remains constant, hence the alternative name “equal-energy-per-octave” noise. In SHM, pink noise arises from environmental loading, temperature fluctuations, and electronic drift when multiple physical processes with different time constants contribute to the signal.

r_{n_p n_p}[h] = \sigma_{n_p}^2 \cdot \frac{1}{|h|^{\alpha} + 1}, \quad 0 < \alpha < 1

(29)

8.3.2Random Walk¶

Red noise (random walk) is the cumulative sum of white noise increments. Each sample equals the accumulated sum of all previous perturbations, creating a non-stationary process with unbounded variance. The generation mechanism is a simple recursion: the current value equals the previous value plus a white noise innovation. In structural monitoring, red noise models foundation settlement, long-term material creep, sensor baseline drift, and accumulated measurement errors. Nearby samples are highly correlated because they share most of their accumulated history; correlation decreases as lag grows and fewer common increments contribute to both samples.

n_r[n] = n_r[n-1] + w[n], \quad w[n] \sim \mathcal{N}(0, \sigma_w^2)

(30)

r_{n_r n_r}[h] = \sigma_{n_r}^2 \left(1 - \frac{|h|}{N}\right), \quad |h| < N

(31)

from matplotlib import rcParams

# Configure publication-quality settings
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Arial', 'Liberation Sans', 'Nimbus Sans', 'Nimbus Sans L', 'Helvetica', 'DejaVu Sans']
rcParams['font.size'] = 11
rcParams['axes.labelsize'] = 12
rcParams['axes.titlesize'] = 13
rcParams['xtick.labelsize'] = 10
rcParams['ytick.labelsize'] = 10
rcParams['legend.fontsize'] = 10
rcParams['figure.dpi'] = 150
rcParams['axes.linewidth'] = 0.8
rcParams['grid.linewidth'] = 0.4
rcParams['lines.linewidth'] = 1.5

<Figure size 2100x1200 with 4 Axes> — Figure 21:Time-domain signals and sample autocorrelation for pink noise and random walk noise.

8.4Using Autocorrelation for Noise Rejection¶

For autocorrelation analysis in the presence of white noise, assume that $s[n]$ is statistically independent of $n[n]$ (i.e., $r_{sn}[h] = 0$ for all lags $h$ ). By linearity of correlation and the white noise property from (26):

r_{xx}[h] = r_{ss}[h] + r_{nn}[h] = r_{ss}[h] + \sigma_n^{2}\,\delta[h]

(32)

This result has three consequences:

Zero-lag contribution: White noise adds only a delta spike at zero lag ( $h = 0$ ), contributing $\sigma_n^2$ to the autocorrelation.
Structure preservation: All non-zero-lag autocorrelation ( $h \neq 0$ ) comes exclusively from the underlying signal $s[n]$ .
Noise immunity: The signal’s correlation structure remains intact despite time-domain corruption.

Even when $x[n]$ is heavily corrupted by noise, analyzing $r_{xx}[h]$ reveals the underlying signal structure, because the noise contribution vanishes at non-zero lags.

import matplotlib.pyplot as plt
import seaborn as sns

# Set publication-quality style
plt.style.use('seaborn-v0_8-whitegrid')
sns.set_palette("husl")

# Configure matplotlib for publication quality
plt.rcParams.update({
    'font.size': 11,
    'axes.labelsize': 12,
    'axes.titlesize': 13,
    'xtick.labelsize': 10,
    'ytick.labelsize': 10,
    'legend.fontsize': 10,
    'figure.titlesize': 14,
    'lines.linewidth': 2.0,
    'axes.linewidth': 1.2,
    'grid.linewidth': 0.8,
    'font.family': 'sans-serif',
    'font.sans-serif': ['Arial', 'Liberation Sans', 'Nimbus Sans', 'Nimbus Sans L', 'Helvetica', 'DejaVu Sans'],
    'mathtext.fontset': 'dejavusans'
})

<Figure size 2250x900 with 2 Axes> — Figure 22:Demonstration of noise rejection using autocorrelation analysis

9Impulse Response and System Characterization¶

9.1Simulation of Dynamic Systems using Impulse Response¶

Consider a single-degree-of-freedom (SDOF) dynamic system subjected to an external force $u(t)$ , with governing equation:

m\ddot{x}(t) + c\dot{x}(t) + kx(t) = u(t),

(33)

where $m$ is the mass, $c$ is the damping coefficient, and $k$ is the stiffness. The solution to this linear time-invariant system is expressed in terms of the impulse response function (IRF) $g(t)$ : the output $x(t)$ is the convolution of $u(t)$ with $g(t)$ :

x(t) \;=\; \int_{-\infty}^{\infty} u(\tau)\,g(t-\tau)\,d\tau \;=\; \int_{0}^{t} u(\tau)\,g(t-\tau)\,d\tau \;=\; (u * g) (t).

(34)

$g(t)$ is the displacement response to a unit impulse $u(t) = \delta(t)$ and fully characterizes the system by encoding its mass, damping, and stiffness properties. In discrete form:

x[n] = \sum_{h=0}^{N-1} u[h]g[n-h], \quad n = 0,1,\ldots, N-1

(35)

For an underdamped SDOF system, $g(t)$ can be obtained analytically:

g(t) \;=\; \frac{1}{m\omega_d}\,e^{-\zeta \omega_n t}\,\sin(\omega_d t), \qquad t \ge 0,

(36)

with parameters:

$\omega_n = \sqrt{k/m}$ , natural frequency,
$\zeta = c/(2\sqrt{km})$ , damping ratio,
$\omega_d = \omega_n \sqrt{1-\zeta^2}$ , damped natural frequency.

$g(t)$ is a decaying sinusoid: the oscillation period is set by $\omega_d$ and the envelope decay rate by $\zeta \omega_n$ .

Once the impulse response is known, the system response to any arbitrary input (wind, traffic, seismic excitation) is obtained by evaluating the convolution integral:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
import numpy as np

# System parameters
m, k, zeta = 1.0, (2*np.pi*2.0)**2, 0.02     # mass, stiffness, damping ratio
wn = np.sqrt(k/m)
wd = wn*np.sqrt(1-zeta**2)

# Time vector
dt, T = 0.01, 20
t = np.arange(0, T, dt)

# Impulse response g(t)
g = (1/(m*wd)) * np.exp(-zeta*wn*t) * np.sin(wd*t)

# Arbitrary input
u = 0.3*np.random.randn(len(t)) + np.sin(2*np.pi*1.5*t)*(t<5)

# Output response via convolution
x = np.convolve(u, g)[:N] * dt

from matplotlib import rcParams

# Set professional plotting parameters
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Arial', 'Liberation Sans', 'Nimbus Sans', 'Nimbus Sans L', 'Helvetica', 'DejaVu Sans']
rcParams['font.size'] = 11
rcParams['axes.linewidth'] = 1.2
rcParams['lines.linewidth'] = 1.8
rcParams['xtick.major.width'] = 1.2
rcParams['ytick.major.width'] = 1.2
rcParams['xtick.minor.width'] = 0.8
rcParams['ytick.minor.width'] = 0.8

<Figure size 1500x900 with 1 Axes> — Figure 23:A mpulse Response for a typical underdamped SDOF system.

9.2Impulse Response Identification Techniques¶

9.2.1Impact Test¶

In practice, $g(t)$ is obtained experimentally by applying an impulse (hammer impact), i.e., $u(t) = \delta(t)$ . Substituting into (34) and using the sifting property of the Dirac delta gives:

x(t) \;=\; \int_{0}^{t} \delta(\tau)\,g(t-\tau)\,d\tau \;=\; g(t).

(37)

9.2.2Ambient Vibration Test¶

Suppose we measure only the output $x(t)$ under ambient excitations. If the ambient input is modeled as zero-mean white Gaussian noise, the output autocorrelation $r_{xx}(\tau)$ contains the same information as the impulse response $g(t)$ .

Starting from the autocorrelation of the output:

r_{xx}(\tau) \;=\; \mathbb{E}[\,x(t)\,x(t+\tau)\,].

(38)

Substituting (34) gives a double convolution:

r_{xx}(\tau) \;=\; \mathbb{E}\!\left[ \int g(\alpha)\,u(t-\alpha)\,d\alpha \; \int g(\beta)\,u(t+\tau-\beta)\,d\beta \right].

(39)

By linearity, $g(\alpha)$ and $g(\beta)$ move outside the expectation:

r_{xx}(\tau) \;=\; \int\!\!\int g(\alpha)\,g(\beta)\,\mathbb{E}[\,u(t-\alpha)\,u(t+\tau-\beta)\,]\,d\alpha\,d\beta.

(40)

The inner expectation is the autocorrelation of the input $u(t)$ :

\mathbb{E}[\,u(t-\alpha)\,u(t+\tau-\beta)\,] \;=\; \mathbb{E}[\,u(t)\,u(t+\tau + \alpha -\beta)\,] \;=\; r_{uu}(\tau + \alpha - \beta).

(41)

For zero-mean white Gaussian noise:

r_{uu}(\tau + \alpha - \beta) \;=\; \sigma_{uu}^2 \delta(\tau + \alpha - \beta).

(42)

Substituting into (40) and applying the sifting property of $\delta$ collapses one integral:

r_{xx}(\tau) \;=\; \sigma_{uu}^2 \int\!\!\int g(\alpha)\,g(\beta)\, \delta(\tau + \alpha - \beta)\,d\alpha\,d\beta \;=\; \sigma_{uu}^2 \int g(\alpha)\,g(\alpha+\tau)\,d\alpha.

(43)

This is the convolution of $g(t)$ with its time-reversal $g(-t)$ :

r_{xx}(\tau) \;=\; \sigma_{uu}^2 \,[g(t) * g(-t)].

(44)

Evaluating the convolution shows that $g_x(t) = g(t) * g(-t)$ retains the form of an exponentially decaying sinusoid: a constant scale factor times $e^{-\zeta \omega_n \tau}$ multiplied by a sinusoid at $\omega_d$ , possibly with a phase shift. Although the amplitude and shape differ from $g(t)$ , the decay rate still reflects the damping ratio $\zeta$ and the oscillation period is still governed by $\omega_d$ .

This result is the foundation of ambient vibration testing: exciting a structure with ambient white-like noise (traffic, wind) and measuring the output autocorrelation indirectly identifies the impulse response and hence the modal properties.

9.3Extracting Structure’s Parameters¶

Once $g(t)$ (or its estimate) is obtained, key system parameters follow:

Natural period: The oscillation period $T_n \approx T_d = 2\pi/\omega_d$ is read from zero-crossings or the dominant frequency of $g(t)$ .
Damping ratio: The exponential envelope decays as $e^{-\zeta \omega_n t}$ ; fitting an exponential to the peak amplitudes yields $\zeta$ .
System characterization: With $\omega_n$ and $\zeta$ known, the system’s response to arbitrary inputs follows by convolution.

The procedure is summarized below.

9.4Python Implementation¶

from matplotlib import rcParams

# Configure matplotlib for journal-quality figures
rcParams.update({
    'font.family': 'sans-serif',
    'font.sans-serif': ['Arial', 'Liberation Sans', 'Nimbus Sans', 'Nimbus Sans L', 'Helvetica', 'DejaVu Sans'],
    'font.size': 11,
    'axes.labelsize': 12,
    'axes.titlesize': 13,
    'legend.fontsize': 10,
    'xtick.labelsize': 10,
    'ytick.labelsize': 10,
    'lines.linewidth': 1.8,
    'axes.linewidth': 1.2,
    'xtick.major.width': 1.2,
    'ytick.major.width': 1.2,
    'xtick.minor.width': 0.8,
    'ytick.minor.width': 0.8,
    'figure.figsize': (12, 8),
    'axes.spines.top': False,
    'axes.spines.right': False,
    'axes.grid': True,
    'grid.alpha': 0.3,
    'grid.linewidth': 0.5
})

<Figure size 1800x1200 with 2 Axes> — Figure 24:Impulse response identification of an SDOF system under ambient load condition.

This implementation identifies the damped period (hence natural frequency) from autocorrelation zero-crossings and estimates damping from the exponential decay of the envelope. The resulting impulse response then predicts future response under arbitrary inputs via convolution.

9.5Impulse Response Functions for MDOF Systems¶

Consider a bridge structure discretized into $N_d$ nodes along its span, each node representing a degree of freedom. The equation of motion has the same form as (33) but in matrix form. An impulse at node $j$ produces a displacement response at node $i$ characterized by $g_{ij}(t)$ , yielding an $N_d \times N_d$ impulse response matrix that completely describes the structural dynamics.

Figure 25:Multi-degree-of-freedom bridge structure showing impulse response function setup.

The complete impulse response matrix $\mathbf{g}(t)$ contains all possible input-output relationships:

\mathbf{g}(t) = \begin{bmatrix} g_{11}(t) & g_{12}(t) & \cdots & g_{1N_d}(t) \\ g_{21}(t) & g_{22}(t) & \cdots & g_{2N_d}(t) \\ \vdots & \vdots & \ddots & \vdots \\ g_{N_d1}(t) & g_{Nd2}(t) & \cdots & g_{N_d N_d}(t) \end{bmatrix}

(45)

Maxwell’s reciprocity theorem, arising from the symmetry of the mass and stiffness matrices, gives:

g_{ij}(t) = g_{ji}(t)

(46)

This reduces the number of independent impulse response functions from $N_d^2$ to $N_d(N_d+1)/2$ , simplifying experimental characterization.

9.5.1Response to Arbitrary Loading¶

When the structure is subjected to arbitrary forces $u_1(t), u_2(t), \ldots, u_{N_d}(t)$ at nodes 1 through $N_d$ , the response at node $i$ is:

x_i(t) = \sum_{j=1}^{N_d} \int_{0}^{\infty} u_j(\tau) g_{ij}(t-\tau) d\tau = \sum_{j=1}^{N_d} (u_j * g_{ij})(t)

(47)

This superposition enables response prediction for complex loading scenarios from the impulse response matrix alone. The complete IRF for an MDOF system is:

g_{ij}(t) = \sum_{r=1}^{N_d} \frac{\phi_{ir}\phi_{jr}}{M_r\omega_{dr}} e^{-\zeta_r\omega_{nr}t} \sin(\omega_{dr}t), \quad t \geq 0

(48)

where $\phi_{ir}$ and $\phi_{jr}$ are the $r$ -th mode shape ordinates at nodes $i$ and $j$ , $M_r$ is the modal mass, $\omega_{nr}$ is the natural frequency, $\zeta_r$ is the damping ratio, and $\omega_{dr} = \omega_{nr}\sqrt{1-\zeta_r^2}$ is the damped natural frequency for mode $r$ . Each IRF is a superposition of $N_d$ modal contributions; each modal term contains:

Modal scaling factor: $\frac{\phi_{ir}\phi_{jr}}{M_r\omega_{dr}}$ determines the relative contribution magnitude
Exponential decay: $e^{-\zeta_r\omega_{nr}t}$ governs the rate of amplitude reduction
Oscillatory component: $\sin(\omega_{dr}t)$ provides the frequency content

9.5.2Practical Implications¶

Experimental modal analysis: The expression provides the theoretical basis for extracting modal parameters ( $\omega_{nr}$ , $\zeta_r$ , $\boldsymbol{\phi}_{r}$ ) from measured impulse response functions by curve fitting.

Damage detection: Structural damage alters modal parameters and therefore the IRFs; comparing baseline and current IRFs reveals these changes.

Model reduction: For lightly damped structures, only the first few modes contribute substantially to the impulse response, permitting reduced-order modeling.

The following figure shows how multiple modes combine to form the complete impulse response function. Multi-modal superposition complicates parameter identification because the measured IRF at each node contains overlapping frequencies with varying amplitudes and decay rates. Frequency-domain analysis is generally required to extract individual modal parameters from the composite response.

from matplotlib import rcParams

# Set professional plotting parameters
rcParams['font.family'] = 'sans-serif'
rcParams['font.sans-serif'] = ['Arial', 'Liberation Sans', 'Nimbus Sans', 'Nimbus Sans L', 'Helvetica', 'DejaVu Sans']
rcParams['font.size'] = 11
rcParams['axes.linewidth'] = 1.2
rcParams['lines.linewidth'] = 1.8

10Exercises¶

The following exercises reinforce the concepts covered in this chapter. They include both conceptual questions and coding tasks. (NumPy/Pandas are assumed available for coding exercises.)

Solution to Exercise 4 #

import numpy as np
import matplotlib.pyplot as plt

# time vector
t = np.linspace(0, 100, 1001)
# components
trend = 0.1 * t  # linear trend
season = np.sin(2*np.pi*0.1 * t)  # sinusoid with period 10 time units (freq 0.1 Hz)
noise = 0.5 * np.random.randn(len(t))
x = trend + season + noise

# Detrend by removing linear fit:
coefs = np.polyfit(t, x, 1)        # fit a line (degree 1 polynomial)
trend_fit = np.polyval(coefs, t)   # evaluate the line
x_detrended = x - trend_fit

# Plot original and detrended
plt.figure()
plt.subplot(2,1,1)
plt.plot(t, x, label='Original')
plt.plot(t, trend_fit, label='Fitted Trend', linestyle='--')
plt.legend()
plt.subplot(2,1,2)
plt.plot(t, x_detrended, label='Detrended')
plt.legend()
plt.show()

The original plot shows an upward trend with oscillation and noise; the detrended plot oscillates around zero with the linear growth removed. The sinusoidal component and noise remain.

Solution to Exercise 5 #

import numpy as np
import matplotlib.pyplot as plt

# generate noisy signal
t = np.linspace(0, 20, 201)
x = np.sin(0.5 * t) + 0.5 * np.random.randn(len(t))

# 5-point moving average filter
N = 5
kernel = np.ones(N) / N
x_smooth = np.convolve(x, kernel, mode='same')

plt.plot(t, x, label='Original')
plt.plot(t, x_smooth, label='5-pt Moving Avg', linewidth=2)
plt.legend(); plt.show()

The output shows a jagged original curve and a smoother version. The moving average smooths rapid changes over approximately 5 data points. The mode='same' argument preserves the output length (with some boundary effects).