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:

• Mdof EMA-based methods ;

• Stochastic Subspace Identification methods.

1 EMA-based methods for OMA

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 as¹²: \[ 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 expression³: \[ 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.

1.2 Modal parameters extraction

Once the half-spectrum is computed, the existing Mdof EMA methods in StructuralVibration.jl can be directly applied to extract the poles from the half-spectrum data. More details on the usage of these methods can be found in the Mdof EMA documentation. Regarding the extraction of mode shapes in OMA using EMA-based methods, a procedure similar to the one used in EMA can’t be followed since the input forces are unknown. In this case, it is generally assumed that the rank of the residues is equal to 1, allowing to extract the mode shapes from the residues using a Singular Value Decomposition (SVD). Specifically, for a given mode \(i\), the residue matrix \({}_i\mathbf{R}\) can be decomposed as: \[ {}_i\mathbf{R} = \mathbf{U}_i\, \mathbf{S}_i\, \mathbf{V}_i^H, \]

where \(\mathbf{U}_i\) and \(\mathbf{V}_i\) are the matrices containing the left and right singular vectors, respectively, and \(\mathbf{S}_i\) is the diagonal matrix of singular values. The mode shape \(\boldsymbol{\phi}_i\) is then approximated by the first left singular vector corresponding to the highest singular value: \[ \boldsymbol{\phi}_i \approx \mathbf{U}_{i,1} \]

where \(\mathbf{U}_{i,1}\) is the first column of \(\mathbf{U}_i\). However, this assumption may not always hold true, especially in cases where the input excitation is not perfectly white noise or when multiple modes are closely spaced. In this case, more than one singular singular value may be significant, and the mode shape extraction may require additional considerations. In StructuralVibration.jl, the following procedure is adopted to extract the mode shapes from the residues:

1. For each mode \(i\), compute the SVD of the residue matrix \({}_i\mathbf{R}\) to obtain \(\mathbf{U}_i\), \(\mathbf{S}_i\).

2. Filter the singular values in \(\mathbf{S}_i\) based on a predefined threshold to determine the number \(j\) of significant singular values. Here, the threshold is set to 1% of the largest singular value.

3. For the first mode, select \(\boldsymbol{\phi}_i = \mathbf{U}_{i,j}\) as the mode shape, based on energy considerations. The selected singular vector corresponds to the highest energy spread. Generally, this will be the first singular vector.

4. For subsequent modes, select the singular vector \(\mathbf{U}_{i,j}\) that minimizes the Modal Assurance Criterion (MAC) with respect to the previously extracted mode shapes, ensuring orthogonality and consistency among the extracted modes.

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 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.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 \]

2. 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} \]

3. 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.

4. 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 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.

2. 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 \]

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

4. 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, ...)

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::OMAModalExtraction: 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, 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::OMAModalExtraction: 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, ismodal = true)
y = solve(prob).u

4.2 EMA-based OMA

4.2.1 Poles extraction using LSCF

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

# 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.2.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.2.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.3 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.4 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. R. Brincker and L. Zhang. “Frequency domain decomposition revisited”. In PProceedings of the 3rd International Operational Modal Analysis conference (IOMAC), Portonovo, Italy, 2009.

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. Ventura. “Introduction to Operational Modal Analysis”. Wiley, 2015.