Coverage for src / sdynpy / doc / sdynpy_vibration_test.py: 9%
611 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-17 21:21 +0000
1# -*- coding: utf-8 -*-
2"""
3Data container and automated plotting tools for random vibration tests.
5Public API
6----------
7RandomVibTest
8 Container for a random vibration test data set. Holds time history data,
9 cross-power spectral densities, and specification PSDs. Provides methods
10 to compute CPSDs/FRFs and produce a standard suite of diagnostic figures.
11optimal_subset
12 Identify the time window with the most stable RMS amplitude.
13make_multi_figures
14 Create paginated grids of subplots for multi-channel data.
15dynamic_barh
16 Horizontal bar chart with automatic multi-column layout.
17response_outside_mask
18 Boolean mask of frequency lines where a response exceeds a limit.
20Copyright 2022 National Technology & Engineering Solutions of Sandia,
21LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S.
22Government retains certain rights in this software.
24This program is free software: you can redistribute it and/or modify
25it under the terms of the GNU General Public License as published by
26the Free Software Foundation, either version 3 of the License, or
27(at your option) any later version.
29This program is distributed in the hope that it will be useful,
30but WITHOUT ANY WARRANTY; without even the implied warranty of
31MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
32GNU General Public License for more details.
34You should have received a copy of the GNU General Public License
35along with this program. If not, see <https://www.gnu.org/licenses/>.
36"""
38from ..core.sdynpy_data import (TransferFunctionArray, CoherenceArray, MultipleCoherenceArray, PowerSpectralDensityArray, TimeHistoryArray, NDDataArray)
39from ..core.sdynpy_geometry import (Geometry,GeometryPlotter,ShapePlotter)
40from ..core.sdynpy_coordinate import CoordinateArray, coordinate_array as sd_coordinate_array
41from ..modal.sdynpy_signal_processing_gui import SignalProcessingGUI
42from ..fileio.sdynpy_rattlesnake import RattlesnakeData, RattlesnakeRandomEnvironmentData
43import matplotlib.pyplot as plt
44from matplotlib.patches import Patch
45from matplotlib.ticker import LogLocator
46import os
47import netCDF4 as nc4
48import numpy as np
49import pandas as pd
50import sdynpy as sdpy
51import warnings
52import math
53from scipy.stats import kurtosis
55# Define Functions
56def optimal_subset(time_data: TimeHistoryArray,num_subset_samples: int, amplitude_factor: float=0.5):
57 """Find the start index of the time subset with most stable amplitude.
59 Searches for a contiguous window of ``num_subset_samples`` samples that
60 minimises the standard deviation of RMS amplitude across channels while
61 remaining above a minimum amplitude threshold.
63 Parameters
64 ----------
65 time_data : TimeHistoryArray
66 Time-domain data whose ordinate has shape ``(n_channels, n_samples)``.
67 num_subset_samples : int
68 Length of the sliding window in samples.
69 amplitude_factor : float, optional
70 Fraction of the peak-to-trough mean range used as the minimum
71 acceptable window mean amplitude. Windows below this threshold are
72 penalised. Default is 0.5.
74 Returns
75 -------
76 int
77 Sample index where the optimal subset begins.
78 """
79 ordinate_rms = np.sqrt(np.mean(time_data.ordinate**2,axis=0))
81 # Compute Cumulative Sums
82 cumsum1 = np.concatenate([[0],np.cumsum(ordinate_rms)])
83 cumsum2 = np.concatenate([[0],np.cumsum(ordinate_rms**2)])
85 # Compute Variance and Standard Deviation
86 window_mean = (cumsum1[num_subset_samples:] - cumsum1[:-num_subset_samples]) / num_subset_samples
87 window_variance = (cumsum2[num_subset_samples:] - cumsum2[:-num_subset_samples]) / num_subset_samples - window_mean**2
88 window_std = np.sqrt(window_variance)
90 # Add threashold for mean data
91 threshold = (max(window_mean) - min(window_mean))*amplitude_factor + min(window_mean)
92 score = np.where(window_mean > threshold, window_std, np.inf)
94 # Define Window Start
95 start_index = np.argmin(score)
97 return start_index
99def make_multi_figures(total_plots, max_cols=3, max_subplots_per_fig=18,xlabel='Frequency (Hz)',ylabel='|PSD| (EU$^2$/Hz)',fig_size_x=8.5,fig_size_y=11):
100 """Create one or more figures, each containing a grid of subplots.
102 Parameters
103 ----------
104 total_plots : int
105 Total number of subplots to create across all figures.
106 max_cols : int, optional
107 Maximum number of subplot columns per figure. Default is 3.
108 max_subplots_per_fig : int, optional
109 Maximum subplots per figure before a new figure is started. Default
110 is 18.
111 xlabel : str, optional
112 Label for the x-axis of each subplot. Default is ``'Frequency (Hz)'``.
113 ylabel : str, optional
114 Label for the y-axis of each subplot. Default is ``'|PSD| (EU$^2$/Hz)'``.
115 fig_size_x : float, optional
116 Figure width in inches. Default is 8.5.
117 fig_size_y : float, optional
118 Figure height in inches. Default is 11.
120 Returns
121 -------
122 figs : list of matplotlib.figure.Figure
123 List of created figures.
124 all_axes : list of list of matplotlib.axes.Axes
125 Nested list where ``all_axes[i]`` contains the axes for ``figs[i]``.
126 """
127 figs = []
128 all_axes = []
130 start_idx = 0
132 while start_idx < total_plots:
133 # Determine how many plots in the figure
134 end_idx = min(start_idx + max_subplots_per_fig,total_plots)
135 n_plots = end_idx - start_idx
137 cols = min(max_cols,n_plots)
138 rows = math.ceil(n_plots/cols)
139 max_rows = math.ceil(max_subplots_per_fig/max_cols)
141 fig = plt.figure(figsize=(fig_size_x,fig_size_y))
142 gs = fig.add_gridspec(max_rows, max_cols)
144 axes = []
145 global_head = None
147 # Create subplots into the GridSpec
148 for i in range(n_plots):
149 r = i // cols
150 c = i % cols
152 # First Subplot becomes the global head
153 if global_head is None:
154 ax = fig.add_subplot(gs[r,c])
155 global_head = ax
156 else:
157 ax = fig.add_subplot(gs[r,c], sharex=global_head, sharey=global_head)
159 axes.append(ax)
161 # Set global axis labels on bottom left subplot
162 if c == 0:
163 ax.set_ylabel(ylabel)
164 ax.tick_params(labelleft=True)
165 else:
166 ax.tick_params(labelleft=False)
167 # ax.set_yticklabels([])
169 if r == rows - 1:
170 ax.set_xlabel(xlabel)
171 ax.tick_params(labelbottom=True)
172 else:
173 ax.tick_params(labelbottom=False)
175 for i in range(n_plots, max_rows * max_cols):
176 fig.add_subplot(gs[i // max_cols, i % max_cols]).set_visible(False)
178 figs.append(fig)
179 all_axes.append(axes)
181 start_idx += n_plots
183 return figs, all_axes
185def dynamic_barh(values, labels=None, max_per_col=35):
186 """Create a horizontal bar chart with automatic multi-column layout.
188 Parameters
189 ----------
190 values : array_like
191 Numeric values to plot as horizontal bars.
192 labels : CoordinateArray or None, optional
193 Labels for the y-axis ticks. If a ``CoordinateArray`` is provided the
194 tick labels are set accordingly. Default is ``None``.
195 max_per_col : int, optional
196 Maximum number of bars per column before a new column is added.
197 Default is 35.
199 Returns
200 -------
201 fig : matplotlib.figure.Figure
202 The created figure.
203 axes : list of matplotlib.axes.Axes
204 List of axes, one per column.
205 """
206 n = len(values)
208 # ---- Determine number of columns ----
209 cols = math.ceil(n / max_per_col)
210 rows = 1
212 # ---- Split values into chunks for each column ----
213 chunks = [
214 values[i * max_per_col : (i + 1) * max_per_col]
215 for i in range(cols)
216 ]
218 if isinstance(labels,CoordinateArray):
219 label_chunks = [
220 labels[i * max_per_col : (i + 1) * max_per_col]
221 for i in range(cols)
222 ]
223 else:
224 label_chunks = [None] * cols
226 # ---- Automatic figure height based on the tallest column ----
227 h_per_bar = 0.35
228 base_height = 1.0
229 tallest_col = max(len(chunk) for chunk in chunks)
230 fig_height = base_height + h_per_bar * tallest_col
232 # Wider figure for more columns
233 fig_width = 6 * cols
235 fig, axes = plt.subplots(
236 rows, cols,
237 figsize=(fig_width, fig_height),
238 sharex=True
239 )
241 if cols == 1:
242 axes = [axes]
243 else:
244 axes = axes.flatten()
246 # ---- Plot each column ----
247 for ax, chunk, lbl_chunk in zip(axes, chunks, label_chunks):
248 y = range(len(chunk))
249 ax.barh(y, chunk,color='deepskyblue',edgecolor='black', linewidth=1)
251 if isinstance(lbl_chunk,CoordinateArray):
252 ax.set_yticks(y)
253 ax.set_yticklabels(lbl_chunk)
254 else:
255 ax.set_yticks(y)
257 ax.invert_yaxis()
258 ax.margins(y=0)
259 ax.set_axisbelow(True)
260 ax.xaxis.grid(True)
262 # Apply shared y-limits
263 ax.set_ylim(axes[0].get_ylim())
265 # Remove unused axes (shouldn't happen but safe)
266 for ax in axes[len(chunks):]:
267 ax.remove()
269 fig.tight_layout()
270 return fig, axes
272def response_outside_mask(response, response_limit, mask_type: str = 'above'):
273 """Compute a boolean mask of where response is outside the specified limit.
275 Parameters
276 ----------
277 response : NDDataArray
278 The measured response data array.
279 response_limit : NDDataArray
280 The limit data array to compare against.
281 mask_type : str, optional
282 Direction of comparison: 'above' to flag where response exceeds the
283 limit, 'below' to flag where response falls below the limit.
284 Default is 'above'.
286 Returns
287 -------
288 np.ndarray
289 Boolean array that is True where the response is outside the limit.
291 Raises
292 ------
293 ValueError
294 If mask_type is not 'above' or 'below'.
295 """
296 response_in_limit_mask = np.isin(response.abscissa, response_limit.abscissa)
297 limit_in_response_mask = np.isin(response_limit.abscissa, response.abscissa)
298 response_exsists_mask = np.invert(np.isnan(response.ordinate))
299 limit_exists_mask = np.invert(np.isnan(response_limit.ordinate))
301 if mask_type == 'above':
302 outside_mask = np.logical_and(np.where(response_in_limit_mask, np.abs(np.abs(response.ordinate)), 0) > np.where(limit_in_response_mask, np.abs(np.abs(response_limit.ordinate)), 0), np.where(response_in_limit_mask, np.abs(response.ordinate), 0) > np.where(limit_in_response_mask, np.abs(response_limit.ordinate), 0))
303 elif mask_type == 'below':
304 outside_mask = np.logical_and(np.where(response_in_limit_mask, np.abs(np.abs(response.ordinate)), 0) < np.where(limit_in_response_mask, np.abs(np.abs(response_limit.ordinate)), 0), np.where(response_in_limit_mask, np.abs(response.ordinate), 0) < np.where(limit_in_response_mask, np.abs(response_limit.ordinate), 0))
305 else:
306 raise ValueError(f"mask_type must be 'above' or 'below', got '{mask_type}'")
308 outside_mask = np.logical_and.reduce((outside_mask, response_in_limit_mask, response_exsists_mask, limit_exists_mask))
310 return outside_mask
312def _shade_out_of_limit(ax, freqs, mask, limit_ordinate, ylim_bound, color, alpha, label):
313 """Shade from the abort limit curve to the plot edge for out-of-limit regions.
315 For each contiguous run of True values in *mask*, builds an x array that
316 extends half the local frequency spacing beyond each edge so that even a
317 single isolated spectral line produces a visible shaded band. The y extent
318 runs from the abort limit ordinate (extended at constant value to the padded
319 edges) to *ylim_bound*, preserving the original fill-between appearance.
320 Works correctly for non-uniform (e.g. octave-band) frequency spacing.
322 Parameters
323 ----------
324 ax : matplotlib.axes.Axes
325 The axes on which to draw the shading.
326 freqs : np.ndarray, shape (N,)
327 Frequency abscissa values corresponding to *mask*.
328 mask : np.ndarray of bool, shape (N,)
329 True at each frequency line that is outside the limit.
330 limit_ordinate : np.ndarray, shape (N,)
331 Abort limit ordinate values at each frequency in *freqs*.
332 ylim_bound : float
333 The plot-edge y value to shade toward (``min(ylim)`` for below,
334 ``max(ylim)`` for above).
335 color : str
336 Patch fill color.
337 alpha : float
338 Patch transparency.
339 label : str
340 Legend label applied to the first patch; subsequent patches use
341 ``'_nolegend_'`` to avoid duplicate legend entries.
342 """
343 indices = np.where(mask)[0]
344 if len(indices) == 0:
345 return
347 # Split into contiguous groups
348 breaks = np.where(np.diff(indices) > 1)[0] + 1
349 groups = np.split(indices, breaks)
351 labeled = False
352 for group in groups:
353 i0, i1 = int(group[0]), int(group[-1])
355 # Left edge: midpoint to the previous line; use right-side spacing at boundary
356 if i0 > 0:
357 x0 = freqs[i0] - (freqs[i0] - freqs[i0 - 1]) / 2
358 else:
359 x0 = freqs[i0] - (freqs[i0 + 1] - freqs[i0]) / 2 if len(freqs) > 1 else freqs[i0]
361 # Right edge: midpoint to the next line; use left-side spacing at boundary
362 if i1 < len(freqs) - 1:
363 x1 = freqs[i1] + (freqs[i1 + 1] - freqs[i1]) / 2
364 else:
365 x1 = freqs[i1] + (freqs[i1] - freqs[i1 - 1]) / 2 if len(freqs) > 1 else freqs[i1]
367 # Build x array with padded endpoints; extend limit ordinate at constant value
368 x_vals = np.concatenate([[x0], freqs[i0:i1 + 1], [x1]])
369 y_limit = np.abs(np.concatenate([[limit_ordinate[i0]], limit_ordinate[i0:i1 + 1], [limit_ordinate[i1]]]))
371 ax.fill_between(x_vals, y_limit, ylim_bound, color=color, alpha=alpha,
372 edgecolor='none', label=label if not labeled else '_nolegend_')
373 labeled = True
375# class TransientVibTest:
377# class RandomVibTest:
378class RandomVibTest:
379 """Container for a random vibration test data set with automated plotting.
381 Holds time history data, cross-power spectral densities, specifications,
382 and data for a random vibration test environment. Provides methods to
383 compute CPSDs/FRFs and produce a standard suite of diagnostic plots.
385 Parameters
386 ----------
387 source_data : RattlesnakeRandomEnvironmentData, optional
388 Parsed Rattlesnake environment data for the test environment.
389 geometry : str or Geometry, optional
390 Path to a geometry file or a ``Geometry`` object for 3-D visualization.
391 coordinate : CoordinateArray, optional
392 Full channel coordinate array (all channels).
393 control_coordinate : CoordinateArray, optional
394 Subset of coordinates for control/reference channels.
395 response_coordinate : CoordinateArray, optional
396 Subset of coordinates for response channels.
397 excitation_coordinate : CoordinateArray, optional
398 Subset of coordinates for excitation (drive) channels.
399 time_data : list of TimeHistoryArray, optional
400 List of time history arrays, one per data set.
401 cpsd : list of PowerSpectralDensityArray or None, optional
402 Pre-computed cross-power spectral densities. Pass ``None`` (default)
403 to have them computed on demand.
404 specification_cpsd : PowerSpectralDensityArray, optional
405 Target PSD specification.
406 specification_warning_psd : PowerSpectralDensityArray, optional
407 Warning-level tolerance bands around the specification.
408 specification_abort_psd : PowerSpectralDensityArray, optional
409 Abort-level tolerance bands around the specification.
410 averages : int, optional
411 Number of spectral averages used when computing CPSDs.
412 overlap : float, optional
413 Fractional overlap between frames (0 to 1).
414 fft_lines : int, optional
415 Number of FFT lines.
416 window : str, optional
417 Name of the windowing function (e.g. ``'hann'``).
418 samples_per_frame : int, optional
419 Number of samples per FFT frame.
420 start_time : float or None, optional
421 Start time (seconds) for the analysis window within each data set.
422 ``None`` triggers automatic selection via ``optimal_subset``.
423 units : np.ndarray, optional
424 String array of engineering-unit labels, one per channel.
425 frf : list of TransferFunctionArray or None, optional
426 Pre-computed frequency response functions.
427 oct_order : int, optional
428 Octave band order for octave-averaged results.
429 """
430 def __init__(self,
431 source_data: RattlesnakeRandomEnvironmentData = None,
432 geometry: str | Geometry = None,
433 coordinate: CoordinateArray = None,
434 control_coordinate: CoordinateArray = None,
435 response_coordinate: CoordinateArray = None,
436 excitation_coordinate: CoordinateArray = None,
437 time_data: list[TimeHistoryArray] = None,
438 cpsd: list[PowerSpectralDensityArray] = None,
439 specification_cpsd: PowerSpectralDensityArray = None,
440 specification_warning_psd: PowerSpectralDensityArray = None,
441 specification_abort_psd: PowerSpectralDensityArray = None,
442 averages: int = None,
443 overlap: float = None,
444 fft_lines: int = None,
445 window: str = None,
446 samples_per_frame: int = None,
447 start_time: float = None,
448 units: np.ndarray = None,
449 frf: list[TransferFunctionArray] = None,
450 oct_order: int = None,
451 ):
452 """Store all test data attributes.
454 Every named parameter is assigned directly to ``self``. See the class
455 docstring for parameter descriptions. ``cpsd`` and ``frf`` default to
456 empty lists when ``None`` is passed, and the specification CPSD is
457 sanitised so that zero values are replaced with ``NaN``.
458 """
459 self.source_data = source_data
460 self._geometry = None
461 self.geometry = geometry
462 self.coordinate = coordinate
463 self.control_coordinate = control_coordinate
464 self.response_coordinate = response_coordinate
465 self.excitation_coordinate = excitation_coordinate
466 self._time_data_truncated = None
467 self._time_data = None
468 self._start_time = None
469 self.time_data = time_data
470 self.start_time = start_time
471 self.cpsd = cpsd if cpsd is not None else []
472 self.specification_cpsd = specification_cpsd
473 self.specification_warning_psd = specification_warning_psd
474 self.specification_abort_psd = specification_abort_psd
475 self.averages = averages
476 self.overlap = overlap
477 self.fft_lines = fft_lines
478 self.window = window
479 self.samples_per_frame = samples_per_frame
480 self.units = units
481 self.frf = frf if frf is not None else []
482 self.oct_order = oct_order
484 # Sanitize the specification so 0 values go to NaN
485 if self.specification_cpsd is not None:
486 zero_mask = np.abs(self.specification_cpsd.ordinate)==0
487 self.specification_cpsd.ordinate[zero_mask] = np.nan
489 @property
490 def geometry(self):
491 """Geometry object for 3-D visualisation.
493 Returns
494 -------
495 Geometry or None
496 The loaded geometry object, or ``None`` if none has been set.
497 """
498 return self._geometry
500 @geometry.setter
501 def geometry(self, arg):
502 """Set the geometry from a file path, ``Geometry`` object, or ``None``.
504 Parameters
505 ----------
506 arg : str, Geometry, or None
507 If a string, the geometry is loaded from that file path. If a
508 ``Geometry`` instance it is stored directly. Any other value
509 (including ``None``) is stored as-is.
510 """
511 if isinstance(arg, str):
512 self._geometry = Geometry.load(arg)
513 elif isinstance(arg, Geometry):
514 self._geometry = arg
515 else:
516 self._geometry = arg
518 @property
519 def time_data(self):
520 """List of full-length ``TimeHistoryArray`` objects, one per data set.
522 Returns
523 -------
524 list of TimeHistoryArray or None
525 """
526 return self._time_data
528 @time_data.setter
529 def time_data(self, value):
530 """Set time data and invalidate any cached truncated arrays.
532 Parameters
533 ----------
534 value : list of TimeHistoryArray or None
535 New time data list.
536 """
537 self._time_data = value
538 self._time_data_truncated = None
540 @property
541 def start_time(self):
542 """List of analysis-window start times in seconds, one per data set.
544 A value of ``None`` for any element triggers automatic selection via
545 :func:`optimal_subset` when the truncated data is first requested.
547 Returns
548 -------
549 list of float or None
550 """
551 return self._start_time
553 @start_time.setter
554 def start_time(self, value):
555 """Set start times and invalidate any cached truncated arrays.
557 Parameters
558 ----------
559 value : list of float or None
560 New start-time list.
561 """
562 self._start_time = value
563 self._time_data_truncated = None
565 @property
566 def time_data_truncated(self):
567 """Time data windowed to the analysis subset used for CPSD computation.
569 Lazily computed on first access and cached. Each element is extracted
570 from the corresponding full ``TimeHistoryArray`` beginning at
571 ``start_time`` and covering exactly
572 ``samples_per_frame * averages - overlap_samples * (averages - 1)``
573 samples. When ``start_time`` is ``None`` for a data set, the optimal
574 start index is found via :func:`optimal_subset` and the corresponding
575 ``start_time`` entry is updated in-place.
577 Returns
578 -------
579 list of TimeHistoryArray
580 One truncated ``TimeHistoryArray`` per data set.
581 """
582 if self._time_data_truncated is None:
583 result = []
584 if self.start_time is None:
585 self.start_time = [None]*len(self.time_data)
586 for index, (start_time, time_data) in enumerate(zip(self.start_time, self.time_data)):
587 avg_dt_per_channel = np.mean(np.diff(time_data.abscissa, axis=1), axis=1)
588 dt = np.mean(avg_dt_per_channel)
589 samples_overlap = int(self.samples_per_frame * self.overlap)
590 samples_total = self.samples_per_frame * self.averages - samples_overlap * (self.averages - 1)
591 if start_time is None:
592 start_index = optimal_subset(time_data, samples_total)
593 self._start_time[index] = start_index * dt
594 result.append(time_data.extract_elements_by_abscissa(self._start_time[index], self._start_time[index] + samples_total * dt))
595 self._time_data_truncated = result
596 return self._time_data_truncated
598 def set_tolerance_limit_psd(self, dB=6):
599 """Set the abort-tolerance PSD from the specification using a dB tolerance.
601 Computes lower and upper bounds as
602 ``specification_asd / 10^(dB/10)`` and
603 ``specification_asd * 10^(dB/10)`` respectively, then stores them
604 as ``self.specification_abort_psd``. Only the abort band is updated;
605 ``specification_warning_psd`` is not modified.
607 Parameters
608 ----------
609 dB : float, optional
610 Tolerance in decibels above and below the specification.
611 Default is 6 dB.
613 Raises
614 ------
615 TypeError
616 If ``specification_cpsd`` is not a ``PowerSpectralDensityArray``.
617 """
618 if isinstance(self.specification_cpsd, PowerSpectralDensityArray):
619 lower_limit = self.specification_cpsd.get_asd() / 10**(dB / 10)
620 upper_limit = self.specification_cpsd.get_asd() * 10**(dB / 10)
621 self.specification_abort_psd = np.concatenate((lower_limit[np.newaxis, :], upper_limit[np.newaxis, :]))
622 else:
623 raise TypeError("specification_cpsd must be a PowerSpectralDensityArray to set tolerance limits")
625 def compute_cpsd(self):
626 """Compute cross-power spectral densities for each time data set.
628 Computes ASDs for each truncated time data segment and appends them to
629 ``self.cpsd``. Already-computed entries (``PowerSpectralDensityArray``
630 instances) are left unchanged; only ``None`` or non-array entries are
631 recomputed.
632 """
633 # Ensure the list is long enough
634 while len(self.cpsd) < len(self.time_data):
635 self.cpsd.append(None)
637 for index, time_data_truncated in enumerate(self.time_data_truncated):
638 if not isinstance(self.cpsd[index], PowerSpectralDensityArray):
639 self.cpsd[index] = time_data_truncated.cpsd(
640 samples_per_frame=self.samples_per_frame,
641 overlap=self.overlap,
642 window=self.window.lower(),
643 averages_to_keep=self.averages,
644 only_asds=True
645 )
647 def plot_cpsd_time_subset(self, save_path=None):
648 """Plot the selected time subset used for CPSD computation.
650 Parameters
651 ----------
652 save_path : str or None, optional
653 If provided, the figure is saved to this path. If there are
654 multiple time data sets the data set number is appended to the
655 filename stem. Default is ``None`` (figure is not saved).
656 """
657 for index in range(len(self.time_data)):
658 with plt.rc_context({'path.simplify': True, 'path.simplify_threshold': 1.0,
659 'agg.path.chunksize': 10000}):
660 fig, ax = plt.subplots(figsize=(12, 6))
661 ordinate_rms = np.sqrt(np.mean(self.time_data[index].ordinate**2, axis=0))
662 ordinate_rms_truncated = np.sqrt(np.mean(self.time_data_truncated[index].ordinate**2, axis=0))
663 ax.plot(self.time_data[index][0].abscissa, ordinate_rms, label='Full Time Set', color='darkgrey')
664 ax.plot(self.time_data_truncated[index][0].abscissa, ordinate_rms_truncated, label='Selected Time Subset', color='deepskyblue')
665 ax.set_title('Subsection of Time Data')
666 ax.set_xlabel('Time (s)')
667 ax.set_ylabel('RMS Amplitude Across All Channels')
668 ax.legend()
669 plt.tight_layout()
671 if save_path is not None:
672 if not os.path.exists(os.path.dirname(save_path)):
673 os.makedirs(os.path.dirname(save_path))
674 if len(self.time_data) == 1:
675 save_name = save_path
676 else:
677 save_name = os.path.splitext(save_path)[0] + ' (Data Set ' + str(index + 1) + ')' + os.path.splitext(save_path)[-1]
678 fig.savefig(save_name, dpi=600, bbox_inches='tight')
680 def compute_frf(self):
681 """Compute FRFs for each time data set and append them to ``self.frf``.
683 For each truncated time data segment, computes H1 frequency response
684 functions between the excitation (reference) channels and all response
685 channels using the same windowing parameters as CPSD computation.
686 Each resulting ``TransferFunctionArray`` is appended to ``self.frf``.
688 Excitation channels are identified by a non-empty ``feedback_channel``
689 entry in the source channel table, matching the logic used in
690 ``RattlesnakeData.excitation_coordinate``.
691 """
692 for index, time_data_truncated in enumerate(self.time_data_truncated):
693 reference_time_data = time_data_truncated[self.excitation_coordinate[:, np.newaxis]]
694 response_time_data = time_data_truncated[self.response_coordinate[:, np.newaxis]]
695 self.frf.append(sdpy.data.frf_from_time_data(
696 reference_data=reference_time_data,
697 response_data=response_time_data,
698 samples_per_average=self.samples_per_frame,
699 overlap=self.overlap,
700 method='H1',
701 window=self.window.lower()
702 ))
704 def plot_control_vs_spec(self, save_path: None|str = None):
705 """Plot measured control PSD against the specification for each control channel.
707 For each data set, creates a grid of subplots — one per control channel
708 — showing the abort-limit band (filled grey), warning-limit band
709 (light grey), specification (dashed black), and measured response
710 (blue). Out-of-limit regions are shaded and a percentage-out-of-limit
711 annotation is added to each subplot. The in-band RMS level is shown
712 as a text annotation.
714 Parameters
715 ----------
716 save_path : str or None, optional
717 Base file path for saving figures. The figure index and (when
718 multiple data sets exist) data set index are appended to the stem.
719 Default is ``None`` (figures are not saved).
721 Returns
722 -------
723 figures : list of list of matplotlib.figure.Figure
724 ``figures[i]`` is the list of figures for data set *i*.
725 axes : list of list of list of matplotlib.axes.Axes
726 ``axes[i]`` is the list of per-figure axes groups for data set *i*.
727 ``axes[i][j]`` is the list of ``Axes`` objects on figure *j*.
728 """
729 figures = [None]*len(self.time_data)
730 axes = [None]*len(self.time_data)
732 # If the CPSD is not defined, compute it automatically
733 if any(not isinstance(cpsd, PowerSpectralDensityArray) for cpsd in self.cpsd):
734 self.compute_cpsd()
736 control_coordinates = np.concatenate((self.control_coordinate[:,np.newaxis],self.control_coordinate[:,np.newaxis]),axis=1)
737 for index,cpsd in enumerate(self.cpsd):
738 xlim = NDDataArray.get_abscissa_limits([cpsd[control_coordinates],self.specification_cpsd.get_drive_points(),self.specification_warning_psd,self.specification_abort_psd])
739 ylim = NDDataArray.get_ordinate_limits([cpsd[control_coordinates],self.specification_cpsd.get_drive_points(),self.specification_warning_psd,self.specification_abort_psd],xlim)
741 np.max(np.abs(cpsd[control_coordinates].extract_elements_by_abscissa(min(xlim),max(xlim)).ordinate))
743 # Plot CPSD Results for Control Channels
744 all_figs, all_axes = make_multi_figures(len(self.control_coordinate))
746 coordinate_index = 0
747 figure_index = 0
748 for fig, axes_group in zip(all_figs,all_axes):
749 for ax in axes_group:
750 coordinate = self.control_coordinate[coordinate_index]
751 handles = []
752 labels = []
754 # Plot Abort Limit
755 if isinstance(self.specification_abort_psd,PowerSpectralDensityArray) and np.any(~np.isnan(self.specification_abort_psd.ordinate)):
756 abscissa = self.specification_abort_psd[0][coordinate].extract_elements_by_abscissa(min(xlim),max(xlim)).abscissa
757 ordinate1 = self.specification_abort_psd[-1][coordinate].extract_elements_by_abscissa(min(xlim),max(xlim)).ordinate
758 ordinate2 = self.specification_abort_psd[-0][coordinate].extract_elements_by_abscissa(min(xlim),max(xlim)).ordinate
759 ax.fill_between(abscissa,np.abs(ordinate1),np.abs(ordinate2), alpha=1, facecolor='lightgrey',edgecolor='darkgrey',linewidth=1)
760 ax.collections[-1].set_label('Abort Limit')
762 # Plot Warning Limit
763 if isinstance(self.specification_warning_psd,PowerSpectralDensityArray) and np.any(~np.isnan(self.specification_warning_psd.ordinate)):
764 abscissa = self.specification_warning_psd[0][coordinate].extract_elements_by_abscissa(min(xlim),max(xlim)).abscissa
765 ordinate1 = self.specification_warning_psd[-1][coordinate].extract_elements_by_abscissa(min(xlim),max(xlim)).ordinate
766 ordinate2 = self.specification_warning_psd[-0][coordinate].extract_elements_by_abscissa(min(xlim),max(xlim)).ordinate
767 ax.fill_between(abscissa,np.abs(ordinate1),np.abs(ordinate2), alpha=0.2, facecolor='lightgrey',edgecolor='darkgrey',linewidth=1)
768 ax.collections[-1].set_label('Warning Limit')
770 # Plot Specification
771 if isinstance(self.specification_cpsd,PowerSpectralDensityArray) and np.any(~np.isnan(self.specification_cpsd.ordinate)):
772 abscissa = self.specification_cpsd[coordinate].extract_elements_by_abscissa(min(xlim),max(xlim)).abscissa
773 ordinate = self.specification_cpsd[coordinate].extract_elements_by_abscissa(min(xlim),max(xlim)).ordinate
774 ax.plot(abscissa,np.abs(ordinate),color='black',linestyle='--',linewidth=1)
775 ax.get_lines()[-1].set_label('Specification')
777 # Plot Control
778 if isinstance(cpsd,PowerSpectralDensityArray):
779 abscissa = cpsd[coordinate].extract_elements_by_abscissa(min(xlim),max(xlim)).abscissa
780 ordinate = cpsd[coordinate].extract_elements_by_abscissa(min(xlim),max(xlim)).ordinate
781 ax.plot(abscissa,np.abs(ordinate),color='deepskyblue',linewidth=1)
782 ax.get_lines()[-1].set_label('Control')
784 # Shade Area Outside the Abort Limit if Response Exceeds the Abort Limit
785 if isinstance(cpsd,PowerSpectralDensityArray) and isinstance(self.specification_abort_psd,PowerSpectralDensityArray) and np.any(~np.isnan(self.specification_abort_psd.ordinate)):
786 response_too_low_mask = response_outside_mask(cpsd[coordinate],self.specification_abort_psd[0][coordinate],mask_type='below')
787 response_too_high_mask = response_outside_mask(cpsd[coordinate],self.specification_abort_psd[1][coordinate],mask_type='above')
789 _shade_out_of_limit(ax, cpsd[coordinate].abscissa, response_too_low_mask,
790 self.specification_abort_psd[0][coordinate].ordinate, min(ylim),
791 color='blue', alpha=0.2, label='Response Below Tolerance')
792 _shade_out_of_limit(ax, cpsd[coordinate].abscissa, response_too_high_mask,
793 self.specification_abort_psd[-1][coordinate].ordinate, max(ylim),
794 color='red', alpha=0.2, label='Response Above Tolerance')
796 response_inside_bandwidth_mask = np.logical_and(cpsd[coordinate].abscissa>=min(xlim),cpsd[coordinate].abscissa<=max(xlim))
797 response_outside_tolerance_mask = np.logical_or(response_too_low_mask,response_too_high_mask)
798 response_inside_bandwith_and_outside_tolerance = np.logical_and(response_inside_bandwidth_mask,response_outside_tolerance_mask)
799 percent_out = sum(response_inside_bandwith_and_outside_tolerance)/sum(response_inside_bandwidth_mask)*100
801 ax.text(0.99, 0.01, 'Out: ' + str(percent_out.round(2)) + '%',transform=ax.transAxes,ha="right", va="bottom", fontsize=8)
803 # Compute RMS Level of Response within the x-limits
804 rms_level = cpsd[coordinate].extract_elements_by_abscissa(min(xlim),max(xlim)).rms(self.oct_order)
806 # Set Plot Properties
807 ax.set_yscale('log')
809 ax.set_title(str(coordinate))
810 units = self.units[np.intersect1d(self.coordinate,coordinate,return_indices=True)[1]][0]
811 ax.text(0.01, 0.01, 'RMS: ' + str(rms_level.round(2)) + ' ' + units,transform=ax.transAxes,ha="left", va="bottom", fontsize=8)
812 ax.set_xlim(min(xlim),max(xlim))
813 ax.set_ylim(min(ylim),max(ylim))
814 ax.set_axisbelow(True)
816 # ax.xaxis.set_minor_locator(LogLocator(base=10.0, subs='auto', numticks=12))
817 ax.yaxis.set_minor_locator(LogLocator(base=10.0, subs='auto', numticks=12))
818 ax.grid(True,which='major',ls='-',color='grey',zorder=0, linewidth=0.2)
819 ax.grid(True,which='minor',ls=':',color='grey',zorder=0, linewidth=0.2)
821 coordinate_index += 1
823 h, l = ax.get_legend_handles_labels()
824 handles.extend(h)
825 labels.extend(l)
827 unique = dict(zip(labels,handles))
829 fig.legend(unique.values(),unique.keys(),loc='lower center',ncol=min([len(unique),3]))
830 fig.suptitle('Control vs. Specification PSD')
832 # Save Picture
833 if save_path is not None:
834 if not os.path.exists(os.path.dirname(save_path)):
835 os.makedirs(os.path.dirname(save_path))
836 if len(self.time_data) == 1:
837 save_name = os.path.splitext(save_path)[0] + ' (Figure ' + str(figure_index + 1) + ')' + os.path.splitext(save_path)[-1]
838 else:
839 save_name = os.path.splitext(save_path)[0] + ' (Data Set ' + str(index + 1) + ', Figure ' + str(figure_index + 1) + ')' + os.path.splitext(save_path)[-1]
840 fig.savefig(save_name,dpi=600 ,bbox_inches='tight')
842 figure_index += 1
843 figures[index] = all_figs
844 axes[index] = all_axes
845 return figures, axes
847 def plot_response(self, save_path: None|str = None, one_figure=False):
848 """Plot the measured response PSD for each response channel.
850 For each data set, creates a grid of subplots — one per response
851 channel — showing the measured ASD (blue) on a log scale. The
852 in-band RMS level is annotated on each subplot.
854 Parameters
855 ----------
856 save_path : str or None, optional
857 Base file path for saving figures. The figure index and (when
858 multiple data sets exist) data set index are appended to the stem.
859 Default is ``None`` (figures are not saved).
861 one_figure : bool, optional
862 Whether successive datasets should be overlayed on the same grid of
863 figures, or new grids should be made for each dataset
864 Default is `False`
866 Returns
867 -------
868 figures : list of list of matplotlib.figure.Figure
869 ``figures[i]`` is the list of figures for data set *i*.
870 axes : list of list of list of matplotlib.axes.Axes
871 ``axes[i]`` is the nested list of axis groups for data set *i*.
872 """
873 figures = [None]*len(self.time_data)
874 axes = [None]*len(self.time_data)
876 # If the CPSD is not defined, compute it automatically
877 if all(not isinstance(cpsd, PowerSpectralDensityArray) for cpsd in self.cpsd):
878 self.compute_cpsd()
880 response_coordinates = np.concatenate((self.response_coordinate[:,np.newaxis],self.response_coordinate[:,np.newaxis]),axis=1)
881 for index,cpsd in enumerate(self.cpsd):
882 xlim = NDDataArray.get_abscissa_limits(cpsd[response_coordinates])
883 ylim = NDDataArray.get_ordinate_limits(cpsd[response_coordinates],xlim)
885 # Plot CPSD Results for Response Channels
886 if (index == 0 and one_figure) or (not one_figure):
887 all_figs, all_axes = make_multi_figures(len(self.response_coordinate))
889 coordinate_index = 0
890 figure_index = 0
891 for fig, axes_group in zip(all_figs,all_axes):
892 for ax in axes_group:
893 coordinate = self.response_coordinate[coordinate_index]
894 # for fig in all_figs:
895 handles = []
896 labels = []
898 # Plot Response
899 if isinstance(cpsd,PowerSpectralDensityArray):
900 abscissa = cpsd[coordinate].extract_elements_by_abscissa(min(xlim),max(xlim)).abscissa
901 ordinate = cpsd[coordinate].extract_elements_by_abscissa(min(xlim),max(xlim)).ordinate
902 ax.plot(abscissa,np.abs(ordinate),color=None if one_figure else 'deepskyblue',linewidth=1)
904 # Compute RMS Level of Response within the x-limits
905 rms_level = cpsd[coordinate].extract_elements_by_abscissa(min(xlim),max(xlim)).rms(self.oct_order)
907 # Set Plot Properties
908 ax.set_yscale('log')
909 ax.set_title(str(coordinate))
910 units = self.units[np.intersect1d(self.coordinate,coordinate,return_indices=True)[1]][0]
911 if not one_figure:
912 ax.text(0.01, 0.01, 'RMS: ' + str(rms_level.round(2)) + ' ' + units,transform=ax.transAxes,ha="left", va="bottom")
914 ax.set_xlim(min(xlim),max(xlim))
915 ax.set_ylim(min(ylim),max(ylim))
916 ax.set_axisbelow(True)
917 ax.yaxis.set_minor_locator(LogLocator(base=10.0, subs='auto', numticks=12))
918 ax.grid(True,which='major',ls='-',color='grey',zorder=0, linewidth=0.2)
919 ax.grid(True,which='minor',ls=':',color='grey',zorder=0, linewidth=0.2)
921 coordinate_index += 1
923 fig.suptitle('Response PSD')
925 # Save Picture
926 if save_path is not None:
927 if (one_figure and index==len(self.cpsd)-1) or (not one_figure):
928 if not os.path.exists(os.path.dirname(save_path)):
929 os.makedirs(os.path.dirname(save_path))
930 if len(self.time_data) == 1 or one_figure:
931 save_name = os.path.splitext(save_path)[0] + ' (Figure ' + str(figure_index + 1) + ')' + os.path.splitext(save_path)[-1]
932 else:
933 save_name = os.path.splitext(save_path)[0] + ' (Data Set ' + str(index + 1) + ', Figure ' + str(figure_index + 1) + ')' + os.path.splitext(save_path)[-1]
934 fig.savefig(save_name,dpi=600 ,bbox_inches='tight')
936 figure_index += 1
937 figures[index] = all_figs
938 axes[index] = all_axes
939 return figures, axes
941 def plot_percent_lines_out(self, save_path: None|str = None):
942 """Plot the percentage of spectral lines out of tolerance for each control channel.
944 For each data set, computes the fraction of frequency lines within the
945 analysis bandwidth where the measured response falls outside the abort
946 tolerance limits, then displays the result as a horizontal bar chart.
947 Bars exceeding a 10 % threshold are coloured red. The fraction of
948 channels that exceed the threshold is shown in the figure title.
950 Parameters
951 ----------
952 save_path : str or None, optional
953 Base file path for saving figures. The data set index is appended
954 to the stem when multiple data sets exist. Default is ``None``
955 (figures are not saved).
957 Returns
958 -------
959 figures : list of matplotlib.figure.Figure
960 One figure per data set.
961 axes : list of list of matplotlib.axes.Axes
962 One list of axes per data set.
963 """
964 figures = [None]*len(self.time_data)
965 axes = [None]*len(self.time_data)
967 # If the CPSD is not defined, compute it automatically
968 if all(not isinstance(cpsd, PowerSpectralDensityArray) for cpsd in self.cpsd):
969 self.compute_cpsd()
971 control_coordinates = np.concatenate((self.control_coordinate[:,np.newaxis],self.control_coordinate[:,np.newaxis]),axis=1)
973 for index,cpsd in enumerate(self.cpsd):
974 xlim = NDDataArray.get_abscissa_limits([cpsd[control_coordinates],self.specification_cpsd.get_drive_points(),self.specification_warning_psd,self.specification_abort_psd])
976 # Plot Percent Lines Out
977 response_too_low_mask = response_outside_mask(cpsd[control_coordinates],self.specification_abort_psd[0][control_coordinates],mask_type='below')
978 response_too_high_mask = response_outside_mask(cpsd[control_coordinates],self.specification_abort_psd[1][control_coordinates],mask_type='above')
980 response_inside_bandwidth_mask = np.logical_and(cpsd[control_coordinates].abscissa>=min(xlim),cpsd[control_coordinates].abscissa<=max(xlim))
981 response_outside_tolerance_mask = np.logical_or(response_too_low_mask,response_too_high_mask)
982 response_inside_bandwith_and_outside_tolerance = np.logical_and(response_inside_bandwidth_mask,response_outside_tolerance_mask)
983 percent_out = sum(response_inside_bandwith_and_outside_tolerance.T)/sum(response_inside_bandwidth_mask.T)*100
985 threshold = 10
987 percent_channels_out = sum(percent_out>threshold)/len(percent_out)*100
989 fig,axes_group = dynamic_barh(values=percent_out,labels=self.control_coordinate)
990 for ax in axes_group:
991 bars = ax.patches
992 for bar in bars:
993 if bar.get_width() > threshold:
994 bar.set_facecolor('red')
995 ax.axvspan(0, threshold, alpha=0.6, facecolor='lightgrey',edgecolor='darkgrey',linewidth=1,zorder=0)
996 ax.set_xlabel('Percent Lines Out')
997 fig.suptitle('Percent Lines Out (' + str(percent_channels_out.round(2)) + r'% of Channels)')
998 fig.tight_layout()
1000 # Save Picture
1001 if save_path is not None:
1002 if not os.path.exists(os.path.dirname(save_path)):
1003 os.makedirs(os.path.dirname(save_path))
1004 if len(self.time_data) == 1:
1005 save_name = save_path
1006 else:
1007 save_name = os.path.splitext(save_path)[0] + ' (Data Set ' + str(index + 1) + ')' + os.path.splitext(save_path)[-1]
1008 fig.savefig(save_name,dpi=600 ,bbox_inches='tight')
1010 figures[index] = fig
1011 axes[index] = axes_group
1012 return figures, axes
1014 def plot_rms_level(self, save_path: None|str = None):
1015 """Plot the broadband RMS level for each response channel.
1017 RMS values are computed from the CPSD restricted to the analysis
1018 bandwidth defined by the specification, warning, and abort PSDs.
1019 Results are displayed as a horizontal bar chart.
1021 Parameters
1022 ----------
1023 save_path : str or None, optional
1024 Base file path for saving figures. The data set index is appended
1025 to the stem when multiple data sets exist. Default is ``None``
1026 (figures are not saved).
1028 Returns
1029 -------
1030 figures : list of matplotlib.figure.Figure
1031 One figure per data set.
1032 axes : list of list of matplotlib.axes.Axes
1033 One list of axes per data set.
1034 """
1035 figures = [None]*len(self.time_data)
1036 axes = [None]*len(self.time_data)
1038 # If the CPSD is not defined, compute it automatically
1039 if all(not isinstance(cpsd, PowerSpectralDensityArray) for cpsd in self.cpsd):
1040 self.compute_cpsd()
1042 response_coordinates = np.concatenate((self.response_coordinate[:,np.newaxis],self.response_coordinate[:,np.newaxis]),axis=1)
1043 for index,cpsd in enumerate(self.cpsd):
1044 xlim = NDDataArray.get_abscissa_limits([cpsd,self.specification_cpsd.get_drive_points(),self.specification_warning_psd,self.specification_abort_psd])
1046 # Get Response Response Level
1047 rms_response = cpsd[response_coordinates].extract_elements_by_abscissa(min(xlim),max(xlim)).rms(oct_order=self.oct_order)
1049 # Plot RMS Response Level
1050 fig,axis = dynamic_barh(values=rms_response,labels=self.response_coordinate)
1052 for ax in axis:
1053 ax.set_xlabel('RMS Level (EU)')
1054 fig.suptitle('RMS Response Level')
1055 fig.tight_layout()
1057 # Save Picture
1058 if save_path is not None:
1059 if not os.path.exists(os.path.dirname(save_path)):
1060 os.makedirs(os.path.dirname(save_path))
1061 if len(self.time_data) == 1:
1062 save_name = save_path
1063 else:
1064 save_name = os.path.splitext(save_path)[0] + ' (Data Set ' + str(index + 1) + ')' + os.path.splitext(save_path)[-1]
1065 fig.savefig(save_name,dpi=600 ,bbox_inches='tight')
1067 figures[index] = fig
1068 axes[index] = axis
1069 return figures,axes
1071 def plot_rms_error(self, save_path: None|str = None):
1072 """Plot the broadband RMS error between response and specification for each control channel.
1074 The RMS error is computed as the difference between the measured RMS
1075 level and the specification RMS level, both restricted to the analysis
1076 bandwidth. Results are displayed as a horizontal bar chart.
1078 Parameters
1079 ----------
1080 save_path : str or None, optional
1081 Base file path for saving figures. The data set index is appended
1082 to the stem when multiple data sets exist. Default is ``None``
1083 (figures are not saved).
1085 Returns
1086 -------
1087 figures : list of matplotlib.figure.Figure
1088 One figure per data set.
1089 axes : list of list of matplotlib.axes.Axes
1090 One list of axes per data set.
1091 """
1092 figures = [None]*len(self.time_data)
1093 axes = [None]*len(self.time_data)
1095 control_coordinates = np.concatenate((self.control_coordinate[:,np.newaxis],self.control_coordinate[:,np.newaxis]),axis=1)
1097 # If the CPSD is not defined, compute it automatically
1098 if all(not isinstance(cpsd, PowerSpectralDensityArray) for cpsd in self.cpsd):
1099 self.compute_cpsd()
1101 for index,cpsd in enumerate(self.cpsd):
1102 xlim = NDDataArray.get_abscissa_limits([cpsd[control_coordinates],self.specification_cpsd.get_drive_points(),self.specification_warning_psd,self.specification_abort_psd])
1104 # Get Response PSDs
1105 response_psds = cpsd[control_coordinates].extract_elements_by_abscissa(min(xlim),max(xlim))
1106 specification_psds = self.specification_cpsd[control_coordinates].extract_elements_by_abscissa(min(xlim),max(xlim))
1108 # Integrate the PSDs over the frequency band to get power in EU
1109 rms_response = response_psds.rms(self.oct_order)
1110 rms_spec = specification_psds.rms(self.oct_order)
1112 # Compute the RMS error between the two signals
1113 rms_error = rms_response - rms_spec
1115 # Plot RMS Error
1116 fig,axis = dynamic_barh(values=rms_error,labels=self.control_coordinate)
1118 for ax in axis:
1119 ax.set_xlabel('RMS Error (EU)')
1120 fig.suptitle('RMS Error Per Channel')
1121 fig.tight_layout()
1123 # Save Picture
1124 if save_path is not None:
1125 if not os.path.exists(os.path.dirname(save_path)):
1126 os.makedirs(os.path.dirname(save_path))
1127 if len(self.time_data) == 1:
1128 save_name = save_path
1129 else:
1130 save_name = os.path.splitext(save_path)[0] + ' (Data Set ' + str(index + 1) + ')' + os.path.splitext(save_path)[-1]
1131 fig.savefig(save_name,dpi=600 ,bbox_inches='tight')
1133 figures[index] = fig
1134 axes[index] = axis
1135 return figures,axes
1137 def plot_kurtosis(self, save_path: None|str = None):
1138 """Plot the kurtosis of each response channel time history.
1140 Kurtosis is computed on the truncated (analysis-window) portion of
1141 each data set using Fisher's definition disabled (Pearson kurtosis,
1142 where Gaussian has kurtosis = 3). Bars that fall outside the range
1143 ``[3 - 1, 3 + 1]`` are coloured red.
1145 Parameters
1146 ----------
1147 save_path : str or None, optional
1148 Base file path for saving figures. The data set index is appended
1149 to the stem when multiple data sets exist. Default is ``None``
1150 (figures are not saved).
1152 Returns
1153 -------
1154 figures : list of matplotlib.figure.Figure
1155 One figure per data set.
1156 axes : list of list of matplotlib.axes.Axes
1157 One list of axes per data set.
1158 """
1159 figures = [None]*len(self.time_data)
1160 axes = [None]*len(self.time_data)
1161 for index,time_data_truncated in enumerate(self.time_data_truncated):
1162 # Remove Excitation Coordinate Data to leave just response Data
1163 control_indices = time_data_truncated.response_coordinate.find_indices(self.excitation_coordinate)[0][0]
1164 time_data_truncated_indices = np.arange(len(time_data_truncated.response_coordinate))
1165 response_time_data = time_data_truncated[np.setdiff1d(time_data_truncated_indices,control_indices)]
1167 # Plot Kurtosis
1168 response_kurtosis = kurtosis(response_time_data.ordinate.T,fisher=False,bias=True)
1170 fig,axis = dynamic_barh(values=response_kurtosis,labels=self.response_coordinate)
1172 figures[index] = fig
1173 axes[index] = axis
1175 tolerance = 1
1177 for ax in axis:
1178 bars = ax.patches
1179 for bar in bars:
1180 if bar.get_width() > 3+tolerance or bar.get_width() < 3-tolerance:
1181 bar.set_facecolor('red')
1182 ax.axvspan(3-tolerance, 3+tolerance, alpha=0.6, facecolor='lightgrey',edgecolor='darkgrey',linewidth=1,zorder=0)
1183 ax.set_xlabel('Kurtosis')
1184 ax.set_xlim([0,3+tolerance])
1185 fig.suptitle('Kurtosis')
1186 fig.tight_layout()
1188 # Save Picture
1189 if save_path is not None:
1190 if not os.path.exists(os.path.dirname(save_path)):
1191 os.makedirs(os.path.dirname(save_path))
1192 if len(self.time_data) == 1:
1193 save_name = save_path
1194 else:
1195 save_name = os.path.splitext(save_path)[0] + ' (Data Set ' + str(index + 1) + ')' + os.path.splitext(save_path)[-1]
1196 fig.savefig(save_name,dpi=600 ,bbox_inches='tight')
1198 return figures,axes
1200 def plot_signal_to_noise(self, save_path: None|str = None):
1201 """Plot the signal-to-noise ratio for a two-dataset (noise + signal) test.
1203 Requires exactly two data sets: the first is treated as the noise
1204 measurement and the second as the signal measurement. SNR is computed
1205 as ``10 * log10(signal_power / noise_power)`` where power is the
1206 integral of the ASD over the full frequency range. Bars with SNR
1207 below 10 dB are coloured red.
1209 Parameters
1210 ----------
1211 save_path : str or None, optional
1212 File path for saving the figure. Default is ``None`` (figure is
1213 not saved).
1215 Returns
1216 -------
1217 fig : matplotlib.figure.Figure or None
1218 The figure, or ``None`` if ``time_data`` does not contain exactly
1219 two data sets.
1220 axis : list of matplotlib.axes.Axes or None
1221 The axes list, or ``None`` as above.
1222 """
1223 if len(self.time_data) == 2:
1224 # If the CPSD is not defined, compute it automatically
1225 if all(not isinstance(cpsd, PowerSpectralDensityArray) for cpsd in self.cpsd):
1226 self.compute_cpsd()
1228 # Calculate spacing along the frequency axis (last axis)
1229 diffs = np.diff(self.cpsd[-1].abscissa, axis=-1) # Shape: (M, N-1)
1231 # Create the inner widths (averages of adjacent gaps)
1232 # Shape: (M, N-2)
1233 inner_df = (diffs[:, :-1] + diffs[:, 1:]) / 2
1235 # Concatenate the boundary conditions
1236 # First col: same as first diff
1237 # Last col: same as last diff
1238 df = np.column_stack([
1239 diffs[:, [0]], # Twice half-distance to next
1240 inner_df, # Centered inner points
1241 diffs[:, [-1]] # Twice half-distance from previous
1242 ])
1244 # Calculate Signal Power
1245 signal_power = np.sum(self.cpsd[-1].ordinate.real*df,axis=1)
1246 noise_power = np.sum(self.cpsd[0].ordinate.real*df,axis=1)
1248 # Calculate Signal to Noise Ratio and Convert to dB
1249 signal_to_noise_ratios = 10*np.log10(signal_power / noise_power)
1251 fig,axis = dynamic_barh(values=signal_to_noise_ratios,labels=self.coordinate)
1253 for ax in axis:
1254 bars = ax.patches
1255 for bar in bars:
1256 if bar.get_width() < 10:
1257 bar.set_facecolor('red')
1258 ax.axvspan(10, max(signal_to_noise_ratios), alpha=0.6, facecolor='lightgrey',edgecolor='darkgrey',linewidth=1,zorder=0)
1259 ax.set_xlabel('Signal to Noise Ratio (dB)')
1260 ax.set_xlim(0,max(signal_to_noise_ratios))
1261 fig.suptitle('Signal to Noise Ratio')
1262 fig.tight_layout()
1264 # Save Picture
1265 if save_path is not None:
1266 if not os.path.exists(os.path.dirname(save_path)):
1267 os.makedirs(os.path.dirname(save_path))
1268 save_name = save_path
1269 fig.savefig(save_name,dpi=600 ,bbox_inches='tight')
1271 return fig,axis
1273 def plot_time_histories(self, save_path: None|str = None):
1274 """Plot the truncated time history for each channel.
1276 For each data set, creates a grid of subplots — one per channel in
1277 the truncated ``TimeHistoryArray`` — showing the time-domain signal.
1278 Path simplification is enabled for performance on large records.
1280 Parameters
1281 ----------
1282 save_path : str or None, optional
1283 Base file path for saving figures. The figure index and (when
1284 multiple data sets exist) data set index are appended to the stem.
1285 Default is ``None`` (figures are not saved).
1287 Returns
1288 -------
1289 figures : list of list of matplotlib.figure.Figure
1290 ``figures[i]`` is the list of figures for data set *i*.
1291 axes : list of list of list of matplotlib.axes.Axes
1292 ``axes[i]`` is the nested list of axis groups for data set *i*.
1293 """
1294 figures = [None]*len(self.time_data)
1295 axes = [None]*len(self.time_data)
1297 for index,time_data_truncated in enumerate(self.time_data_truncated):
1298 xlim = NDDataArray.get_abscissa_limits(time_data_truncated)
1299 ylim = NDDataArray.get_ordinate_limits(time_data_truncated,xlim)
1301 # Plot Time Data for Response Channels — enable path simplification for large signals
1302 with plt.rc_context({'path.simplify': True, 'path.simplify_threshold': 1.0,
1303 'agg.path.chunksize': 10000}):
1304 all_figs, all_axes = make_multi_figures(len(time_data_truncated),xlabel='Time (s)',ylabel='EU')
1306 coordinate_index = 0
1307 figure_index = 0
1308 for fig, axes_group in zip(all_figs,all_axes):
1309 for ax in axes_group:
1310 coordinate = np.squeeze(time_data_truncated.coordinate)[coordinate_index]
1311 handles = []
1312 labels = []
1314 # Plot Response
1315 ax.plot(time_data_truncated[coordinate].abscissa,time_data_truncated[coordinate].ordinate,color='deepskyblue',linewidth=1)
1317 units = self.units[np.intersect1d(self.coordinate,coordinate,return_indices=True)[1]][0]
1318 ax.text(0.99, 0.01, 'EU: ' + units,transform=ax.transAxes,ha="right", va="bottom")
1320 # Set Plot Properties
1321 ax.set_title(str(coordinate))
1322 ax.set_xlim(min(xlim),max(xlim))
1323 ax.set_ylim(-max(np.abs(ylim)),max(np.abs(ylim)))
1325 coordinate_index += 1
1327 fig.suptitle('Time Histories')
1329 # Save Picture
1330 if save_path is not None:
1331 if not os.path.exists(os.path.dirname(save_path)):
1332 os.makedirs(os.path.dirname(save_path))
1333 if len(self.time_data) == 1:
1334 save_name = os.path.splitext(save_path)[0] + ' (Figure ' + str(figure_index + 1) + ')' + os.path.splitext(save_path)[-1]
1335 else:
1336 save_name = os.path.splitext(save_path)[0] + ' (Data Set ' + str(index + 1) + ', Figure ' + str(figure_index + 1) + ')' + os.path.splitext(save_path)[-1]
1337 fig.savefig(save_name,dpi=600 ,bbox_inches='tight')
1339 figure_index += 1
1340 figures[index] = all_figs
1341 axes[index] = all_axes
1342 return figures, axes
1344 def plot_multiple_coherence(self, save_path: None|str = None):
1345 """Plot multiple coherence for each response channel within the analysis bandwidth.
1347 For each data set, computes the multiple coherence between every
1348 response channel and all excitation channels (using the same windowing
1349 parameters as CPSD computation), truncates to the specification
1350 bandwidth, then produces a grid of subplots with one panel per
1351 response channel. The y-axis is fixed to ``[0, 1]``.
1353 Parameters
1354 ----------
1355 save_path : str or None, optional
1356 Base file path for saving figures. The figure index and (when
1357 multiple data sets exist) data set index are appended to the stem.
1358 Default is ``None`` (figures are not saved).
1360 Returns
1361 -------
1362 figures : list of list of matplotlib.figure.Figure
1363 ``figures[i]`` is the list of figures for data set *i*.
1364 axes : list of list of list of matplotlib.axes.Axes
1365 ``axes[i]`` is the nested list of axis groups for data set *i*.
1366 """
1368 # If the CPSD is not defined, compute it automatically
1369 if all(not isinstance(cpsd, PowerSpectralDensityArray) for cpsd in self.cpsd):
1370 self.compute_cpsd()
1372 figures = [None]*len(self.time_data)
1373 axes = [None]*len(self.time_data)
1375 for index,(time_data_truncated,cpsd) in enumerate(zip(self.time_data_truncated,self.cpsd)):
1376 if self.specification_cpsd is not None:
1377 xlim = NDDataArray.get_abscissa_limits([self.specification_cpsd.get_drive_points(),self.specification_warning_psd,self.specification_abort_psd])
1378 else:
1379 xlim = NDDataArray.get_abscissa_limits(self.cpsd)
1381 # Get Time Data
1382 time_data_response = time_data_truncated[self.response_coordinate[:,np.newaxis]]
1383 time_data_reference = time_data_truncated[self.excitation_coordinate[:,np.newaxis]]
1385 # Compute Coherence
1386 multiple_coherence = sdpy.data.MultipleCoherenceArray.from_time_data(response_data=time_data_response,samples_per_average=self.samples_per_frame,overlap=self.overlap,window=self.window.lower(),reference_data=time_data_reference)
1388 # Truncate Frequency to Bandwith Analyzed
1389 if not any(lim is None for lim in xlim):
1390 multiple_coherence = multiple_coherence.extract_elements_by_abscissa(min(xlim),max(xlim))
1392 # Plot Multiple Coherence Results for Response Channels
1393 all_figs, all_axes = make_multi_figures(len(self.response_coordinate),ylabel='Coherence')
1395 coordinate_index = 0
1396 figure_index = 0
1397 for fig, axes_group in zip(all_figs,all_axes):
1398 for ax in axes_group:
1399 coordinate = self.response_coordinate[coordinate_index]
1401 # Plot Multiple Coherence
1402 abscissa = multiple_coherence[coordinate].abscissa
1403 ordinate = multiple_coherence[coordinate].ordinate
1404 ax.plot(abscissa,ordinate,color='deepskyblue',linewidth=1)
1406 # Set Plot Properties
1407 ax.set_title(str(coordinate))
1408 if any(lim is None for lim in xlim):
1409 ax.set_xlim(min(abscissa),max(abscissa))
1410 else:
1411 ax.set_xlim(min(xlim),max(xlim))
1412 ax.set_ylim(0,1)
1414 coordinate_index += 1
1416 fig.suptitle('Multiple Coherence')
1418 # Save Picture
1419 if save_path is not None:
1420 if not os.path.exists(os.path.dirname(save_path)):
1421 os.makedirs(os.path.dirname(save_path))
1422 if len(self.time_data) == 1:
1423 save_name = os.path.splitext(save_path)[0] + ' (Figure ' + str(figure_index + 1) + ')' + os.path.splitext(save_path)[-1]
1424 else:
1425 save_name = os.path.splitext(save_path)[0] + ' (Data Set ' + str(index + 1) + ', Figure ' + str(figure_index + 1) + ')' + os.path.splitext(save_path)[-1]
1426 fig.savefig(save_name,dpi=600 ,bbox_inches='tight')
1428 figure_index += 1
1429 figures[index] = all_figs
1430 axes[index] = all_axes
1431 return figures, axes
1433 def create_all_plots(self, figure_root_path='Figures'):
1434 """Generate and save all standard diagnostic plots to disk.
1436 Computes CPSDs if they have not already been computed, saves a time-
1437 subset plot, then iterates over all standard plot methods, saving each
1438 figure to a labelled subdirectory under *figure_root_path*. All
1439 matplotlib figures are closed after each method to free memory.
1441 The plot methods called (and their subfolder names) are:
1443 - ``plot_control_vs_spec`` -> ``Control vs Specification PSD Plot``
1444 - ``plot_response`` -> ``Response PSD Plot``
1445 - ``plot_percent_lines_out`` -> ``Percent Lines Out Plot``
1446 - ``plot_rms_level`` -> ``RMS Plot``
1447 - ``plot_rms_error`` -> ``RMS Error Plot``
1448 - ``plot_kurtosis`` -> ``Kurtosis Plot``
1449 - ``plot_time_histories`` -> ``Time Histories Plot``
1450 - ``plot_multiple_coherence`` -> ``Multiple Coherence Plot``
1451 - ``plot_signal_to_noise`` -> ``Signal to Noise Plot``
1453 Parameters
1454 ----------
1455 figure_root_path : str, optional
1456 Root directory under which all figure subdirectories are created.
1457 Default is ``'Figures'``.
1458 """
1459 # If the CPSD is not defined, compute it automatically
1460 if all(not isinstance(cpsd, PowerSpectralDensityArray) for cpsd in self.cpsd):
1461 self.compute_cpsd()
1462 self.plot_cpsd_time_subset(save_path=figure_root_path + os.sep + 'Analyzed Section of Time' + os.sep + 'Analyzed Section of Time Plot')
1464 # Specify Folder Names for Each Set of Plots
1465 # Keys are Method Names
1466 # Values are a sub-dictionary with keys 'plot_folder' and 'plot_name'
1467 path_and_plot_names = {
1468 'plot_control_vs_spec':{'plot_folder':'Control vs Specification PSD Plot','plot_name':'Control vs Specification PSD Plot'},
1469 'plot_response':{'plot_folder':'Response PSD Plot','plot_name':'Response PSD Plot'},
1470 'plot_percent_lines_out':{'plot_folder':'Percent Lines Out Plot','plot_name':'Percent Lines Out Plot'},
1471 'plot_rms_level':{'plot_folder':'RMS Plot','plot_name':'RMS Plot'},
1472 'plot_rms_error':{'plot_folder':'RMS Error Plot','plot_name':'RMS Error Plot'},
1473 'plot_kurtosis':{'plot_folder':'Kurtosis Plot','plot_name':'Kurtosis Plot'},
1474 'plot_time_histories':{'plot_folder':'Time Histories Plot','plot_name':'Time Histories Plot'},
1475 'plot_multiple_coherence':{'plot_folder':'Multiple Coherence Plot','plot_name':'Multiple Coherence Plot'},
1476 'plot_signal_to_noise':{'plot_folder':'Signal to Noise Plot','plot_name':'Signal to Noise Plot'},
1477 }
1479 # Create Plots for Each Plot Method
1480 for method_name in path_and_plot_names:
1481 method = getattr(self,method_name)
1482 method(save_path = figure_root_path + os.sep + path_and_plot_names[method_name]['plot_folder'] + os.sep + path_and_plot_names[method_name]['plot_name'])
1483 plt.close('all')
1485 def convert_to_octave(self, oct_order = 6):
1486 """Convert all stored PSDs to octave-band averages.
1488 Restricts the analysis to the specification bandwidth, computes nominal
1489 octave-band centres for the given *oct_order*, then replaces the stored
1490 CPSDs, specification CPSD, warning PSD, and abort PSD with their
1491 octave-band equivalents using ``bandwidth_average``. Sets
1492 ``self.oct_order`` so that subsequent RMS computations use octave
1493 integration.
1495 Parameters
1496 ----------
1497 oct_order : int, optional
1498 Octave-band order (e.g. 3 for third-octave, 6 for sixth-octave).
1499 Default is 6.
1500 """
1501 self.oct_order = oct_order
1503 # Compute CPSDs if they don't already exsist
1504 if all(not isinstance(cpsd, PowerSpectralDensityArray) for cpsd in self.cpsd):
1505 self.compute_cpsd()
1507 # Get Analysis Bandwidth
1508 xlim = NDDataArray.get_abscissa_limits([self.specification_cpsd.get_drive_points(),self.specification_warning_psd,self.specification_abort_psd])
1510 # Get Frequency Spacing
1511 nominal_band_centers,band_lb,band_ub,band_centers = sdpy.cpsd.nth_octave_freqs(freq=xlim,oct_order=oct_order)
1513 # Convert to Octave Space
1514 for index,cpsd in enumerate(self.cpsd):
1515 self.cpsd[index] = cpsd.extract_elements_by_abscissa(min(xlim),max(xlim)).bandwidth_average(band_lb,band_ub)
1517 self.specification_cpsd = self.specification_cpsd.extract_elements_by_abscissa(min(xlim),max(xlim)).bandwidth_average(band_lb,band_ub)
1518 self.specification_warning_psd = self.specification_warning_psd.extract_elements_by_abscissa(min(xlim),max(xlim)).bandwidth_average(band_lb,band_ub)
1519 self.specification_abort_psd = self.specification_abort_psd.extract_elements_by_abscissa(min(xlim),max(xlim)).bandwidth_average(band_lb,band_ub)
1521 @classmethod
1522 def load_rattlesnake_streaming_data(cls, file: os.PathLike | RattlesnakeData,
1523 environment_name: str = None,
1524 # time_array_index: int = 0,
1525 geometry: os.PathLike | Geometry = None):
1526 """Construct a ``RandomVibTest`` from a Rattlesnake streaming nc4 file.
1528 Reads (or reuses) a ``RattlesnakeData`` object, selects the
1529 requested environment, and populates a new ``RandomVibTest`` with all
1530 time data, coordinate arrays, specification PSDs, and spectral
1531 processing parameters drawn from the file.
1533 Parameters
1534 ----------
1535 file : os.PathLike or RattlesnakeData
1536 Path to a Rattlesnake nc4 file, or an already-parsed
1537 ``RattlesnakeData`` object.
1538 environment_name : str, optional
1539 Name of the environment group to load. When ``None`` (default)
1540 the first environment in the file is used. A warning is emitted
1541 when the file contains more than one environment and none is
1542 specified.
1543 geometry : os.PathLike or Geometry, optional
1544 Path to a geometry file or a ``Geometry`` object for 3-D
1545 visualisation. Default is ``None``.
1547 Returns
1548 -------
1549 RandomVibTest
1550 Fully initialised ``RandomVibTest`` instance ready for CPSD
1551 computation and plotting.
1553 Warns
1554 -----
1555 UserWarning
1556 If the file contains multiple environments and *environment_name*
1557 is not specified.
1558 """
1559 if not isinstance(file, RattlesnakeData):
1560 data = RattlesnakeData.read_rattlesnake_nc4(file)
1561 else:
1562 data = file
1564 if environment_name is None:
1565 environment_name = next(iter(data.environments))
1566 if len(data.environments) > 1:
1567 warnings.warn('There are multiple environments in the Rattlesnake data and none were specified. Using the first environmnet named: ' + environment_name + ' .')
1569 return cls(
1570 source_data = data.environments[environment_name],
1571 geometry = geometry,
1572 coordinate = data.get_coordinate(env_name=environment_name),
1573 control_coordinate = data.environments[environment_name].control_coordinate,
1574 response_coordinate = data.response_coordinate,
1575 excitation_coordinate = data.excitation_coordinate,
1576 time_data = data.get_time_data(env_name=environment_name),
1577 start_time = [None]*len(data.get_time_data(env_name=environment_name)),
1578 specification_cpsd = data.environments[environment_name].specification_cpsd,
1579 specification_warning_psd = data.environments[environment_name].specification_warning_psd,
1580 specification_abort_psd = data.environments[environment_name].specification_abort_psd,
1581 averages = data.environments[environment_name].sysid_averages,
1582 overlap = data.environments[environment_name].cpsd_overlap,
1583 fft_lines = data.environments[environment_name].fft_lines,
1584 window = data.environments[environment_name].cpsd_window,
1585 samples_per_frame = data.environments[environment_name].samples_per_frame,
1586 cpsd = [],
1587 units = data.units(env_name=environment_name),
1588 oct_order = None,
1589 )