reduce pvss computation via find_peaks

endaq/batch/analyzer.py

@@ -321,7 +321,7 @@ def _PVSSResultantData(self):
             damp=0.05,
             mode="pvss",
             aggregate_axes=True,
-        )
+        )[['Resultant']]
 
         return pv
 

endaq/calc/shock.py

@@ -13,8 +13,7 @@
 import pandas as pd
 import scipy.signal
 
-from endaq.calc.stats import L2_norm
-from endaq.calc import utils
+from endaq.calc import utils, stats
 
 
 def _absolute_acceleration_coefficients(omega, Q, T):
@@ -410,11 +409,14 @@ def shock_spectrum(
     bins_per_octave: float = 12.0,
     damp: float = 0.05,
     mode: typing.Literal["srs", "pvss"] = "srs",
+    max_time: float = 2.0,
+    peak_threshold: float = 0.1,
     two_sided: bool = False,
     aggregate_axes: bool = False,
 ) -> pd.DataFrame:
     """
-    Calculate the shock spectrum of an acceleration signal.
+    Calculate the shock spectrum of an acceleration signal. Note this defaults to first find peak events, then compute
+    the spectrum on those peaks to speed up processing time.
 
     :param accel: the absolute acceleration `y"`
     :param freqs: the natural frequencies across which to calculate the spectrum,
@@ -428,6 +430,12 @@ def shock_spectrum(
 
         *  `'srs'` (default) specifies the Shock Response Spectrum (SRS)
         *  `'pvss'` specifies the Pseudo-Velocity Shock Spectrum (PVSS)
+    :param max_time: the maximum duration in seconds to compute the shock spectrum for, if the time duration is greater
+        than :py:func:`~endaq.calc.stats.find_peaks()` is used to find peak locations, then the shock spectrums at each
+        peak is calculated with :py:func:`~endaq.calc.shock.rolling_shock_spectrum()` with `max_time` defining the
+        `slice_width`. Set `max_time` to `None` to force the function to not do the peak finding.
+    :param peak_threshold: if the duration is greater than `max_time` all peaks that are greater than `peak_threshold`
+        will be calculated, and the aggregate max per frequency will be reported
     :param two_sided: whether to return for each frequency:
         both the maximum negative and positive shocks (`True`),
         or simply the maximum absolute shock (`False`; default)
@@ -451,6 +459,39 @@ def shock_spectrum(
     freqs = np.asarray(freqs)
     if freqs.ndim != 1:
         raise ValueError("target frequencies must be in a 1D-array")
+
+    # Check for total time, if too long compute only for peaks
+    dt = utils.sample_spacing(accel)
+    if max_time is None:
+        max_time = dt * len(accel) * 2
+    if dt * len(accel) > max_time:
+        rolling_srs = rolling_shock_spectrum(
+            accel,
+            damp=damp,
+            mode=mode,
+            add_resultant=aggregate_axes,
+            freqs=freqs,
+            indexes=stats.find_peaks(accel, time_distance=max_time, threshold_multiplier=peak_threshold),
+            slice_width=max_time
+        )
+        var_name = 'variable'
+        if 'axis' in rolling_srs.columns:
+            var_name = 'axis'
+        srs = pd.DataFrame(
+            index=pd.Series(rolling_srs['frequency (Hz)'].unique(), name='frequency (Hz)'),
+            columns=rolling_srs[var_name].unique()
+        )
+        for var in srs.columns:
+            pivot = rolling_srs[rolling_srs[var_name] == var].copy()
+            pivot = pivot.pivot_table(
+                columns='timestamp',
+                index='frequency (Hz)',
+                values='value'
+            )
+            srs[var] = pivot.max(axis=1)
+        srs.columns.name = accel.columns.name
+        return srs
+
     omega = 2 * np.pi * freqs
 
     if mode == "srs":
@@ -465,7 +506,6 @@ def shock_spectrum(
         dtype=np.float64,
     )
 
-    dt = utils.sample_spacing(accel)
     T_padding = 1 / (
         freqs.min() * np.sqrt(1 - damp ** 2)
     )  # uses lowest damped frequency
@@ -489,8 +529,8 @@ def shock_spectrum(
         )
 
         if aggregate_axes:
-            rd = L2_norm(rd, axis=-1, keepdims=True)
-            rd_padding = L2_norm(rd_padding, axis=-1, keepdims=True)
+            rd = stats.L2_norm(rd, axis=-1, keepdims=True)
+            rd_padding = stats.L2_norm(rd_padding, axis=-1, keepdims=True)
 
         results[(0,) + i_nd] = -np.minimum(rd.min(axis=0), rd_padding.min(axis=0))
         results[(1,) + i_nd] = np.maximum(rd.max(axis=0), rd_padding.max(axis=0))
@@ -499,7 +539,7 @@ def shock_spectrum(
         return pd.DataFrame(
             np.maximum(results[0], results[1]),
             index=pd.Series(freqs, name="frequency (Hz)"),
-            columns=(["resultant"] if aggregate_axes else accel.columns),
+            columns=(["Resultant"] if aggregate_axes else accel.columns),
         )
 
     return namedtuple("PseudoVelocityResults", "neg pos")(
@@ -517,6 +557,7 @@ def rolling_shock_spectrum(
         damp: float = 0.05,
         mode: typing.Literal["srs", "pvss"] = "srs",
         add_resultant: bool = True,
+        freqs: np.ndarray = None,
         init_freq: float = 0.5,
         bins_per_octave: float = 12.0,
         num_slices: int = 100,
@@ -536,9 +577,9 @@ def rolling_shock_spectrum(
         *  `'srs'` (default) specifies the Shock Response Spectrum (SRS)
         *  `'pvss'` specifies the Pseudo-Velocity Shock Spectrum (PVSS)
     :param add_resultant: if `True` (default) the column-wise resultant will
-        also be computed
-    :param add_resultant: whether to calculate the column-wise resultant (`True`)
-        or calculate spectra along each column independently (`False`; default)
+        also be computed and returned with the spectra per-column
+    :param freqs: the natural frequencies across which to calculate the spectrum,
+        if `None` (the default) it uses `init_freq` and `bins_per_octave` to define them
     :param init_freq: the initial frequency in the sequence; if `None`,
         use the frequency corresponding to the data's duration, default is 0.5 Hz
     :param bins_per_octave: the number of frequencies per octave, default is 12    
@@ -548,8 +589,6 @@ def rolling_shock_spectrum(
     :param index_values: the index values of each peak event to quantify (slower but more intuitive than using `indexes`)
     :param slice_width: the time in seconds to center about each slice index,
         if none is provided it will calculate one based upon the number of slices
-    :param add_resultant: if `True` the root sum of squares of each shock spectrum column will
-        also be computed, default is `False`
     :param disable_warnings: if `True` (default) it disables the warnings on the initial frequency
     :return: a dataframe containing all the shock spectrums, stacked on each other
 
@@ -558,6 +597,8 @@ def rolling_shock_spectrum(
     Surface plots, and Animations
 
     """
+    if freqs is None:
+        freqs = utils.logfreqs(df, init_freq=init_freq, bins_per_octave=bins_per_octave)
 
     indexes, slice_width, num, length = utils._rolling_slice_definitions(
         df,
@@ -579,8 +620,7 @@ def rolling_shock_spectrum(
                 df.iloc[window_start:window_end],
                 mode=mode,
                 damp=damp,
-                init_freq=init_freq,
-                bins_per_octave=bins_per_octave,
+                freqs=freqs,
             )
         if add_resultant:
             with warnings.catch_warnings():
@@ -590,10 +630,9 @@ def rolling_shock_spectrum(
                     df.iloc[window_start:window_end],
                     mode=mode,
                     damp=damp,
-                    init_freq=init_freq,
-                    bins_per_octave=bins_per_octave,
+                    freqs=freqs,
                     aggregate_axes=True
-                )['resultant']
+                )['Resultant']
 
         slice_srs = slice_srs.reset_index().melt(id_vars=slice_srs.index.name)
         slice_srs['timestamp'] = df.index[i]

endaq/calc/stats.py

@@ -214,7 +214,7 @@ def shock_vibe_metrics(
         pvss = shock.shock_spectrum(df, freqs=freqs, damp=damp, mode='pvss')
         if include_resultant:
             pvss['Resultant'] = shock.shock_spectrum(df, freqs=freqs, damp=damp, mode='pvss', aggregate_axes=True)[
-                                    'resultant']
+                                    'Resultant']
 
         if display_plots:
             px.line(pvss, log_x=True, log_y=True).update_layout(yaxis_title_text='Pseudo Velocity',
@@ -412,22 +412,23 @@ def rolling_metrics(
 
     Here's a continuation of the example shown in :py:func:`~endaq.calc.stats.find_peaks()`::
 
-    .. code:: python3
+    .. code-block:: python
 
-        #Calculate for all Peak Event Indexes
+        # Calculate for all Peak Event Indexes
         metrics = endaq.calc.stats.rolling_metrics(accel, indexes=indexes, slice_width=2.0)
 
-        #Calculate for 3 Specific Times
+        # Calculate for 3 Specific Times
         import pandas as pd
         metrics = endaq.calc.stats.rolling_metrics(
             accel,
             index_values = pd.DatetimeIndex(['2020-03-13 23:40:13', '2020-03-13 23:45:00', '2020-03-13 23:50:00'],tz='UTC'),
             slice_width=5.0)
 
-        #Calculate for 50 Equally Spaced & Sized Slices, Turning off Pseudo Velocity (Only Recommended for Smaller Time Slices)
+        # Calculate for 50 Equally Spaced & Sized Slices, Turning off Pseudo Velocity (Only Recommended for Smaller Time Slices)
         metrics = endaq.calc.stats.rolling_metrics(
             accel, num_slices=50, highpass_cutoff=2,
             tukey_percent=0.0, include_pseudo_velocity=False)
+
     """
     indexes, slice_width, num, length = utils._rolling_slice_definitions(
         df,