pyOMA – Single-Setup Operational Modal Analysis#

This notebook walks through a complete Operational Modal Analysis (OMA) workflow using a single measurement setup. It uses the example data bundled with the repository so you can run it immediately after installation.

Workflow overview

Step

Class / function

Purpose

1

GeometryProcessor

Load structural geometry (nodes, lines)

2

PreProcessSignals

Load measurement data and assign metadata

3

PreProcessSignals

Decimate, filter, compute correlations

4

BRSSICovRef (or other)

Identify modal parameters

5

StabilCalc / StabilCluster

Compute stabilisation diagram

6

Results & export

Inspect frequencies, damping, mode shapes

Tip: Run cells top-to-bottom with Shift+Enter. Each section begins with a short explanation of what it does and why.

0 Imports and display setup#

from pathlib import Path
import numpy as np
import matplotlib

# Use the interactive widget backend when available, fall back to inline.
try:
    import ipympl  # noqa: F401 – presence check only
    %matplotlib widget
except ImportError:
    %matplotlib inline

import matplotlib.pyplot as plt

# --- pyOMA imports -----------------------------------------------------------
from pyOMA.core import (
    GeometryProcessor,
    PreProcessSignals,
    SignalPlot,
    BRSSICovRef,
    SSIData,
    PLSCF,
    VarSSIRef,
    StabilCalc,
    StabilCluster,
    StabilPlot,
    ModeShapePlot,
)

print('pyOMA imported successfully')

1 Paths#

The cell below locates the example data that ships with the repository (tests/files/). When you apply this workflow to your own project, replace EXAMPLE_DATA with the folder that contains your measurement files and set SETUP_DIR to the folder for the specific setup you want to analyse.

import pyOMA

# Root of the repository — works whether pyOMA is installed or run in-place
REPO_ROOT    = Path(pyOMA.__file__).parent.parent
EXAMPLE_DATA = REPO_ROOT / 'tests' / 'files'
SETUP_DIR    = EXAMPLE_DATA / 'measurement_1'

assert EXAMPLE_DATA.exists(), f'Example data not found at {EXAMPLE_DATA}'
print(f'Example data : {EXAMPLE_DATA}')
print(f'Setup folder : {SETUP_DIR}')

2 Structural geometry#

GeometryProcessor stores the 3-D coordinates of measurement nodes, the connectivity lines between them, and optional parent–child relationships used when sensors do not sit directly on structural nodes.

File

Content

Required

grid.txt

Node names and (x, y, z) coordinates

Yes

lines.txt

Pairs of connected node names

No

parent_child_assignments.txt

Sensor-to-node offsets

No

See Input File Formats in the documentation for the exact column layout.

geometry_data = GeometryProcessor.load_geometry(
    nodes_file=EXAMPLE_DATA / 'grid.txt',
    lines_file=EXAMPLE_DATA / 'lines.txt',
    parent_childs_file=EXAMPLE_DATA / 'parent_child_assignments.txt',
)

print(f'{len(geometry_data.nodes)} nodes, {len(geometry_data.lines)} lines loaded')

3 Loading measurement signals#

3.2 Inspect raw signals and spectra#

plot_signals shows time histories and auto-spectral densities for each channel. Use this to check for offsets, spikes, and to identify the frequency band that contains structural modes.

plot_signals = SignalPlot(prep_signals)
plot_signals.plot_signals(per_channel_axes=True, n_lines=512)
plt.suptitle('Raw signals', fontsize=12)
plt.tight_layout()

3.3 Geometry with sensor positions#

ModeShapePlot can visualise the structural geometry together with the channel-DOF arrows showing each sensor’s measurement direction before any identification is run.

# Nodes (grey dots), structural lines and channel-DOF direction arrows
geo_plot = ModeShapePlot(geometry_data=geometry_data, prep_signals=prep_signals)
geo_plot.reset_view()
geo_plot.draw_nodes()
geo_plot.draw_lines()
geo_plot.draw_chan_dofs()
geo_plot.refresh_parent_childs(False)
geo_plot.subplot.view_init(elev=20, azim=110)
display(geo_plot.fig)

3.4 Correlation functions and power spectral densities#

After calling prep_signals.correlation() you can inspect the estimated correlation functions and their PSDs — useful to verify m_lags is large enough and the frequency content looks reasonable.

# Run after pre-processing (Step 4) to inspect correlation functions
sig_plot2 = SignalPlot(prep_signals)
sig_plot2.plot_signals(per_channel_axes=True, timescale='lags', psd_scale='db')
plt.suptitle('Correlation functions and PSDs', fontsize=12)
plt.show()

4 Signal pre-processing#

4.1 Decimation#

Reduces the sampling rate by an integer factor. An anti-aliasing filter is applied automatically before down-sampling.

The example data was recorded at 256 Hz. Decimating twice by 3 gives ≈ 28.4 Hz — sufficient for modes below ~12 Hz.

prep_signals.decimate_signals(3)
prep_signals.decimate_signals(3)

print(f'After decimation: {prep_signals.sampling_rate:.2f} Hz  '
      f'({prep_signals.total_time_steps} samples  '
      f'{prep_signals.duration:.1f} s)')

4.2 Correlation functions#

Covariance-driven SSI and PLSCF need the cross-correlation matrix between all channels and the reference channels. pyOMA provides two estimators:

Method

Function

Notes

Welch

corr_welch

Fast; slight damping bias from windowing

Blackman–Tukey

corr_blackman_tukey

Better resolution; recommended for SSI-cov

m_lags must exceed num_block_columns + num_block_rows (set in Step 5).

prep_signals.corr_blackman_tukey(m_lags=200)

print(f'Correlation matrix shape: {prep_signals.corr_matrix.shape}'
      f'  (channels x refs x lags)')

5 System identification#

pyOMA implements several classical OMA algorithms. Pick the class and call the two core methods — no config file required.

Class

Algorithm

Call sequence

BRSSICovRef

Covariance-driven SSI

build_toeplitz_cov → compute_modal_params

SSIData

Data-driven SSI

build_block_hankel → compute_modal_params

VarSSIRef

SSI + uncertainty

build_subspace_mat → compute_state_matrices → prepare_sensitivities → compute_modal_params

PLSCF

Poly-ref. LS Complex Freq.

build_half_spectra → compute_modal_params

Key parameters

  • num_block_columns — number of block columns in the Toeplitz / Hankel matrix. Must satisfy num_block_columns + num_block_rows < m_lags (here < 200). Typical range: 50–150.

  • max_model_order — modal parameters are estimated at every order 1 … max. Set to roughly 2× the number of expected modes.

# ── Choose a method by uncommenting one block ─────────────────────────────────

# --- BRSSICovRef (Covariance-driven SSI — most common) -----------------------
modal_data = BRSSICovRef(prep_signals)
modal_data.build_toeplitz_cov(num_block_columns=100)
modal_data.compute_modal_params(max_model_order=40)

# --- SSIData (Data-driven SSI) -----------------------------------------------
# modal_data = SSIData(prep_signals)
# modal_data.build_block_hankel(num_block_rows=100)
# modal_data.compute_modal_params(max_model_order=40)

# --- PLSCF (frequency-domain) -------------------------------------------------
# modal_data = PLSCF(prep_signals)
# modal_data.build_half_spectra(
#     nperseg=200,           # same as m_lags above
#     begin_frequency=0.0,   # Hz
#     end_frequency=12.0,    # Hz
# )
# modal_data.compute_modal_params(max_model_order=40)
# ─────────────────────────────────────────────────────────────────────────────

print(f'Method   : {type(modal_data).__name__}')
print(f'Orders   : 1 - {modal_data.max_model_order}')
print(f'Channels : {modal_data.num_analised_channels}')

stabil_calc = StabilCluster(modal_data)

stabil_calc.calculate_stabilization_masks(
    order_range=(0, 1, modal_data.max_model_order),
    d_range=(0, 0.10),   # damping 0–10 %
    df_max=0.01,          # max relative frequency change between orders
    dd_max=0.05,          # max relative damping change between orders
    dmac_max=0.05,        # max MAC difference between orders
)

n_stable = stabil_calc.masks['mask_stable'].sum()
print(f'{n_stable} stable poles found across all model orders')

6.1 Static stabilisation diagram#

Marker colours indicate stability status:

  • Light grey dots — all poles

  • Coloured markers — poles stable in individual criteria

  • Black stars / filled markers — stable in all active criteria simultaneously

stabil_plot = StabilPlot(stabil_calc)
stabil_plot.plot_stabil('plot_stable')
plt.tight_layout
display(stabil_plot.fig)

6.2 Interactive stabilisation GUI (requires pyOMA[jupyter])#

Install the optional Jupyter dependencies with:

pip install "pyOMA[jupyter]"

then run the cell below. The widget lets you click on poles to select them and see the corresponding mode shape animated in 3-D.

try:
    from pyOMA.GUI.JupyterGUI import StabilGUIWeb
    stabil_widget, snap_cursor = StabilGUIWeb(stabil_plot)
    display(stabil_widget)
except ImportError:
    print('ipywidgets / ipympl not installed.  '
          'Run: pip install "pyOMA[jupyter]"')

7 Automated mode selection and result export#

automatic_clearing removes poles that are clearly spurious by applying soft statistical criteria. For a fully headless workflow you can follow it with automatic_classification and automatic_selection (see the API reference).

export_results writes a tab-separated summary table with frequencies, damping ratios, and quality metrics for all selected modes.

# Automated clearing
try:
    stabil_calc.automatic_clearing()
    print(f'Poles remaining after clearing: '
          f'{stabil_calc.masks["mask_stable"].sum()}')
except Exception as exc:
    print(f'Automatic clearing skipped ({exc})')

# Export to a text file
import tempfile, os
with tempfile.NamedTemporaryFile(suffix='.txt', delete=False, mode='w') as tf:
    export_path = tf.name

stabil_calc.export_results(export_path)
print(f'\nResults written to {export_path}')

# Print the first lines
with open(export_path) as f:
    for line in f.readlines()[:20]:
        print(line, end='')
os.unlink(export_path)

8 Mode shape visualisation#

ModeShapePlot animates selected mode shapes on the 3-D geometry.

msh_plot = ModeShapePlot(
    geometry_data=geometry_data,
    stabil_calc=stabil_calc,
    modal_data=modal_data,
    prep_signals=prep_signals,
    amplitude=20,
)

# Interactive Jupyter version (requires pyOMA[jupyter])
try:
    from pyOMA.GUI.JupyterGUI import PlotMSHWeb
    display(PlotMSHWeb(msh_plot))
except ImportError:
    print('Jupyter widget not available.  '
          'For the desktop GUI: from pyOMA.GUI.PlotMSHGUI import start_msh_gui; '
          'start_msh_gui(msh_plot)')

9 Saving and resuming work#

Every major object serialises to a compressed NumPy archive. Checkpoint between stages to avoid recomputing expensive steps.

import tempfile

with tempfile.TemporaryDirectory() as tmpdir:
    tmp = Path(tmpdir)

    # --- save ----------------------------------------------------------------
    prep_signals.save_state(tmp / 'prep_signals.npz')
    modal_data.save_state(tmp / 'modal_data.npz')
    stabil_calc.save_state(tmp / 'stabil_data.npz')

    # --- reload --------------------------------------------------------------
    prep2   = PreProcessSignals.load_state(tmp / 'prep_signals.npz')
    modal2  = BRSSICovRef.load_state(tmp / 'modal_data.npz', prep2)
    stabil2 = StabilCalc.load_state(tmp / 'stabil_data.npz', modal2)

    import numpy.testing as npt
    npt.assert_array_equal(prep2.signals, prep_signals.signals)
    npt.assert_array_equal(modal2.modal_frequencies,
                           modal_data.modal_frequencies)
    print('Save / load round-trip: OK')

Next steps#

  • Multi-setup analysis — see PogerSSICovRef and MergePoSER for the PoGer / PoSER methodology when all channels cannot be measured simultaneously.

  • Uncertainty quantification — VarSSIRef adds frequency and damping uncertainty estimates to the identified modes.

  • Desktop GUI — from pyOMA.GUI.StabilGUI import start_stabil_gui after pip install "pyOMA[gui]".

  • API reference — https://py-oma.readthedocs.io