Operational Modal Analysis

Operational Modal Analysis (OMA) is a technique used to extract modal parameters of a structure based solely on its response measurements, without requiring knowledge of the input excitation forces. This is particularly useful in scenarios where it is difficult or impossible to measure the input forces, such as in ambient vibration testing of large structures.

From an algorithmic perspective, OMA can be performed using various methods, including EMA-based algorithms or specific OMA techniques such as Stochastic Subspace Identification techniques. In StructuralVibration.jl, the following methods are implemented for OMA:

1 Sdof method

The Frequency-Spatial Domain Decomposition (FSDD) method is a single-degree-of-freedom (Sdof) technique used in Operational Modal Analysis (OMA) to extract modal parameters from output-only data1. The FSDD method relies on the fact that in the vicinity of a resonance frequency \(\omega_k\), the cross-spectral density (CSD) matrix of the output responses can be approximated by a rank-1 matrix, which corresponds to the dominant mode shape of the structure at that frequency. In particular, this means that the CSD matrix \(\mathbf{G}_{yy}(\omega_k)\) can be expressed as: \[ \mathbf{G}_{yy}(\omega_k) \approx \mathbf{u}_1(\omega_k)\, s_1(\omega_k) \mathbf{u}_1(\omega_k)^H, \] where \(\mathbf{u}_1(\omega_k)\) is the first singular vector of the CSD matrix, and \(s_1(\omega_k)\) is the corresponding singular value.

By noting that the CSD matrix near the resonance frequency can also be expressed in terms of the mode shape \(\boldsymbol{\Phi}\) and the power spectral density of the modal coordinate \({G}_{qq}(\omega)\) as: \[ \mathbf{G}_{yy}(\omega_k) \approx \boldsymbol{\Phi}\boldsymbol{G}_{qq}(\omega_k) \boldsymbol{\Phi}^H, \] where \(\boldsymbol{G}_{qq}(\omega)\) is dominated by the contribution of the \(k\)-th mode, we can establish a relationship between the singular vectors and the mode shapes. More specifically, the mode shape \(\boldsymbol{\Phi}_k\) at frequency \(\omega_k\) is given by the first singular vector: \[ \boldsymbol{\Phi}_k \approx \mathbf{u}_1(\omega_k). \]

Once the mode shapes are extracted, the FSDD utilizes the unitary property of the singular vectors to derive an enhanced output PSD \(\hat{G}_{yy}(\omega)\), which is a real quantity, in the vicinity of the resonance frequency, such that: \[ \hat{G}_{yy}(\omega) = \mathbf{u}_1(\omega_k)^H\, \mathbf{G}_{yy}(\omega)\, \mathbf{u}_1(\omega_k) \approx \frac{d_k}{\sigma_k^2 + (\omega - \Omega_k)^2}, \] where \(d_k\) is a constant related to the mode shape normalization, \(\sigma_k = -\xi_k\omega_k\) , and \(\Omega_k\) is the damped natural angular frequency of the \(k\)-th mode. By fitting a Sdof model to the enhanced PSD \(\hat{G}_{yy}(\omega)\) around the resonance frequency \(\omega_k\) using a least-squares approach, the modal parameters (natural frequency and damping ratio) can be estimated.

Note

The FSDD method is particularly effective for structures with well-separated modes and low damping, as it relies on the dominance of a single mode in the CSD matrix near resonance frequencies. However, it may face challenges when dealing with closely spaced modes or high damping scenarios, where the rank-1 approximation may not hold true.

To perform a robust FSDD, StructuralVibration.jl implements two additional steps during the estimation process:

  1. Because the signals are often contaminated with noise, the singular vector \(\mathbf{u}_1(\omega_k)\) is estimated by averaging the first left-singlular vectors within a frequency band around each resonance frequency. This aims to reduce the effects of noise on the mode shape estimates. During this process, a MAC threshold is used to select only the singular vectors that are sufficiently similar to the reference singular vector at the peak frequency. Said differently, only the singular vectors, for which the first singular value is approximately constant within the frequency band, are considered for the averaging.

    In doing so, we define the so-called Sdof-bell corresponding to the frequencies around each peak where the singular vectors are considered for averaging. It results that the left-singular vector \(\mathbf{u}_1(\omega_k)\) is estimated as: \[ \hat{\mathbf{u}}_1(\omega_k) = \frac{\sum_{i=1}^{N_k} w_i\, \mathbf{u}_1(\omega_i)}{\sum_{i=1}^{N_k} w_i}, \]

    where \(N_k\) is the number of selected candidates within the bell and \(w_i\) are the weights assigned to each candidate. StructuralVibration.jl proposes two types of weights:

    • Uniform weights: \(w_i = 1\) corresponding to a linear average of the selected candidates.

    • Singular value-based weights: \(w_i = s_1(\omega_i)\). This weighting scheme gives more importance to candidates that are more likely to represent the true mode shape and less importance to those that may be influenced by noise or other factors.

    Knowing the Sdof-bell also allows to define the frequency range around each resonance frequency where the Sdof model is fitted to the enhanced PSD. To do so, the fitting equation is expressed as: \[ \begin{bmatrix} \hat{G}_{yy}(\omega) & -2\omega \hat{G}_{yy}(\omega) & -1 \\ \end{bmatrix} \begin{bmatrix} \omega_k^2 \\ \Omega_k \\ d_k \end{bmatrix} = -\omega^2 \hat{G}_{yy}(\omega) \]

    By applying this equation over all frequencies \(\omega\) within the Sdof-bell, a least-squares solution can be obtained for the unknowns \(\omega_k\), \(\Omega_k\), and \(d_k\). From these estimates, the poles \(p_k\) can be computed as: \[ p_k = -\sqrt{\omega_k^2 - \Omega_k^2} + j \Omega_k, \]

  2. For lightly damped structures, the least-squares fitting may failed to provide a physical solution due to numerical inaccuracies. In this case, one typically observes that \(\Omega_k > \omega_k\), leading to imaginary values for the damping ratio. To overcome this issue, StructuralVibration.jl uses the estimation inspired by the Enhanced Frequency Domain Decomposition (EFDD) method2.

    The estimation consists in taking the inverse Fourier transform of the enhanced PSD \(\hat{G}_{yy}(\omega)\) to obtain the corresponding Autocorrelation function \(R(\tau)\) of the Sdof system: \[ R(\tau) \approx C\, e^{-\sigma_k \tau} \cos(\Omega_k \tau). \]

    Then, we compute the envelope \(E(\tau)\) of the signal using the Hilbert transform, which gives: \[ E(\tau) = C\, e^{-\sigma_k \tau}. \] By fitting a linear model to the logarithm of the envelope, we can estimate \(\sigma_k\) from a least-squares regression.

    Assuming that the damped natural frequency \(\Omega_k\) was properly estimated in the initial estimation, the pole \(p_k\) can then be estimated as: \[ p_k = -\sigma_k + j \Omega_k. \]

2 Mdof methods

2.1 EMA-based methods for OMA

2.1.1 General principle

The EMA-based methods for OMA in StructuralVibration.jl leverage the existing Mdof EMA algorithms to extract modal parameters from response-only data. For applying these methods, one has to replace the Frequency Response Function by the so-called positive half-spectrum of the response signals.

To obtain the half-spectrum, we start from the full cross-spectral density (CSD) matrix between the output responses \(\mathbf{y}\), which is given by: \[ \mathbf{G}_{yy}(\omega) = \mathbf{H}(\omega) \mathbf{G}_{ff}(\omega) \mathbf{H}^H(\omega), \] where \(\mathbf{H}(\omega)\) is the Frequency Response Function matrix of the system, \(\mathbf{G}_{ff}(\omega)\) is the cross-spectral density matrix of the input forces, and \(^H\) denotes the Hermitian transpose.

Generally, OMA relies on the assumption that the input excitation are random white noise processes, implying that \(\mathbf{G}_{ff}(\omega) = \mathbf{I}\). Under this assumption, a component of the cross-spectral density matrix \(G_{pq}(\omega)\) can be expressed under the poles-residues form as34: \[ G_{pq}(\omega) = \sum_{i=1}^{n} \left( \frac{{}_iR_{pq}}{j\omega - \lambda_i} + \frac{{}_iR_{pq}^\ast}{j\omega - \lambda_i^\ast} + \frac{{}_iS_{pq}}{-j\omega - \lambda_i} + \frac{{}_iS_{pq}^\ast}{-j\omega - \lambda_i^\ast} \right), \] where \(\lambda_i\) and (\({}_iR_{pq}\), \({}_iS_{pq}\)) are the system poles and residues, respectively.

From this definition, the half-spectrum \(G_{pq}^+(\omega)\) is defined as the sum of the first two terms in the above expression5: \[ G_{pq}^+(\omega) = \sum_{i=1}^{n} \left( \frac{{}_iR_{pq}}{j\omega - \lambda_i} + \frac{{}_iR_{pq}^\ast}{j\omega - \lambda_i^\ast} \right). \]

In StructuralVibration.jl, the half-spectrum can be computed using the half_csd function, which can proceed in two different ways. These approaches can be summarized as follows:

  • First approach - Post-processing of the full CSD matrix:

    1. Compute the correlation matrix \(\mathbf{R}_{yy}(\tau)\) of the response signals by applying an inverse Fourier transform to the CSD matrix \(\mathbf{G}_{yy}(\omega)\).

    2. Keep only the positive time lags of the correlation matrix \(\mathbf{R}_{yy}^+(\tau)\).

    3. Multiply \(\mathbf{R}_{yy}^+(\tau)\) by a flat triangular window and zero-pad the negative time lags of the correlation matrix to obtain a full matrix of the same size as the original correlation matrix (See Ref. [3]). The resulting signal is noted \(\tilde{\mathbf{R}}_{yy}^+(\tau)\).

    4. Compute the half-spectrum \(\mathbf{G}_{yy}^+(\omega)\) by applying a Fourier transform to \(\tilde{\mathbf{R}}_{yy}^+(\tau)\).

  • Second approach - Direct computation from time-domain signals:

    1. Compute the discrete directe correlation matrix \(\mathbf{R}^+_k\) of the response signals directly from the time-domain data, considering only the positive time lags, using the following formula: \[ \mathbf{R}^+_k = \frac{1}{N - k} \sum_{n=1}^{N-k} \mathbf{y}_n \mathbf{y}_{n+k}^T \] where \(N\) is the total number of samples.

    2. Multiply \(\mathbf{R}_k^+\) by a flat triangular window and zero-pad the signal to obtain \(\tilde{\mathbf{R}}_k^+\).

    3. Compute the half-spectrum \(\mathbf{G}_{yy}^+(\omega)\) by applying a Fourier transform to \(\tilde{\mathbf{R}}_k^+\).

Note

If the length of the time signals are longer than the block size used for the FFT computation, the half-spectrum is computed for each block and then averaged over all blocks.

2.1.3 Half-spectrum reconstruction

From the identified modal parameters (poles and residues), it is possible to reconstruct the half_spectrum of the system under study using the half_csd_reconstruction function. This function implements the following relation for a number \(M\) of identified modes:

\[ G_{pq}^+(\omega) = \sum_{i = 1}^{M} \left[\frac{{}_i R_{pq}}{(j\omega - \lambda_i)} + \frac{{}_i R_{pq}^\ast}{(j\omega - \lambda_i^\ast)}\right], \]

However, as for the FRF reconstruction in EMA, the reconstructed half-spectrum may not perfectly match the original half-spectrum due to modeling errors and noise in the data. To To account for this, it is common to include lower and upper residuals in the FRF reconstruction. This means that the admittance matrix can be expressed as: \[ G_{pq}^+(\omega) = \sum_{i = 1}^{M} \left[\frac{{}_i R_{pq}}{(j\omega - \lambda_i)} + \frac{{}_i R_{pq}^\ast}{(j\omega - \lambda_i^\ast)}\right] + \frac{L_{pq}}{j\omega} + j\omega\, U_{pq}, \]

where \(L_{pq}\) and \(U_{pq}\) are the lower and upper residuals, respectively. These residuals can be computed using the compute_residuals function.

Note

As indicated by Peeters and Van der Auweraer in Ref. [2], the lower and upper residuals are independent of the output quantity. More precisely, this means that the residuals terms are always of the form: \[ \text{Res}_{pq} = \frac{L_{pq}}{j\omega} + j\omega\, U_{pq}. \]

2.2 Stochastic Subspace Identification methods

Stochastic Subspace Identification (SSI) methods are a class of techniques used in Operational Modal Analysis (OMA) to extract modal parameters from output-only data. These methods are based on the concept of subspace identification, which involves constructing a state-space model of the system from the measured response data. From a mathematical standpoint, the goal is to identify the system matrices \((\mathbf{A}, \mathbf{C})\) of the state-space representation: \[ \begin{aligned} \dot{\mathbf{x}}_{k+1} &= \mathbf{A}\, \mathbf{x}_k + \mathbf{w}_k \\ \mathbf{y}_k &= \mathbf{C}\, \mathbf{x}_k + \mathbf{v}_k \end{aligned} \] where \(\mathbf{x}_k\) is the state vector, \(\mathbf{y}_k\) is the output vector, and \(\mathbf{w}_k\) and \(\mathbf{v}_k\) are process and measurement noise, respectively. Here, the system matrix \(\mathbf{A}\) is the discrete-time state transition matrix that contains information about the system dynamics. If a ZOH discretization is considered, the relationship between the continuous-time system matrix \(\mathbf{A}_c\) and the discrete-time system matrix \(\mathbf{A}\) is given by: \[ \mathbf{A} = e^{\mathbf{A}_c T_s} \] where \(T_s\) is the sampling period. The eigenvalues of \(\mathbf{A}\) can be used to extract the continuous-time poles \(\lambda_i\) of the system using the following relation: \[ \lambda_i = \frac{\ln(\mu_i)}{T_s} \] where \(\mu_i\) are the eigenvalues of \(\mathbf{A}\).

The mode shapes \(\boldsymbol{\phi}_i\) can be obtained from the output matrix \(\mathbf{C}\) and the eigenvectors of \(\mathbf{A}\) using the following relation: \[ \boldsymbol{\phi}_i = \mathbf{C}\, \boldsymbol{\psi}_i \] where \(\boldsymbol{\psi}_i\) are the eigenvectors of \(\mathbf{A}\).

2.2.1 Covariance-driven Stochastic Subspace Identification

Covariance-driven Stochastic Subspace Identification (Cov-SSI) is a method that utilizes the covariance of the output data to identify the system matrices. The key steps involved in Cov-SSI are as follows:

  1. Compute the covariance matrices \(\mathbf{R}_i\) of the output data at different time lags \(i\): \[ \mathbf{R}_i = \frac{1}{N - i} \sum_{n = 0}^{N-i} \mathbf{y}_n\, \mathbf{y}_{n+i}^T \]
Note

If a set of reference signals \(\mathbf{y}_{ref}\) is provided, the cross-covariance matrices between the output signals and the reference signals are computed instead: \[ \mathbf{R}_i = \frac{1}{N - i} \sum_{n=0}^{N-i} \mathbf{y}_n\, {\mathbf{y}_{n+i}^\text{ref}}^T \]

  1. Construct a block Toeplitz matrix from the covariance matrices:

\[ \mathbf{T} = \begin{bmatrix} \mathbf{R}_i & \mathbf{R}_{i-1} & \cdots & \mathbf{R}_1 \\ \mathbf{R}_{i+1} & \mathbf{R}_i & \ddots & \mathbf{R}_2 \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{R}_{2i-1} & \mathbf{R}_{2i-2} & \cdots & \mathbf{R}_i \end{bmatrix} \]

  1. Compute the observability matrix \(\mathbf{O}\) from the SVD of \(\mathbf{T}\): \[ \mathbf{O} = \begin{bmatrix} \mathbf{C} \\ \mathbf{C} \mathbf{A} \\ \mathbf{C} \mathbf{A}^2 \\ \vdots \\ \mathbf{C} \mathbf{A}^{i-1} \end{bmatrix} = \mathbf{U}_n\, \mathbf{S}_n^{1/2}, \] where \(\mathbf{U}_n\) and \(\mathbf{S}_n\) are the matrices containing the first \(n\) left singular vectors and singular values, respectively.

  2. Extract the system matrices \(\mathbf{A}\) and \(\mathbf{C}\) from the Observability matrix

    • The output matrix \(\mathbf{C}\) is obtained from the first block row of \(\mathbf{O}\).
    • The system matrix \(\mathbf{A}\) is obtained by solving the following equation: \[ \mathbf{O}_\text{up}\, \mathbf{A} = \mathbf{O}_\text{down} \] where \(\mathbf{O}_\text{up}\) and \(\mathbf{O}_\text{down}\) are derived from the observability matrix removing the last and first block rows, respectively.

2.2.2 Data-driven Stochastic Subspace Identification

Data-driven Stochastic Subspace Identification (Data-SSI) is a method that directly utilizes the output data to identify the system matrices without explicitly computing the covariance matrices. The key steps involved in Data-SSI are as follows:

  1. Construct a block Hankel matrix from the output data: \[ \mathbf{H} = \frac{1}{\sqrt{j}}\begin{bmatrix} \mathbf{y}_0 & \mathbf{y}_1 & \cdots & \mathbf{y}_{j-1} \\ \mathbf{y}_1 & \mathbf{y}_2 & \cdots & \mathbf{y}_j \\ \cdots & \cdots & \cdots & \cdots \\ \mathbf{y}_{i - 1} & \mathbf{y}_i & \cdots & \mathbf{y}_{i+j-2} \\ \hline \mathbf{y}_i & \mathbf{y}_{i+1} & \cdots & \mathbf{y}_{i+j-1} \\ \mathbf{y}_{i+1} & \mathbf{y}_{i+2} & \cdots & \mathbf{y}_{i+j} \\ \cdots & \cdots & \cdots & \cdots \\ \mathbf{y}_{2i - 1} & \mathbf{y}_{2i} & \cdots & \mathbf{y}_{2i+j-2} \end{bmatrix} = \begin{bmatrix} \mathbf{Y}_p \\ \hline \mathbf{Y}_f \end{bmatrix} \]

    where \(i\) is the number of block rows and \(j\) is the number of columns of the Hankel matrix and \(\mathbf{Y}_p\) and \(\mathbf{Y}_f\) are the past and future output data matrices. respectively.

Note

If a set of reference signals \(\mathbf{y}_{ref}\) is provided, the past data matrix \(\mathbf{Y}_p\) is constructed using the reference signals instead of the output signals.

  1. Construct the projection of the future output data onto the past output data: \[ \mathbf{P} = \mathbf{Y}_f\, \mathbf{Y}_p^T\, (\mathbf{Y}_p\, \mathbf{Y}_p^T)^+\, \mathbf{Y}_p \]

  2. Compute the observability matrix \(\mathbf{O}\) from the SVD of \(\mathbf{P}\).

  3. Extract the system matrices \(\mathbf{A}\) and \(\mathbf{C}\) from the observability matrix as described in the Cov-SSI method.

3 API

Data types

OMAProblem(Gyy, freq; frange, type_data)
OMAProblem(y, yref, t, fs, bs; frange, type_data, win, overlap)
OMAProblem(y, t, fs, bs; frange, type_data, win, overlap)

Data structure defining the inputs for Operational Modal Analysis (OMA) methods.

Constructor parameters

  • Gyy::Array{Complex, 3}: Cross-spectral density matrix of size (no, nref, nf),

  • freq::AbstractArray{Real}: Frequency vector corresponding to the CSD matrix

  • frange::AbstractVector{Real}: Frequency range to consider for analysis (default: full range)

Alternative constructor parameters

  • y::AbstractMatrix{Real}: Matrix of measured outputs (no x nt)

  • yref::AbstractMatrix{Real}: Matrix of reference outputs (nref x nt)

  • t::AbstractArray{Real}: Time vector corresponding to the measurements

  • fs::Real: Sampling frequency (Hz)

  • bs::Int: Block size for CSD estimation

  • frange::AbstractVector{Real}: Frequency range to consider for analysis (default: [0., fs/2.56])

  • win::Function: Window function for CSD estimation (default: hanning)

  • overlap::Real: Overlap ratio for CSD estimation (default: 0.5)

Fields

  • y::AbstractMatrix{Real}: Matrix of measured outputs (no x nt)

  • yref::AbstractMatrix{Real}: Matrix of reference outputs (nref x nt)

  • t::AbstractArray{Real}: Time vector corresponding to the measurements

  • fullspec::Array{Complex, 3}: Full cross-spectral density matrix

  • halfspec::Array{Complex, 3}: Half spectral density matrix

  • freq::AbstractArray{Real}: Frequency vector corresponding to the half-spectrum

Note

  • OMAProblem(y, t, ...) = OMAProblem(y, y, t, ...)

FSDD

Frequency-Spatial Domain Decomposition method for OMA modal extraction

Fields

  • nothing

Reference

[1] L. Zhang, T. Wang and Y. Tamura. "A frequency–spatial domain decomposition (FSDD) method for operational modal analysis". Mechanical Systems and Signal Processing, 24: 1227-1239, 2010.

CovSSI

Covariance-driven Stochastic Subspace Identification method for OMA modal extraction

Fields

  • nothing

Reference

[1] C. Rainieri and G. Fabbrocino. "Operational Modal Analysis of Civil Engineering Structures: An Introduction and Guide for Applications". Springer, 2014.

[2] P. Peeters and G. De Roeck. "Reference-based stochastic subspace identification for output-only modal analysis". Mechanical Systems and Signal Processing, 13(6):855-878, 1999.

[3] L. Hermans and H. Van der Auweraer. "Modal testing and analysis of structures under operational conditions: Industrial applications". Mechanical Systems and Signal Processing, 13(2):193-216, 1999.

DataSSI

Data-driven Stochastic Subspace Identification method for OMA modal extraction

Fields

  • nothing

Reference

[1] C. Rainieri and G. Fabbrocino. "Operational Modal Analysis of Civil Engineering Structures: An Introduction and Guide for Applications". Springer, 2014.

[2] P. Peeters and G. De Roeck. "Reference-based stochastic subspace identification for output-only modal analysis". Mechanical Systems and Signal Processing, 13(6):855-878, 1999.

[3] L. Hermans and H. Van der Auweraer. "Modal testing and analysis of structures under operational conditions: Industrial applications". Mechanical Systems and Signal Processing, 13(2):193-216, 1999.

Related functions

xcorr

xcorr(y, yref, N)
xcorr(y, N)

Compute the cross-correlation matrix R of a signal y up to lag N-1.

Inputs

  • y::VecOrMat{Real}: Input signal

  • yref::VecOrMat{Real}: Reference input signal

  • N: Number of lags

Output

  • R: Cross-correlation matrix of size (ny, nyref, N), where ny is the number of signals (ny = 1 if y is a vector)

Note

  • xcorr(y, N) = xcorr(y, y, N)

half_csd

half_psd(Gyy, freq)
half_psd(y, yref, freq, fs, bs)
half_psd(y, freq, fs, bs)

Compute the half power spectral density matrix Gyy_half given the
spectral density matrix Gyy.

Inputs

  • Gyy: Power spectral density matrix of size (ni, no, nf), where ni is the number of inputs, no is the number of outputs, and nf is the number of frequency points.

  • freq::AbstractVector: Frequency vector (Hz)

Alternative inputs

  • y::VecOrMat{Real}: Output signal

  • yref::VecOrMat{Real}: Reference output signal

  • freq::AbstractVector: Frequency vector (Hz)

  • fs: Sampling frequency (Hz)

  • bs: Block size used to compute the cross-correlation matrix

Output

  • Gyy_half::Array{Complex, 3}: Half-spectrum matrix

References

[1] S. Chauhan. "Parameter estimation algorithms in operational modal analysis". In Proceedings of the 6th International Operational Modal Analysis Conference (IOMAC), Gijon, Spain, 2015.

[2] B. Peeters and H. Van der Auweraer, "PolyMax: A revolution in operational modal analysis". In Proceedings of the 1st International Operational Modal Analysis Conference (IOMAC), Copenhagen, Denmark, 2005.

[3] R. Brincker and C. E. Ventura. "Introduction to operational modal analysis". Wiley, 2015.

csd_from_tf

csd_from_tf(H, Gxx; id_input, id_output)

Compute the output cross-spectral density matrix Gyy given the frequency response function matrix H and the input cross-spectral density matrix Gxx.

Inputs

  • H::Array{Complex, 3}: Frequency response function matrix of size (no, ni, nf),

  • Gxx: Input power spectral density matrix of size (ni, ni, nf)

  • id_input::AbstractVector: Indices of input channels to consider (default: all inputs)

  • id_output::AbstractVector: Indices of output channels to consider (default: id_input)

Output

  • Gyy::Array{Complex, 3}: Output cross-spectral density matrix of size (ni, no, nf)

References

[1] B. Peeters and H. Van der Auweraer, "PolyMax: A revolution in operational modal analysis". In Proceedings of the 1st International Operational Modal Analysis Conference (IOMAC), Copenhagen, Denmark, 2005.

convert_csd

convert_csd(Gyy, freq; type = :dis)

Convert a cross-spectral density matrix to displacement type

Inputs

  • Gyy: Cross-spectral density matrix of size (ni, no, nf), where ni is the number of inputs, no is the number of outputs, and nf is the number of frequency points.

  • freq::AbstractVector: Frequency vector (Hz)

  • type::Symbol: Type of Gyy

    • :dis: displacement (default) - no conversion

    • :vel: velocity

    • :acc: acceleration

Output

  • Gyy_conv: Converted cross-spectral density matrix

poles_extraction

poles_extraction(prob, order, method; stabdiag)

Extract poles using Stochastic System Identification methods.

Inputs

  • prob::OMAProblem: Structure containing half-spectrum data and frequency vector

  • order::Int: Order of the system to identify

  • method::MdofOMA: OMA method to use for pole extraction

    • CovSSI: Covariance-based SSI (default)

    • DataSSI: Data-based SSI

  • stabdiag::Bool: Whether to compute stabilization diagram (default: false)

Output

  • poles::Vector{Complex}: Extracted poles

References [1] C. Rainieri and G. Fabbrocino. "Operational Modal Analysis of Civil Engineering Structures: An Introduction and Guide for Applications". Springer, 2014.

[2] P. Peeters and G. De Roeck. "Reference-based stochastic subspace identification for output-only modal analysis". Mechanical Systems and Signal Processing, 13(6):855-878, 1999.

[3] L. Hermans and H. Van der Auweraer. "Modal testing and analysis of structures under operational conditions: Industrial applications". Mechanical Systems and Signal Processing, 13(2):193-216, 1999.

modes_extraction

modes_extraction(prob, alg::FSDD; width, min_prom, max_prom , pks_indices, min_mac, avg_alg)

Extract modes from an Operational Modal Analysis problem using the Full Spectrum Decomposition and Decay (FSDD) method

Inputs

  • prob::OMAProblem: OMA problem containing the CSD matrix and frequency lines.

  • alg::SdofOMA: OMA method to use for pole extraction

    • FSDD: Frequency-spatial Domain Decomposition)

  • width::Int = 1: Width parameter for peak detection in the singular value spectrum

  • min_prom::Float64 = 0.: Minimum peak prominence for peak detection

  • max_prom::Float64 = Inf: Maximum peak prominence for peak detection

  • pks_indices::Vector{Int} = Int[]: Predefined indices of peaks

  • min_mac::Float64 = 0.9: Minimum MAC value to consider neighboring frequencies for mode shape estimation

  • avg_alg::Symbol = :sv: Averaging algorithm to use for mode shape estimation

    • :sv: Weighted averaging using singular values

    • :lin: Linear averaging

Outputs

  • poles::Vector{ComplexF64}: Extracted poles

  • ms::Array{ComplexF64, 2}: Extracted mode shapes

modes_extraction(prob, order, method; stabdiag)

Extract modes using Stochastic System Identification methods.

Inputs

  • prob::OMAProblem: Structure containing half-spectrum data and frequency vector

  • order::Int: Order of the system to identify

  • method: OMA method to use for pole extraction

    • CovSSI: Covariance-based SSI (default)

    • DataSSI: Data-based SSI

  • stabdiag::Bool: Whether to compute stabilization diagram (default: false)

Output

  • poles::Vector{Complex}: Extracted poles

  • ms::Matrix{Complex}: Extracted modal parameters

References [1] C. Rainieri and G. Fabbrocino. "Operational Modal Analysis of Civil Engineering Structures: An Introduction and Guide for Applications". Springer, 2014.

[2] P. Peeters and G. De Roeck. "Reference-based stochastic subspace identification for output-only modal analysis". Mechanical Systems and Signal Processing, 13(6):855-878, 1999.

[3] L. Hermans and H. Van der Auweraer. "Modal testing and analysis of structures under operational conditions: Industrial applications". Mechanical Systems and Signal Processing, 13(2):193-216, 1999.

half_csd_reconstruction

half_csd_reconstruction(res, poles, freq; lr, ur, type)

Reconstruct a half-spectrum from its residues and poles.

Inputs

  • res::Array{Complex, 3}: Residues corresponding to each pole

  • poles::Vector{Complex}: Poles extracted from the half-spectrum

  • freq::Vector{Float64}: Frequency vector

  • lr::Matrix{Complex, 2}: Lower residuals (default: zeros)

  • ur::Matrix{Complex, 2}: Upper residuals (default: zeros)

Output

  • H_rec::Array{Complex, 3}: Reconstructed FRF

compute_residuals

compute_residuals(prob, res, poles)

Compute lower and upper residuals of a frequency response function (FRF) given its residues and poles.

Inputs

  • prob::EMAProblem: Structure containing FRF data and frequency vector

  • res::Array{Complex, 3}: Residues corresponding to each pole

  • poles::Vector{Complex}: Poles extracted from the FRF

Outputs

  • lr::Matrix{Complex, 2}: Lower residuals

  • ur::Matrix{Complex, 2}: Upper residuals

compute_residuals(prob, res, poles)

Compute lower and upper residuals of a half-spectrum given its residues and poles.

Inputs

  • prob::OMAProblem: Structure containing half-spectrum data and frequency vector

  • res::Array{Complex, 3}: Residues corresponding to each pole

  • poles::Vector{T}: Poles extracted from the half-spectrum

Outputs

  • lr::Matrix{Complex, 2}: Lower residuals

  • ur::Matrix{Complex, 2}: Upper residuals

4 Example

4.1 Data preparation

# Structure parameters of the beam
L = 1.        # Length
b = 0.03      # Width
h = 0.01      # Thickness
S = b*h       # Cross-section area
Iz = b*h^3/12 # Moment of inertia

# Material parameters
E = 2.1e11  # Young's modulus
ρ = 7850.   # Density
ξ = 0.01    # Damping ratio

# Mesh
x = 0:0.05:L

# Acquisition parameters
sample_rate = 4096
block_size = 4096
fft_params = FFTParameters(sample_rate, block_size)

freq = fft_params.freq
t = fft_params.t
dt = fft_params.dt

# Mode calculation - Simply supported boundary conditions
beam = Beam(L, S, Iz, E, ρ)
ωn, kn = modefreq(beam, 2freq[end])
ms_ref = modeshape(beam, kn, x)

# Chirp excitation
F0 = 10.
chirp = SweptSine(F0, t[1], 0.8t[end], freq[1], freq[end], zero_end = true)
force = zeros(length(x), length(t))
force[2, :] .= excitation(chirp, t)

# Response calculation
prob = ForcedModalTimeProblem(ωn, ms_ref, ξ*ones(length(kn)), ms_ref'force, (zeros(length(x)), zeros(length(x))), t)
y = solve(prob).u

# OMA problem definition
tukeywin(x) = tukey(x, 0.5)
prob_oma = OMAProblem(y, t, sample_rate, block_size, win = tukeywin, frange = [1., 1000.])

4.2 Frequency-spatial Domain Decomposition

pks_indices = [23, 94, 211, 375, 586, 845]
p_fsdd, ms_fsdd = modes_extraction(prob_oma, FSDD(), pks_indices = pks_indices)
ms_fsdd_real = real_normalization(ms_fsdd)

# Mode shape scaling for comparison using the modal scale factor
scaling = msf(ms_fsdd_real, ms_ref[:, 1:size(ms_fsdd_real, 2)])
ms_fsdd_scaled = ms_fsdd_real .* scaling'

4.3 EMA-based OMA

4.3.1 Poles extraction using LSCF

# EMA-MDOF pole stability analysis
sol_stab = stabilization(prob_oma, 30, LSCF())

# Extraction of stable poles
stab_poles = sol_stab.poles[end][.!isnan.(sol_stab.poles[end])]

# Plot stabilization diagram
stabilization_plot(sol_stab)

4.3.2 Mode shapes extraction

# Computation of the mode residues
res = mode_residues(prob_oma, stab_poles)

# Extraction of the mode shapes
ms_ema = modeshape_extraction(res, stab_poles, LSCF(), modetype = :oma)[1]

# Convert to real mode shapes
ms_ema_real = real_normalization(ms_ema)

# Mode shape scaling for comparison using the modal scale factor
scaling = msf(ms_ema_real, ms_ref[:, 1:size(ms_ema_real, 2)])
ms_ema_scaled = ms_ema_real .* scaling'

4.3.3 Half-spectrum reconstruction

# Half-spectrum reconstruction - without residuals
Syy = half_csd_reconstruction(res, stab_poles, prob_oma.freq)

# Half-spectrum reconstruction - with residuals
lr, ur = compute_residuals(prob_oma, res, stab_poles)
Syy2 = half_csd_reconstruction(res, stab_poles, prob_oma.freq, lr = lr, ur = ur)

4.4 Covariance-driven Stochastic Subspace Identification

# Stabilization analysis
sol_stab_cssi = stabilization(prob_oma, 30, CovSSI())

# Plot stabilization diagram
stabilization_plot(sol_stab_cssi)
# Extraction of stable modes
p_cssi, ms_cssi = modes_extraction(prob_oma, 30, CovSSI())
ms_cssi_real = real_normalization(ms_cssi)

# Mode shape scaling for comparison using the modal scale factor
scaling = msf(ms_cssi_real, ms_ref[:, 1:size(ms_cssi_real, 2)])
ms_cssi_scaled = ms_cssi_real .* scaling'

4.5 Data-driven Stochastic Subspace Identification

# Stabilization analysis
sol_stab_dssi = stabilization(prob_oma, 30, DataSSI())

# Plot stabilization diagram
stabilization_plot(sol_stab_dssi)
# Extraction of stable modes
p_dssi, ms_dssi = modes_extraction(prob_oma, 30, DataSSI())
ms_dssi_real = real_normalization(ms_dssi)

# Mode shape scaling for comparison using the modal scale factor
scaling = msf(ms_dssi_real, ms_ref[:, 1:size(ms_dssi_real, 2)])
ms_dssi_scaled = ms_dssi_real .* scaling'

Footnotes

  1. W. Tong, Z. Lingmi and T. Yukio. “An operational modal analysis method in frequency and spatial domain”. Earthquake engineering and engineering vibration, 4(2): 295-300, 2005.↩︎

  2. R. Brincker, C. E. Ventura and P. Andersen. “Damping estimation by Frequency Domain Decomposition”. In Proceedings of IMAC XIX, Orlando, Florida, USA. 2001.↩︎

  3. R. Brincker and L. Zhang. “Frequency domain decomposition revisited”. In PProceedings of the 3rd International Operational Modal Analysis conference (IOMAC), Portonovo, Italy, 2009.↩︎

  4. B. Peeters and H. Van der Auweraer. “Polymax: A revolution in operational modal analysis”. In Proceedings of the 1st International Operational Modal Analysis conference (IOMAC), Copenhagen, Denmark, 2005.↩︎

  5. R. Brincker and C. Ventura. “Introduction to Operational Modal Analysis”. Wiley, 2015.↩︎