In an MDOF mechanical structure, the (linear) vibration dynamics around an equilibrium are modeled by \( \mathbf{q}(t)\in\mathbb{R}^n \) and \( \mathbf{f}(t)\in\mathbb{R}^m \):

\[ \mathbf{M}\ddot{\mathbf{q}}(t)+\mathbf{C}\dot{\mathbf{q}}(t)+\mathbf{K}\mathbf{q}(t)=\mathbf{B}_f\,\mathbf{f}(t), \quad \mathbf{y}(t)=\mathbf{C}_y\mathbf{q}(t)+\mathbf{D}_y\dot{\mathbf{q}}(t)+\mathbf{v}(t) \]

Here \( \mathbf{y}(t) \) is a measured response (displacement/velocity/acceleration), and \( \mathbf{v}(t) \) represents sensor noise. EMA aims to estimate:

• Natural frequencies \( \omega_r \) and damping ratios \( \zeta_r \),
• Mode shapes \( \boldsymbol{\phi}_r \in \mathbb{R}^n \) (up to scaling),
• Optionally, a reduced-order state-space or transfer-function model consistent with measured data.

In the frequency domain, define the transfer matrix (FRF matrix) \( \mathbf{H}(j\omega) \) from forces to measured outputs:

\[ \mathbf{Y}(j\omega)=\mathbf{H}(j\omega)\mathbf{F}(j\omega),\quad \mathbf{H}(j\omega)=\mathbf{G}_{y}(j\omega)\mathbf{G}_{f}(j\omega)^{-1} \]

For the standard displacement FRF (outputs \( \mathbf{q} \), inputs \( \mathbf{f} \)), substituting the harmonic ansatz gives:

\[ \mathbf{Q}(j\omega)=\Big(\mathbf{K}- \omega^2\mathbf{M} + j\omega \mathbf{C}\Big)^{-1}\mathbf{B}_f\,\mathbf{F}(j\omega). \]

flowchart TD
  A["Define objective: modes or reduced model"] --> B["Select excitation: impact / shaker / broadband"]
  B --> C["Measure input f(t) and outputs y(t) (multi-sensor)"]
  C --> D["Preprocess: detrend, window, anti-alias, segment"]
  D --> E["Estimate spectra: G_uu, G_yy, G_yu"]
  E --> F["Compute FRFs and coherence"]
  F --> G["Extract modal parameters: poles, zeta, mode shapes"]
  G --> H["Validate: residuals, coherence, MAC, stability"]
        

Measurements are sampled at period \( T_s \) (sampling frequency \( f_s=1/T_s \)): \( y[k]=y(kT_s) \), \( u[k]=u(kT_s) \). A standard stochastic measurement model is:

\[ y[k] = y_0[k] + v[k], \quad u_m[k] = u[k] + w[k], \quad E[v[k]] = 0,\; E[w[k]] = 0 \]

Finite records introduce spectral leakage. With a window \( w[k] \), the discrete Fourier transform (DFT) of a length-\(N\) record is:

\[ X[\ell]=\sum_{k=0}^{N-1} x[k]\,w[k]\;e^{-j2\pi \ell k/N},\quad \omega_\ell = 2\pi \ell f_s/N. \]

Windowing reduces leakage but increases estimator variance. Hence EMA uses segmenting + averaging (Welch-type procedures) to trade bias/variance systematically.

For input \( u \) and output \( y \) (SISO notation for clarity), define autospectral and cross-spectral densities: \( G_{uu}(\omega) \), \( G_{yy}(\omega) \), \( G_{yu}(\omega) \). The most common estimators are:

\[ \hat{H}_1(\omega)=\frac{\hat{G}_{yu}(\omega)}{\hat{G}_{uu}(\omega)},\qquad \hat{H}_2(\omega)=\frac{\hat{G}_{yy}(\omega)}{\hat{G}_{uy}(\omega)},\qquad \hat{H}_v(\omega)=\sqrt{\hat{H}_1(\omega)\hat{H}_2(\omega)} \]

Interpretation: \( H_1 \) is robust when output noise dominates, while \( H_2 \) is robust when input noise dominates. A fundamental (and exam-relevant) unbiasedness statement is:

Proposition (unbiasedness of \(H_1\) under output noise only). Assume the true relation \( y_0 = H u \) and measured output \( y = y_0 + v \), where \( v \) is zero-mean, wide-sense stationary, and uncorrelated with \( u \). Then (in expectation) \( E[\hat{H}_1(\omega)] = H(\omega) \).

Proof sketch (frequency-domain moment factorization):

\[ G_{yu}(\omega)=G_{(y_0+v)u}(\omega)=G_{y_0u}(\omega)+G_{vu}(\omega) \]

Since \( y_0 = H u \), we have \( G_{y_0u}(\omega)=H(\omega)G_{uu}(\omega) \). Uncorrelatedness implies \( G_{vu}(\omega)=0 \). Therefore:

\[ \hat{H}_1(\omega)=\frac{\hat{G}_{yu}(\omega)}{\hat{G}_{uu}(\omega)} \approx \frac{H(\omega)G_{uu}(\omega)}{G_{uu}(\omega)} = H(\omega). \]

In contrast, if input is noisy \( u_m=u+w \), then \( \hat{G}_{u_mu_m} \) includes the noise spectrum, biasing \(H_1\). This is the practical reason to choose \(H_1\) vs. \(H_2\).

The magnitude-squared coherence is defined by:

\[ \gamma^2(\omega)=\frac{|G_{yu}(\omega)|^2}{G_{uu}(\omega)\,G_{yy}(\omega)}. \]

A key inequality (Cauchy–Schwarz) shows \( 0 \le \gamma^2(\omega) \le 1 \). When \( \gamma^2(\omega)\approx 1 \), output is (nearly) linearly explained by the input at that frequency; low coherence indicates nonlinearities, leakage, unmeasured disturbances, or poor excitation energy.

Proof idea: treat \(G_{yu}(\omega)\) as an inner product of Fourier coefficients and apply Cauchy–Schwarz, yielding \( |G_{yu}|^2 \le G_{uu}G_{yy} \).

Suppose the system is lightly damped and well described by a modal expansion. For a collocated input/output (force-to-displacement), a commonly used representation is:

\[ \mathbf{H}(\omega) \approx \sum_{r=1}^{n} \frac{\boldsymbol{\phi}_r\boldsymbol{\phi}_r^{\mathsf{T}}}{\omega_r^2-\omega^2 + 2j\zeta_r\omega_r\omega} \]

Each mode contributes a rational “peak” centered near \( \omega_r \) whose width scales with \( \zeta_r \). In practice, EMA extracts \((\omega_r,\zeta_r)\) and estimates \( \boldsymbol{\phi}_r \) by fitting the complex FRF data across frequency.

A widely taught introductory method is peak-picking + half-power bandwidth, justified by the SDOF FRF:

\[ H(\omega)=\frac{1}{k}\cdot\frac{1}{1-(\omega/\omega_n)^2 + j\,2\zeta(\omega/\omega_n)}. \]

Near resonance (\( \omega\approx \omega_n \)) and for light damping, the half-power frequencies \(\omega_1,\omega_2\) where \( |H(\omega_{1,2})|=|H(\omega_n)|/\sqrt{2} \) satisfy:

\[ \zeta \approx \frac{\omega_2-\omega_1}{2\omega_n} \quad \text{(valid when } 0 < \zeta < 0.1 \text{ and modes are well separated).} \]

A control-oriented model is often written in state-space form (continuous time):

\[ \dot{\mathbf{x}}(t)=\mathbf{A}\mathbf{x}(t)+\mathbf{B}\mathbf{u}(t)+\mathbf{w}(t),\qquad \mathbf{y}(t)=\mathbf{C}\mathbf{x}(t)+\mathbf{D}\mathbf{u}(t)+\mathbf{v}(t). \]

The transfer matrix is \( \mathbf{G}(s)=\mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B}+\mathbf{D} \), and the FRF is simply \( \mathbf{G}(j\omega) \). Thus, “EMA in the frequency domain” and “system identification in the time domain” are two views of the same LTI object.

System identification selects a parametric family \( \mathbf{G}(s;\boldsymbol{\theta}) \) (or \((\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D})\)) and estimates parameters by minimizing a prediction error:

\[ \hat{\boldsymbol{\theta}}=\arg\min_{\boldsymbol{\theta}} \sum_{k=0}^{N-1}\big\| \varepsilon[k;\boldsymbol{\theta}]\big\|^2, \quad \varepsilon[k;\boldsymbol{\theta}]=y[k]-\hat{y}[k\mid \boldsymbol{\theta}] \]

Under Gaussian noise assumptions, minimizing squared prediction errors is equivalent to maximum likelihood. A central asymptotic result (large \(N\)) is that estimation error covariance is governed by a Fisher-information-like matrix, giving a principled uncertainty quantification for identified modal parameters.

flowchart TD
  A["Collect u[k], y[k] (sampled)"] --> B["Choose model class: state-space / transfer / reduced modal"]
  B --> C["Define criterion: prediction error (least-squares) or likelihood"]
  C --> D["Estimate parameters theta (optimization / subspace / ERA)"]
  D --> E["Validate: residual whiteness, fit on new data, coherence"]
  E --> F["Use model: modal parameters, controller design, simulation"]
        

The following script simulates a 3-DOF chain, estimates the FRF using the H1 estimator with Welch averaging, computes coherence, extracts rough modal estimates by half-power bandwidth, and demonstrates a compact ERA from an impulse response.

Chapter13_Lesson5.py

"""
Chapter 13 - Vibrations and Multi-Degree-of-Freedom (MDOF) Systems
Lesson 5 - Intro to Experimental Modal Analysis and System Identification Concepts

File: Chapter13_Lesson5.py

This script:
1) Simulates a 3-DOF mass-spring-damper system under broadband force excitation.
2) Estimates frequency response functions (FRFs) using the H1 estimator via Welch spectral averages.
3) Computes magnitude-squared coherence.
4) Extracts rough modal estimates (natural frequencies + damping) via peak-picking + half-power bandwidth.
5) Demonstrates a compact Eigensystem Realization Algorithm (ERA) identification from an impulse response.

Dependencies: numpy, scipy, matplotlib
"""

from __future__ import annotations
import numpy as np
from numpy.typing import NDArray
from dataclasses import dataclass
import scipy.signal as sig
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt


# ----------------------------
# 1) Mechanical model (3-DOF)
# ----------------------------
@dataclass
class MDOF:
    M: NDArray[np.float64]
    C: NDArray[np.float64]
    K: NDArray[np.float64]

    def to_state_space(self):
        n = self.M.shape[0]
        Minv = np.linalg.inv(self.M)
        Z = np.zeros((n, n))
        I = np.eye(n)

        A = np.block([
            [Z, I],
            [-Minv @ self.K, -Minv @ self.C],
        ])

        # input: force vector f(t) in physical coordinates
        B = np.vstack([np.zeros((n, n)), Minv])

        # outputs: displacement and acceleration (choose later)
        Cx = np.block([I, Z])                           # displacement
        Cv = np.block([Z, I])                           # velocity
        Ca = np.block([-Minv @ self.K, -Minv @ self.C])  # acceleration (assuming output = qdd)

        D = np.zeros((n, n))
        return A, B, Cx, Cv, Ca, D


def make_chain_3dof() -> MDOF:
    # masses
    m1, m2, m3 = 1.0, 0.9, 0.8
    M = np.diag([m1, m2, m3])

    # springs: grounded at DOF1 and DOF3, coupling springs between masses
    k1, k2, k3, k4 = 800.0, 600.0, 500.0, 700.0
    # K for chain: (ground)-k1-m1-k2-m2-k3-m3-k4-(ground)
    K = np.array([
        [k1 + k2, -k2,      0.0],
        [-k2,     k2 + k3, -k3 ],
        [0.0,     -k3,     k3 + k4],
    ])

    # viscous damping (Rayleigh-like, but specified directly for simplicity)
    c1, c2, c3, c4 = 2.0, 1.8, 1.5, 2.2
    C = np.array([
        [c1 + c2, -c2,      0.0],
        [-c2,     c2 + c3, -c3 ],
        [0.0,     -c3,     c3 + c4],
    ])
    return MDOF(M=M, C=C, K=K)


# ---------------------------------------------
# 2) Time simulation under broadband excitation
# ---------------------------------------------
def simulate(mdof: MDOF, fs: float, T: float, force_dof: int = 0, seed: int = 7):
    """
    Simulate with band-limited white noise force on one DOF.

    Returns:
        t (N,), u (N,), y (N,) where y is acceleration at a measured DOF (same as force_dof by default).
    """
    rng = np.random.default_rng(seed)
    dt = 1.0 / fs
    t = np.arange(0.0, T, dt)
    N = len(t)

    # force: band-limited white noise (shape through lowpass filter)
    u = rng.normal(0.0, 1.0, size=N)
    b, a = sig.butter(4, 0.25)  # cutoff at 0.25*(Nyquist) => 0.125*fs
    u = sig.filtfilt(b, a, u)
    u = u / np.std(u)

    A, B, Cx, Cv, Ca, D = mdof.to_state_space()
    n = mdof.M.shape[0]

    # Input vector f(t) has dimension n; excite one coordinate
    def rhs(ti, x):
        # interpolate input
        k = int(np.clip(np.floor(ti / dt), 0, N - 1))
        f = np.zeros(n)
        f[force_dof] = u[k]
        return (A @ x + B @ f)

    x0 = np.zeros(2 * n)
    sol = solve_ivp(rhs, (t[0], t[-1]), x0, t_eval=t, method="RK45", rtol=1e-7, atol=1e-9)
    X = sol.y.T  # (N, 2n)

    # acceleration output at measured DOF (here: same DOF as force)
    y = (Ca @ X.T).T[:, force_dof]

    # add mild measurement noise
    y = y + 0.02 * rng.normal(size=N)

    return t, u, y


# ----------------------------------------------------
# 3) FRF estimation: Welch averages + H1 + coherence
# ----------------------------------------------------
def estimate_frf_h1(u: NDArray[np.float64], y: NDArray[np.float64], fs: float,
                    nperseg: int = 4096, noverlap: int = 2048):
    """
    H1 estimator:
        H1(w) = G_yu(w) / G_uu(w)

    Coherence:
        gamma^2(w) = |G_yu(w)|^2 / (G_uu(w) G_yy(w))
    """
    f, G_uu = sig.welch(u, fs=fs, nperseg=nperseg, noverlap=noverlap, return_onesided=True)
    f, G_yy = sig.welch(y, fs=fs, nperseg=nperseg, noverlap=noverlap, return_onesided=True)
    f, G_yu = sig.csd(y, u, fs=fs, nperseg=nperseg, noverlap=noverlap, return_onesided=True)

    H1 = G_yu / G_uu
    coh = np.abs(G_yu) ** 2 / (G_uu * G_yy + 1e-30)
    return f, H1, coh


# ----------------------------------------------------
# 4) Peak-picking + half-power damping (rough)
# ----------------------------------------------------
def half_power_damping(f: NDArray[np.float64], Hmag: NDArray[np.float64], peak_idx: int):
    """
    For an SDOF-like isolated peak:
        zeta ≈ (f2 - f1) / (2 fn)
    where f1 and f2 are the half-power frequencies (|H| drops to 1/sqrt(2) of peak).
    """
    fn = f[peak_idx]
    peak = Hmag[peak_idx]
    target = peak / np.sqrt(2.0)

    # search left
    i1 = peak_idx
    while i1 > 1 and Hmag[i1] > target:
        i1 -= 1
    # linear interpolate
    f1 = np.interp(target, [Hmag[i1], Hmag[i1 + 1]], [f[i1], f[i1 + 1]])

    # search right
    i2 = peak_idx
    while i2 < len(f) - 2 and Hmag[i2] > target:
        i2 += 1
    f2 = np.interp(target, [Hmag[i2 - 1], Hmag[i2]], [f[i2 - 1], f[i2]])

    zeta = (f2 - f1) / (2.0 * fn + 1e-30)
    return fn, zeta, f1, f2


def extract_modes(f: NDArray[np.float64], H1: NDArray[np.complex128], coh: NDArray[np.float64],
                  fmax: float = 80.0, min_coh: float = 0.6):
    """
    Simple peak-picking on |H1|, filtering by coherence.
    """
    mask = (f > 0.5) & (f < fmax) & (coh > min_coh)
    f2 = f[mask]
    Hm = np.abs(H1[mask])

    peaks, props = sig.find_peaks(Hm, prominence=np.percentile(Hm, 70), distance=10)
    est = []
    for p in peaks[:3]:  # expect ~3 modes
        fn, zeta, f1, f2hp = half_power_damping(f2, Hm, p)
        est.append((fn, zeta, f1, f2hp))
    return est


# ----------------------------------------------------
# 5) ERA identification from impulse response (SISO)
# ----------------------------------------------------
def simulate_impulse(mdof: MDOF, fs: float, T: float, in_dof: int = 0, out_dof: int = 0):
    """
    Approximate impulse response by applying a narrow pulse input in continuous time.
    Output is acceleration at out_dof.
    """
    dt = 1.0 / fs
    t = np.arange(0.0, T, dt)
    N = len(t)

    # pulse area ~1 to mimic impulse
    u = np.zeros(N)
    u[0] = 1.0 / dt

    A, B, Cx, Cv, Ca, D = mdof.to_state_space()
    n = mdof.M.shape[0]

    def rhs(ti, x):
        k = int(np.clip(np.floor(ti / dt), 0, N - 1))
        f = np.zeros(n)
        f[in_dof] = u[k]
        return (A @ x + B @ f)

    x0 = np.zeros(2 * n)
    sol = solve_ivp(rhs, (t[0], t[-1]), x0, t_eval=t, method="RK45", rtol=1e-8, atol=1e-10)
    X = sol.y.T
    y = (Ca @ X.T).T[:, out_dof]
    return t, y


def era_siso(y: NDArray[np.float64], s: int, r: int, order: int):
    """
    Eigensystem Realization Algorithm (ERA) for SISO impulse response.

    Inputs:
        y[k] = Markov parameters (impulse response samples), k=0..N-1
        s, r : Hankel block sizes
        order: desired reduced order (model order)

    Returns:
        A, B, C (discrete-time), and singular values.
    """
    # Build Hankel matrices H0 and H1
    # H0[i,j] = y[i+j], H1[i,j] = y[i+j+1]
    H0 = np.zeros((s, r))
    H1 = np.zeros((s, r))
    for i in range(s):
        for j in range(r):
            H0[i, j] = y[i + j]
            H1[i, j] = y[i + j + 1]

    U, S, Vt = np.linalg.svd(H0, full_matrices=False)
    U1 = U[:, :order]
    S1 = np.diag(S[:order])
    V1 = Vt[:order, :]

    S1_sqrt = np.sqrt(S1)
    S1_isqrt = np.linalg.inv(S1_sqrt)

    A = S1_isqrt @ (U1.T @ H1 @ V1.T) @ S1_isqrt
    B = S1_sqrt @ V1[:, [0]]    # first column corresponds to first Markov parameter
    C = U1[[0], :] @ S1_sqrt    # first row corresponds to first output
    return A, B, C, S


def main():
    fs = 500.0   # Hz
    T = 40.0     # seconds
    mdof = make_chain_3dof()

    t, u, y = simulate(mdof, fs=fs, T=T, force_dof=0, seed=7)
    f, H1, coh = estimate_frf_h1(u, y, fs=fs, nperseg=4096, noverlap=2048)

    # Extract rough modal estimates
    est = extract_modes(f, H1, coh, fmax=80.0, min_coh=0.6)
    print("Rough modal estimates (fn Hz, zeta, half-power f1,f2):")
    for (fn, zeta, f1, f2hp) in est:
        print(f"  fn={fn:7.3f} Hz, zeta={zeta:7.4f}, f1={f1:7.3f}, f2={f2hp:7.3f}")

    # Plots
    plt.figure()
    plt.semilogy(f, np.abs(H1))
    plt.xlim(0, 120)
    plt.xlabel("Frequency [Hz]")
    plt.ylabel("|H1(jw)| [acc/force]")
    plt.title("Estimated FRF (H1)")

    plt.figure()
    plt.plot(f, coh)
    plt.ylim(0, 1.05)
    plt.xlim(0, 120)
    plt.xlabel("Frequency [Hz]")
    plt.ylabel("Coherence")
    plt.title("Magnitude-squared coherence")

    # ERA from impulse response (SISO)
    t_imp, y_imp = simulate_impulse(mdof, fs=fs, T=10.0, in_dof=0, out_dof=0)
    s = 150
    r = 150
    order = 6  # 2n states for 3-DOF, but measured accel is SISO; try reduced order
    A, B, C, Svals = era_siso(y_imp, s=s, r=r, order=order)

    eigvals = np.linalg.eigvals(A)
    # Convert discrete poles to continuous (approx): lambda_c = ln(z)/dt
    dt = 1.0 / fs
    lambdas = np.log(eigvals) / dt
    wn = np.abs(lambdas) / (2.0 * np.pi)
    zeta = -np.real(lambdas) / (np.abs(lambdas) + 1e-30)

    print("\nERA (discrete-time) identified poles (continuous approx):")
    for i in range(len(lambdas)):
        print(f"  pole {i+1}: lambda={lambdas[i]: .4e}, fn~{wn[i]:.3f} Hz, zeta~{zeta[i]:.4f}")

    plt.figure()
    plt.semilogy(Svals, marker="o")
    plt.xlabel("Index")
    plt.ylabel("Singular value")
    plt.title("ERA Hankel singular values")

    plt.show()


if __name__ == "__main__":
    main()
      

This C++ code is intentionally compact and uses a direct DFT (O(N^2)) to keep the implementation self-contained. For serious work, replace the DFT with FFT (FFTW) and add segmentation/averaging and windowing.

Chapter13_Lesson5.cpp

/*
Chapter 13 - Vibrations and Multi-Degree-of-Freedom (MDOF) Systems
Lesson 5 - Intro to Experimental Modal Analysis and System Identification Concepts

File: Chapter13_Lesson5.cpp

This is a compact, didactic C++ example that:
1) Simulates a 3-DOF state-space model with RK4 integration.
2) Estimates an FRF using the H1 estimator via a simple DFT (O(N^2)).

Notes:
- For production, replace the DFT with FFT (e.g., FFTW) and use robust windowing/averaging.
- Requires Eigen (header-only): https://eigen.tuxfamily.org
Build (example):
  g++ -O2 -std=c++17 Chapter13_Lesson5.cpp -I path/to/eigen -o Chapter13_Lesson5
*/

#include <Eigen/Dense>
#include <complex>
#include <iostream>
#include <vector>
#include <cmath>
#include <random>

using Eigen::MatrixXd;
using Eigen::VectorXd;

static constexpr double PI = 3.14159265358979323846;

struct Spectra {
    std::vector<double> f;
    std::vector<std::complex<double>> Guu;
    std::vector<std::complex<double>> Gyy;
    std::vector<std::complex<double>> Gyu;
};

std::vector<std::complex<double>> dft(const std::vector<double>& x) {
    const int N = (int)x.size();
    std::vector<std::complex<double>> X(N);
    for (int k = 0; k < N; ++k) {
        std::complex<double> s(0.0, 0.0);
        for (int n = 0; n < N; ++n) {
            double ang = -2.0 * PI * k * n / (double)N;
            s += std::complex<double>(std::cos(ang), std::sin(ang)) * x[n];
        }
        X[k] = s;
    }
    return X;
}

Spectra spectra_single_segment(const std::vector<double>& u, const std::vector<double>& y, double fs) {
    // Single-segment (no averaging) autospectra and cross-spectrum using DFT:
    // Guu(k) = (1/N) U(k) conj(U(k)), etc.
    const int N = (int)u.size();
    auto U = dft(u);
    auto Y = dft(y);

    Spectra out;
    const int K = N/2 + 1;
    out.f.resize(K);
    out.Guu.resize(K);
    out.Gyy.resize(K);
    out.Gyu.resize(K);

    for (int k = 0; k < K; ++k) {
        out.f[k] = (fs * k) / (double)N;
        std::complex<double> Uc = std::conj(U[k]);
        std::complex<double> Yc = std::conj(Y[k]);
        out.Guu[k] = (U[k] * Uc) / (double)N;
        out.Gyy[k] = (Y[k] * Yc) / (double)N;
        out.Gyu[k] = (Y[k] * Uc) / (double)N;
    }
    return out;
}

int main() {
    // 3-DOF chain example
    MatrixXd M = MatrixXd::Zero(3,3);
    M(0,0) = 1.0; M(1,1) = 0.9; M(2,2) = 0.8;

    double k1=800, k2=600, k3=500, k4=700;
    MatrixXd K(3,3);
    K << k1+k2, -k2,    0,
         -k2,   k2+k3, -k3,
          0,    -k3,   k3+k4;

    double c1=2.0, c2=1.8, c3=1.5, c4=2.2;
    MatrixXd C(3,3);
    C << c1+c2, -c2,    0,
         -c2,   c2+c3, -c3,
          0,    -c3,   c3+c4;

    // State-space: x=[q; qd], xdot = A x + B f
    MatrixXd Minv = M.inverse();
    MatrixXd Z = MatrixXd::Zero(3,3);
    MatrixXd I = MatrixXd::Identity(3,3);

    MatrixXd A(6,6);
    A << Z, I,
         -Minv*K, -Minv*C;

    MatrixXd B = MatrixXd::Zero(6,3);
    B.block(3,0,3,3) = Minv;

    // Acceleration output: qdd = -Minv*K q - Minv*C qd + Minv f
    MatrixXd Ca(3,6);
    Ca << -Minv*K, -Minv*C;

    const double fs = 500.0;
    const double dt = 1.0 / fs;
    const double T = 10.0;
    const int N = (int)std::floor(T * fs);

    // Input: band-limited noise (very simple smoothing)
    std::mt19937 gen(7);
    std::normal_distribution<double> nd(0.0, 1.0);

    std::vector<double> u(N, 0.0);
    for (int i = 0; i < N; ++i) u[i] = nd(gen);

    // crude lowpass: moving average
    const int W = 9;
    std::vector<double> u_lp(N, 0.0);
    for (int i = 0; i < N; ++i) {
        double s = 0.0;
        int cnt = 0;
        for (int j = -W; j <= W; ++j) {
            int k = i + j;
            if (k >= 0 && k < N) { s += u[k]; cnt++; }
        }
        u_lp[i] = s / (double)cnt;
    }
    u = u_lp;

    // simulate (force at DOF0, measure accel at DOF0)
    VectorXd x = VectorXd::Zero(6);
    std::vector<double> y(N, 0.0);

    auto fvec = [&](int i)->VectorXd{
        VectorXd f = VectorXd::Zero(3);
        f(0) = u[i];
        return f;
    };

    auto rhs = [&](const VectorXd& xk, const VectorXd& fk)->VectorXd{
        return A * xk + B * fk;
    };

    for (int i = 0; i < N; ++i) {
        VectorXd f = fvec(i);

        // RK4
        VectorXd k1 = rhs(x, f);
        VectorXd k2 = rhs(x + 0.5*dt*k1, f);
        VectorXd k3 = rhs(x + 0.5*dt*k2, f);
        VectorXd k4 = rhs(x + dt*k3, f);
        x = x + (dt/6.0)*(k1 + 2.0*k2 + 2.0*k3 + k4);

        // acceleration output (without direct term Minv*f for simplicity; add it)
        VectorXd acc = Ca * x + Minv * f;
        y[i] = acc(0);
    }

    // FRF estimate with H1 (single-segment)
    Spectra sp = spectra_single_segment(u, y, fs);
    const int Kbins = (int)sp.f.size();

    std::cout << "f_Hz, |H1|, coherence\n";
    for (int k = 1; k < Kbins && sp.f[k] <= 120.0; ++k) {
        std::complex<double> H1 = sp.Gyu[k] / (sp.Guu[k] + std::complex<double>(1e-30,0.0));
        double coh = std::norm(sp.Gyu[k]) / ( (std::real(sp.Guu[k])*std::real(sp.Gyy[k]) + 1e-30) );
        std::cout << sp.f[k] << ", " << std::abs(H1) << ", " << coh << "\n";
    }

    std::cout << "\nDone.\n";
    return 0;
}
      

This Java example shows a minimal complex-number DFT and FRF calculation. For realistic modal testing pipelines, you would use FFT libraries and Welch averaging.

Chapter13_Lesson5.java

/*
Chapter 13 - Vibrations and Multi-Degree-of-Freedom (MDOF) Systems
Lesson 5 - Intro to Experimental Modal Analysis and System Identification Concepts

File: Chapter13_Lesson5.java

This Java example demonstrates a minimal FRF (H1) estimation from measured input u[n]
and output y[n] using a direct DFT (O(N^2)). It is intentionally simple and
educational.

For realistic work:
- Use FFT libraries (e.g., JTransforms) and Welch averaging.
- Use linear algebra libraries for MIMO and state-space identification.

Compile & run:
  javac Chapter13_Lesson5.java
  java Chapter13_Lesson5
*/

import java.util.Random;

public class Chapter13_Lesson5 {

    static final double PI = Math.PI;

    static class Complex {
        double re, im;
        Complex(double re, double im) { this.re = re; this.im = im; }
        Complex add(Complex b) { return new Complex(this.re + b.re, this.im + b.im); }
        Complex mul(Complex b) {
            return new Complex(this.re*b.re - this.im*b.im, this.re*b.im + this.im*b.re);
        }
        Complex conj() { return new Complex(this.re, -this.im); }
        double abs() { return Math.hypot(re, im); }
        double norm() { return re*re + im*im; }
        Complex div(Complex b) {
            double den = b.re*b.re + b.im*b.im + 1e-30;
            return new Complex((this.re*b.re + this.im*b.im)/den, (this.im*b.re - this.re*b.im)/den);
        }
        static Complex expj(double ang) { return new Complex(Math.cos(ang), Math.sin(ang)); }
    }

    static Complex[] dft(double[] x) {
        int N = x.length;
        Complex[] X = new Complex[N];
        for (int k = 0; k < N; k++) {
            Complex s = new Complex(0.0, 0.0);
            for (int n = 0; n < N; n++) {
                double ang = -2.0 * PI * k * n / (double)N;
                Complex w = Complex.expj(ang);
                s = s.add(new Complex(x[n], 0.0).mul(w));
            }
            X[k] = s;
        }
        return X;
    }

    public static void main(String[] args) {
        double fs = 500.0;
        double T = 8.0;
        int N = (int)Math.floor(fs * T);

        // Example "measured" signals (for demo only):
        // u: band-limited noise; y: u passed through a lightly damped 2nd-order filter + noise.
        double[] u = new double[N];
        double[] y = new double[N];
        Random rng = new Random(7);

        for (int i = 0; i < N; i++) u[i] = rng.nextGaussian();

        // simple smoothing lowpass
        int W = 9;
        double[] uLP = new double[N];
        for (int i = 0; i < N; i++) {
            double s = 0.0; int cnt = 0;
            for (int j = -W; j <= W; j++) {
                int k = i + j;
                if (0 <= k && k < N) { s += u[k]; cnt++; }
            }
            uLP[i] = s / (double)cnt;
        }
        u = uLP;

        // 2nd-order resonant filter (biquad-like) to mimic a mode
        double fn = 18.0;  // Hz
        double zeta = 0.03;
        double w0 = 2.0 * PI * fn;
        double dt = 1.0 / fs;

        double x1 = 0.0, x2 = 0.0; // states: position and velocity
        for (int i = 0; i < N; i++) {
            double f = u[i];

            // x1dot = x2
            // x2dot = -2 zeta w0 x2 - w0^2 x1 + f
            double k1_1 = x2;
            double k1_2 = -2.0*zeta*w0*x2 - w0*w0*x1 + f;

            double x1m = x1 + 0.5*dt*k1_1;
            double x2m = x2 + 0.5*dt*k1_2;
            double k2_1 = x2m;
            double k2_2 = -2.0*zeta*w0*x2m - w0*w0*x1m + f;

            x1 += dt*k2_1;
            x2 += dt*k2_2;

            double acc = -2.0*zeta*w0*x2 - w0*w0*x1 + f;
            y[i] = acc + 0.02*rng.nextGaussian();
        }

        // DFT-based spectra (single segment)
        Complex[] U = dft(u);
        Complex[] Y = dft(y);
        int K = N/2 + 1;

        System.out.println("f_Hz, |H1|, coherence");
        for (int k = 1; k < K && (fs*k/(double)N) <= 120.0; k++) {
            Complex Guu = U[k].mul(U[k].conj()).mul(new Complex(1.0/N, 0.0));
            Complex Gyy = Y[k].mul(Y[k].conj()).mul(new Complex(1.0/N, 0.0));
            Complex Gyu = Y[k].mul(U[k].conj()).mul(new Complex(1.0/N, 0.0));

            Complex H1 = Gyu.div(Guu);
            double coh = Gyu.norm() / ( (Guu.re*Gyy.re) + 1e-30 );

            double f = fs*k/(double)N;
            System.out.println(f + ", " + H1.abs() + ", " + coh);
        }
        System.out.println("Done.");
    }
}
      

This script uses pwelch and cpsd to estimate spectra and compute the H1 FRF and coherence. It also programmatically generates a simple Simulink model containing a State-Space block driven by band-limited noise.

Chapter13_Lesson5.m

% Chapter 13 - Vibrations and Multi-Degree-of-Freedom (MDOF) Systems
% Lesson 5 - Intro to Experimental Modal Analysis and System Identification Concepts
%
% File: Chapter13_Lesson5.m
%
% This script:
% 1) Builds a 3-DOF mass-spring-damper model.
% 2) Simulates broadband excitation with lsim.
% 3) Estimates FRF using H1 = G_yu / G_uu via Welch averages (pwelch/cpsd).
% 4) Computes coherence, does rough peak-picking + half-power damping.
% 5) (Optional) Auto-builds a simple Simulink model containing the same state-space block.

clear; clc; close all;

%% 1) Model (3-DOF chain)
m1=1.0; m2=0.9; m3=0.8;
M = diag([m1 m2 m3]);

k1=800; k2=600; k3=500; k4=700;
K = [k1+k2  -k2      0;
     -k2    k2+k3   -k3;
      0     -k3     k3+k4];

c1=2.0; c2=1.8; c3=1.5; c4=2.2;
C = [c1+c2  -c2      0;
     -c2    c2+c3   -c3;
      0     -c3     c3+c4];

n = 3;
Minv = inv(M);
Z = zeros(n); I = eye(n);

A = [Z I; -Minv*K -Minv*C];
B = [zeros(n); Minv];    % inputs are forces at each DOF

% output: acceleration at DOF1 (index 1) => qdd = -Minv*K q - Minv*C qd + Minv*f
Ca = [-Minv*K -Minv*C];
Da = Minv;

sysa = ss(A, B, Ca, Da);

%% 2) Excitation + simulation
fs = 500;                 % Hz
T  = 40;                  % seconds
t  = (0:1/fs:T-1/fs)';

rng(7);
u = randn(length(t),1);
[b,a] = butter(4, 0.25);  % normalized cutoff
u = filtfilt(b,a,u);
u = u/std(u);

F = zeros(length(t), n);
F(:,1) = u; % force on DOF1

y = lsim(sysa, F, t);     % y is 3 columns: acceleration at each DOF
y1 = y(:,1) + 0.02*randn(size(y(:,1)));

%% 3) FRF estimation (H1) and coherence
nperseg = 4096;
nover   = 2048;
win     = hanning(nperseg);

[Guu,f] = pwelch(u, win, nover, nperseg, fs, 'onesided');
[Gyy,~] = pwelch(y1, win, nover, nperseg, fs, 'onesided');
[Gyu,~] = cpsd(y1, u, win, nover, nperseg, fs, 'onesided');

H1 = Gyu ./ (Guu + 1e-30);
coh = (abs(Gyu).^2) ./ (Guu.*Gyy + 1e-30);

figure; semilogy(f, abs(H1)); xlim([0 120]);
xlabel('Frequency [Hz]'); ylabel('|H1(jw)|'); title('Estimated FRF (H1)');

figure; plot(f, coh); ylim([0 1.05]); xlim([0 120]);
xlabel('Frequency [Hz]'); ylabel('Coherence'); title('Magnitude-squared coherence');

%% 4) Peak-picking + half-power damping (rough)
fmax = 80;
mask = (f > 0.5) & (f < fmax) & (coh > 0.6);
ff = f(mask);
Hm = abs(H1(mask));

[pks, locs] = findpeaks(Hm, ff, 'MinPeakProminence', prctile(Hm,70), 'MinPeakDistance', 2);
locs = locs(1:min(3, numel(locs)));

fprintf('Rough modal estimates (fn Hz, zeta, f1,f2):\n');
for i = 1:numel(locs)
    fn = locs(i);
    peak = interp1(ff, Hm, fn);
    target = peak/sqrt(2);

    % left
    idx0 = find(ff <= fn, 1, 'last');
    idx1 = idx0;
    while idx1 > 2 && Hm(idx1) > target
        idx1 = idx1 - 1;
    end
    f1 = interp1([Hm(idx1) Hm(idx1+1)], [ff(idx1) ff(idx1+1)], target);

    % right
    idx2 = idx0;
    while idx2 < numel(ff)-1 && Hm(idx2) > target
        idx2 = idx2 + 1;
    end
    f2 = interp1([Hm(idx2-1) Hm(idx2)], [ff(idx2-1) ff(idx2)], target);

    zeta = (f2 - f1)/(2*fn);
    fprintf('  fn=%7.3f Hz, zeta=%7.4f, f1=%7.3f, f2=%7.3f\n', fn, zeta, f1, f2);
end

%% 5) Optional: build a simple Simulink model programmatically
% This makes a model with:
% [Band-Limited White Noise] -> [State-Space] -> [To Workspace]
%
% You can comment this section if you prefer manual modeling.

mdl = 'Chapter13_Lesson5_Simulink';
if bdIsLoaded(mdl); close_system(mdl,0); end
if exist([mdl '.slx'], 'file'); delete([mdl '.slx']); end

new_system(mdl);
open_system(mdl);

add_block('simulink/Sources/Band-Limited White Noise', [mdl '/Noise']);
set_param([mdl '/Noise'], 'SampleTime', num2str(1/fs), 'NoisePower', '1', 'Seed', '7');

add_block('simulink/Continuous/State-Space', [mdl '/Plant']);
set_param([mdl '/Plant'], 'A', 'A', 'B', 'B(:,1)', 'C', 'Ca(1,:)', 'D', 'Da(1,1)');

add_block('simulink/Sinks/To Workspace', [mdl '/y_to_ws']);
set_param([mdl '/y_to_ws'], 'VariableName', 'y_sim', 'SaveFormat', 'Array');

add_line(mdl, 'Noise/1', 'Plant/1');
add_line(mdl, 'Plant/1', 'y_to_ws/1');

set_param(mdl, 'StopTime', num2str(T));
save_system(mdl);

disp(['Simulink model saved: ' mdl '.slx']);
      

The notebook below implements a small 2-DOF simulation and a simple (single-segment) H1 FRF estimate using Fourier transforms. (For variance reduction, extend it to segmentation and averaging, analogous to Welch’s method.)

Chapter13_Lesson5.nb

(* ::Package:: *)

(* Chapter 13 - Vibrations and Multi-Degree-of-Freedom (MDOF) Systems
   Lesson 5 - Intro to Experimental Modal Analysis and System Identification Concepts

   File: Chapter13_Lesson5.nb

   This notebook:
   1) Defines a 2-DOF mass-spring-damper model in state-space.
   2) Simulates broadband excitation and computes a simple H1 FRF estimate using Fourier transforms.
   3) Computes coherence.
*)

Notebook[{

Cell["Chapter 13 - Lesson 5: Experimental Modal Analysis (EMA) and System ID (concepts)", "Title"],

Cell[CellGroupData[{
Cell["1) Define a 2-DOF model", "Section"],
Cell[BoxData[
 RowBox[{
  RowBox[{"m1", "=", "1.0"}], ";",
  RowBox[{"m2", "=", "0.8"}], ";",
  RowBox[{"M", "=", "DiagonalMatrix", "[",
   RowBox[{"{",
    RowBox[{"m1", ",", "m2"}], "}"}], "]"}], ";", "\n",
  RowBox[{"k1", "=", "600"}], ";",
  RowBox[{"k2", "=", "500"}], ";",
  RowBox[{"k3", "=", "700"}], ";",
  RowBox[{"K", "=", 
   RowBox[{"{",
    RowBox[{
     RowBox[{"{",
      RowBox[{"k1", "+", "k2", ",",
       RowBox[{"-", "k2"}]}], "}"}], ",",
     RowBox[{"{",
      RowBox[{
       RowBox[{"-", "k2"}], ",", "k2", "+", "k3"}], "}"}]}], "}"}]}], ";", "\n",
  RowBox[{"c1", "=", "2.0"}], ";",
  RowBox[{"c2", "=", "1.6"}], ";",
  RowBox[{"c3", "=", "2.1"}], ";",
  RowBox[{"C", "=", 
   RowBox[{"{",
    RowBox[{
     RowBox[{"{",
      RowBox[{"c1", "+", "c2", ",",
       RowBox[{"-", "c2"}]}], "}"}], ",",
     RowBox[{"{",
      RowBox[{
       RowBox[{"-", "c2"}], ",", "c2", "+", "c3"}], "}"}]}], "}"}]}], ";", "\n",
  RowBox[{"Minv", "=", 
   RowBox[{"Inverse", "[", "M", "]"}]}], ";"}]], "Input"],

Cell[BoxData[
 RowBox[{
  RowBox[{"n", "=", "2"}], ";",
  RowBox[{"Z", "=", "ConstantArray", "[",
   RowBox[{"0.0", ",",
    RowBox[{"{",
     RowBox[{"n", ",", "n"}], "}"}]}], "]"}], ";",
  RowBox[{"I", "=", "IdentityMatrix", "[", "n", "]"}], ";", "\n",
  RowBox[{"A", "=", 
   RowBox[{"ArrayFlatten", "[",
    RowBox[{"{",
     RowBox[{
      RowBox[{"{",
       RowBox[{"Z", ",", "I"}], "}"}], ",",
      RowBox[{"{",
       RowBox[{
        RowBox[{"-", "Minv", ".", "K"}], ",",
        RowBox[{"-", "Minv", ".", "C"}]}], "}"}]}], "}"}], "]"}]}], ";", "\n",
  RowBox[{"B", "=", 
   RowBox[{"ArrayFlatten", "[",
    RowBox[{"{",
     RowBox[{
      RowBox[{"{", "Z", "}"}], ",",
      RowBox[{"{", "Minv", "}"}]}], "}"}], "]"}]}], ";"}]], "Input"],

Cell["We excite DOF1 by force u(t) and measure acceleration at DOF1.", "Text"],

Cell[BoxData[
 RowBox[{
  RowBox[{"Ca", "=", 
   RowBox[{"ArrayFlatten", "[",
    RowBox[{"{",
     RowBox[{"{",
      RowBox[{
       RowBox[{"-", "Minv", ".", "K"}], ",",
       RowBox[{"-", "Minv", ".", "C"}]}], "}"}], "}"}], "]"}]}], ";",
  RowBox[{"Da", "=", "Minv"}], ";"}]], "Input"]
}, Open  ]],

Cell[CellGroupData[{
Cell["2) Simulate broadband excitation", "Section"],
Cell[BoxData[
 RowBox[{
  RowBox[{"fs", "=", "500.0"}], ";",
  RowBox[{"dt", "=",
   RowBox[{"1.0", "/", "fs"}]}], ";",
  RowBox[{"T", "=", "20.0"}], ";",
  RowBox[{"tgrid", "=",
   RowBox[{"Range", "[",
    RowBox[{"0.0", ",",
     RowBox[{"T", "-", "dt"}], ",", "dt"}], "]"}]}], ";",
  RowBox[{"N", "=",
   RowBox[{"Length", "[", "tgrid", "]"}]}], ";"}]], "Input"],

Cell[BoxData[
 RowBox[{
  RowBox[{"SeedRandom", "[", "7", "]"}], ";",
  RowBox[{"u", "=",
   RowBox[{"RandomVariate", "[",
    RowBox[{
     RowBox[{"NormalDistribution", "[",
      RowBox[{"0", ",", "1"}], "]"}], ",", "N"}], "]"}]}], ";",
  RowBox[{"u", "=",
   RowBox[{"MovingAverage", "[",
    RowBox[{"u", ",", "9"}], "]"}]}], ";",
  RowBox[{"u", "=",
   RowBox[{"u", "/",
    RowBox[{"StandardDeviation", "[", "u", "]"}]}]}], ";"}]], "Input"],

Cell[BoxData[
 RowBox[{
  RowBox[{"x", "=",
   RowBox[{"ConstantArray", "[",
    RowBox[{"0.0", ",", "4"}], "]"}]}], ";",
  RowBox[{"y", "=",
   RowBox[{"ConstantArray", "[",
    RowBox[{"0.0", ",", "N"}], "]"}]}], ";"}]], "Input"],

Cell[BoxData[
 RowBox[{
  RowBox[{"rhs", "[",
   RowBox[{"x_", ",", "f_"}], "]"}], ":=",
  RowBox[{"A", ".", "x", "+",
   RowBox[{"B", ".", "f"}]}]}]], "Input"],

Cell[BoxData[
 RowBox[{
  RowBox[{"Do", "[",
   RowBox[{
    RowBox[{
     RowBox[{"f", "=",
      RowBox[{"{",
       RowBox[{
        RowBox[{"u", "[",
         RowBox[{"[", "i", "]"}], "]"}], ",", "0.0"}], "}"}]}], ";",
     RowBox[{"k1", "=",
      RowBox[{"rhs", "[",
       RowBox[{"x", ",", "f"}], "]"}]}], ";",
     RowBox[{"k2", "=",
      RowBox[{"rhs", "[",
       RowBox[{
        RowBox[{"x", "+",
         RowBox[{"0.5", " ", "dt", " ", "k1"}]}], ",", "f"}], "]"}]}], ";",
     RowBox[{"k3", "=",
      RowBox[{"rhs", "[",
       RowBox[{
        RowBox[{"x", "+",
         RowBox[{"0.5", " ", "dt", " ", "k2"}]}], ",", "f"}], "]"}]}], ";",
     RowBox[{"k4", "=",
      RowBox[{"rhs", "[",
       RowBox[{
        RowBox[{"x", "+",
         RowBox[{"dt", " ", "k3"}]}], ",", "f"}], "]"}]}], ";",
     RowBox[{"x", "=",
      RowBox[{"x", "+",
       RowBox[{
        RowBox[{"dt", "/", "6.0"}],
        RowBox[{"(",
         RowBox[{"k1", "+",
          RowBox[{"2", "k2"}], "+",
          RowBox[{"2", "k3"}], "+", "k4"}], ")"}]}]}]}], ";",
     RowBox[{"acc", "=",
      RowBox[{
       RowBox[{"Ca", ".", "x"}], "+",
       RowBox[{"Da", ".", "f"}]}]}], ";",
     RowBox[{"y", "[",
      RowBox[{"[", "i", "]"}], "]"}], "=",
     RowBox[{"acc", "[",
      RowBox[{"[", "1", "]"}], "]"}]}], ",",
    RowBox[{"{",
     RowBox[{"i", ",", "1", ",", "N"}], "}"}]}], "]"}], ";"}]], "Input"]
}, Open  ]],

Cell[CellGroupData[{
Cell["3) Simple H1 FRF estimate via Fourier transforms", "Section"],
Cell["We use: H1(w) = G_yu(w) / G_uu(w) and coherence gamma^2(w) = |G_yu|^2/(G_uu G_yy).", "Text"],

Cell[BoxData[
 RowBox[{
  RowBox[{"U", "=",
   RowBox[{"Fourier", "[", "u", "]"}]}], ";",
  RowBox[{"Y", "=",
   RowBox[{"Fourier", "[", "y", "]"}]}], ";",
  RowBox[{"Guu", "=",
   RowBox[{
    RowBox[{"U", " ",
     RowBox[{"Conjugate", "[", "U", "]"}]}], "/", "N"}]}], ";",
  RowBox[{"Gyy", "=",
   RowBox[{
    RowBox[{"Y", " ",
     RowBox[{"Conjugate", "[", "Y", "]"}]}], "/", "N"}]}], ";",
  RowBox[{"Gyu", "=",
   RowBox[{
    RowBox[{"Y", " ",
     RowBox[{"Conjugate", "[", "U", "]"}]}], "/", "N"}]}], ";",
  RowBox[{"H1", "=",
   RowBox[{"Gyu", "/",
    RowBox[{"(",
     RowBox[{"Guu", "+",
      RowBox[{"10", "^",
       RowBox[{"-", "30"}]}]}], ")"}]}]}], ";",
  RowBox[{"coh", "=",
   RowBox[{
    RowBox[{"Abs", "[", "Gyu", "]"}], "^", "2", "/",
    RowBox[{"(",
     RowBox[{"Guu", " ", "Gyy", "+",
      RowBox[{"10", "^",
       RowBox[{"-", "30"}]}]}], ")"}]}]}], ";"}]], "Input"],

Cell[BoxData[
 RowBox[{
  RowBox[{"K", "=",
   RowBox[{"Floor", "[",
    RowBox[{"N", "/", "2"}], "]"}]}], ";",
  RowBox[{"freq", "=",
   RowBox[{"Table", "[",
    RowBox[{
     RowBox[{"fs", " ",
      RowBox[{"(",
       RowBox[{"k", "-", "1"}], ")"}], "/", "N"}], ",",
     RowBox[{"{",
      RowBox[{"k", ",", "1", ",",
       RowBox[{"K", "+", "1"}]}], "}"}]}], "]"}]}], ";"}]], "Input"],

Cell[BoxData[
 RowBox[{"ListLinePlot", "[",
  RowBox[{
   RowBox[{"Transpose", "[",
    RowBox[{"{",
     RowBox[{"freq", ",",
      RowBox[{"Abs", "[",
       RowBox[{"H1", "[",
        RowBox[{"[",
         RowBox[{"1", ";;",
          RowBox[{"K", "+", "1"}]}], "]"}], "]"}], "]"}]}], "}"}], "]"}], ",",
   RowBox[{"ScalingFunctions", "->",
    RowBox[{"{",
     RowBox[{"None", ",", "\"Log\""}], "}"}]}], ",",
   RowBox[{"PlotRange", "->",
    RowBox[{"{",
     RowBox[{
      RowBox[{"{",
       RowBox[{"0", ",", "120"}], "}"}], ",", "All"}], "}"}]}], ",",
   RowBox[{"AxesLabel", "->",
    RowBox[{"{",
     RowBox[{"\"Hz\"", ",", "\"|H1|\""}], "}"}]}], ",",
   RowBox[{"PlotLabel", "->", "\"FRF magnitude (H1)\""}]}], "]"}]], "Input"],

Cell[BoxData[
 RowBox[{"ListLinePlot", "[",
  RowBox[{
   RowBox[{"Transpose", "[",
    RowBox[{"{",
     RowBox[{"freq", ",",
      RowBox[{"coh", "[",
       RowBox[{"[",
        RowBox[{"1", ";;",
         RowBox[{"K", "+", "1"}]}], "]"}], "]"}]}], "}"}], "]"}], ",",
   RowBox[{"PlotRange", "->",
    RowBox[{"{",
     RowBox[{
      RowBox[{"{",
       RowBox[{"0", ",", "120"}], "}"}], ",",
      RowBox[{"{",
       RowBox[{"0", ",", "1.05"}], "}"}]}], "}"}]}], ",",
   RowBox[{"AxesLabel", "->",
    RowBox[{"{",
     RowBox[{"\"Hz\"", ",", "\"coherence\""}], "}"}]}], ",",
   RowBox[{"PlotLabel", "->", "\"Magnitude-squared coherence\""}]}], "]"}]], "Input"]
}, Open  ]]

}]
      

Problem 1 (Derive the displacement FRF matrix): Starting from \( \mathbf{M}\ddot{\mathbf{q}}+\mathbf{C}\dot{\mathbf{q}}+\mathbf{K}\mathbf{q}=\mathbf{B}_f\mathbf{f} \), assume a steady-state harmonic input at frequency \( \omega \) and derive the frequency response matrix mapping \( \mathbf{F}(j\omega) \) to \( \mathbf{Q}(j\omega) \).

Solution: Let \( \mathbf{q}(t)=\Re\{\mathbf{Q}e^{j\omega t}\} \) and \( \mathbf{f}(t)=\Re\{\mathbf{F}e^{j\omega t}\} \). Then \( \dot{\mathbf{q}} \mapsto j\omega \mathbf{Q} \) and \( \ddot{\mathbf{q}} \mapsto -\omega^2 \mathbf{Q} \). Substitution yields:

\[ \big(\mathbf{K}-\omega^2\mathbf{M}+j\omega \mathbf{C}\big)\mathbf{Q}=\mathbf{B}_f\mathbf{F} \;\Rightarrow\; \mathbf{Q}(j\omega)=\underbrace{\big(\mathbf{K}-\omega^2\mathbf{M}+j\omega \mathbf{C}\big)^{-1}\mathbf{B}_f}_{\mathbf{H}(j\omega)}\mathbf{F}(j\omega). \]

Problem 2 (Unbiasedness of \(H_1\)): Prove that if \(y=y_0+v\), \(y_0=Hu\), and \(v\) is uncorrelated with \(u\), then \(G_{yu}=HG_{uu}\), hence \(H_1=G_{yu}/G_{uu}=H\).

Solution: By linearity of cross spectra, \( G_{yu}=G_{(y_0+v)u}=G_{y_0u}+G_{vu} \). Since \(y_0=Hu\), the cross-spectrum factorizes: \( G_{y_0u}=H\,G_{uu} \). Uncorrelatedness gives \(G_{vu}=0\). Therefore \(G_{yu}=H G_{uu}\) and \(H_1=H\).

Problem 3 (Coherence bounds): Show that the magnitude-squared coherence satisfies \(0 \le \gamma^2(\omega) \le 1\).

Solution: Consider Fourier coefficients \(U(\omega)\) and \(Y(\omega)\) as random variables. Define an inner product \( \langle Y, U\rangle = E[YU^\ast] \), so that \(G_{yu}=\langle Y,U\rangle\), \(G_{uu}=\langle U,U\rangle\), and \(G_{yy}=\langle Y,Y\rangle\). Cauchy–Schwarz gives:

\[ |G_{yu}|^2 = |\langle Y,U\rangle|^2 \le \langle Y,Y\rangle\langle U,U\rangle = G_{yy}G_{uu}. \]

Divide both sides by \(G_{uu}G_{yy}\) (positive) to obtain \( \gamma^2(\omega) \le 1 \); nonnegativity is obvious.

Problem 4 (Half-power bandwidth formula): For the SDOF FRF \( H(\omega)=\frac{1}{k}\frac{1}{1-r^2 + j2\zeta r} \) with \( r=\omega/\omega_n \), derive (for small \(\zeta\)) the approximation \( \zeta \approx (\omega_2-\omega_1)/(2\omega_n) \) where \(\omega_1,\omega_2\) are half-power frequencies.

Solution: Compute the magnitude squared:

\[ |H(\omega)|^2 = \frac{1}{k^2}\cdot\frac{1}{(1-r^2)^2 + (2\zeta r)^2}. \]

The peak occurs near \(r\approx 1\), where \(|H|^2_{\max}\approx \frac{1}{k^2}\frac{1}{(2\zeta)^2}\). Half-power means \(|H|^2 = |H|^2_{\max}/2\), so:

\[ (1-r^2)^2 + (2\zeta r)^2 \approx 2(2\zeta)^2. \]

For light damping, \(r\approx 1\) at the half-power points, so set \(r\approx 1\) in the \(2\zeta r\) term and solve \((1-r^2)^2 \approx (2\zeta)^2\). Taking square roots gives \(1-r^2 \approx \pm 2\zeta\). Linearize \(r^2\approx 1+2(r-1)\) to obtain \(r-1 \approx \pm \zeta\), hence \(\omega_{1,2}=\omega_n(1\mp \zeta)\) and \(\omega_2-\omega_1\approx 2\zeta\omega_n\).

Problem 5 (Normal equations for a linear-in-parameters predictor): Suppose a predictor is linear in parameters: \( \hat{y}[k\mid \boldsymbol{\theta}]=\boldsymbol{\varphi}[k]^{\mathsf{T}}\boldsymbol{\theta} \), and define the prediction error \( \varepsilon[k;\boldsymbol{\theta}]=y[k]-\boldsymbol{\varphi}[k]^{\mathsf{T}}\boldsymbol{\theta} \). Derive the normal equations for minimizing \( \sum_k \varepsilon[k;\boldsymbol{\theta}]^2 \).

Solution: Let \(J(\boldsymbol{\theta})=\sum_k (y[k]-\boldsymbol{\varphi}[k]^{\mathsf{T}}\boldsymbol{\theta})^2\). Differentiate and set gradient to zero:

\[ \frac{\partial J}{\partial \boldsymbol{\theta}} = -2\sum_k \boldsymbol{\varphi}[k]\big(y[k]-\boldsymbol{\varphi}[k]^{\mathsf{T}}\boldsymbol{\theta}\big)=\mathbf{0} \;\Rightarrow\; \Big(\sum_k \boldsymbol{\varphi}[k]\boldsymbol{\varphi}[k]^{\mathsf{T}}\Big)\boldsymbol{\theta} = \sum_k \boldsymbol{\varphi}[k]y[k]. \]

This is the least-squares solution \( \hat{\boldsymbol{\theta}}=(\Phi^{\mathsf{T}}\Phi)^{-1}\Phi^{\mathsf{T}}y \) when \( \Phi^{\mathsf{T}}\Phi \) is nonsingular.

We formalized the experimental pathway from measured input–output vibration data to modal parameters, derived core FRF estimators (H1/H2) and coherence, linked modal peaks to poles and damping, and introduced the system-identification mindset: choose a model class and estimate parameters by prediction-error / likelihood principles, then validate aggressively. These concepts form the bridge between physical vibration models and data-driven model building used in modern control engineering.