From 2703c56c24b87a3d0cb08e5b2ca67a679eb8ab72 Mon Sep 17 00:00:00 2001 From: shanlyMIDE <35080650+shanlyMIDE@users.noreply.github.com> Date: Thu, 7 Apr 2022 17:11:24 -0400 Subject: [PATCH 01/24] spectrum_over_time --- endaq/plot/plots.py | 300 +++++++++++++++++++++++++++++++++++++------- 1 file changed, 256 insertions(+), 44 deletions(-) diff --git a/endaq/plot/plots.py b/endaq/plot/plots.py index a4f1e53..d0d99eb 100644 --- a/endaq/plot/plots.py +++ b/endaq/plot/plots.py @@ -4,11 +4,13 @@ import pandas as pd import plotly.graph_objects as go import plotly.express as px +from plotly import colors import scipy.spatial.transform from plotly.subplots import make_subplots from scipy import signal, spatial from typing import Optional from collections.abc import Container +import datetime from endaq.calc import sample_spacing from endaq.calc.psd import to_octave, welch @@ -275,8 +277,8 @@ def octave_spectrogram(df: pd.DataFrame, window: float, bins_per_octave: int = 3 max_freq: float = float('inf'), db_scale: bool = True, log_scale_y_axis: bool = True ) -> go.Figure: """ - Produces an octave spectrogram of the given data. - + Produces an octave spectrogram of the given data, this is a wrapper around + :py:func:`~endaq.calc.psd.rolling_psd()` and :py:func:`~endaq.plot.plots.spectrum_over_time()` :param df: The dataframe of sensor data. This must only have 1 column. :param window: The time window for each of the columns in the spectrogram :param bins_per_octave: The number of frequency bins per octave @@ -285,55 +287,18 @@ def octave_spectrogram(df: pd.DataFrame, window: float, bins_per_octave: int = 3 :param db_scale: If the spectrogram should be log-scaled for visibility (with 10*log10(x)) :param log_scale_y_axis: If the y-axis of the plot should be log scaled :return: a tuple containing: - - - the frequency bins - - the time bins - - the spectrogram data + - dataframe of the spectrogram data - the corresponding plotly figure """ if len(df.columns) != 1: raise ValueError("The parameter 'df' must have only one column of data!") - - ary = df.values.squeeze() - - fs = 1/sample_spacing(df)#(len(df) - 1) / (df.index[-1] - df.index[0]) - N = int(fs * window) #Number of points in the fft - w = signal.blackman(N, False) - freqs, bins, Pxx = signal.spectrogram(ary, fs, window=w, nperseg=N, noverlap=0) + num_slices = int(len(df)*utils.sample_spacing(df) / window) - time_dfs = [pd.DataFrame({bins[j]: Pxx[:, j]}, index=freqs) for j in range(len(bins))] - - octave_dfs = list(map(lambda x: to_octave(x, freq_start, octave_bins=bins_per_octave, agg='sum'), time_dfs)) - - combined_df = pd.concat(octave_dfs, axis=1) + df_psd = rolling_psd(df, num_slices=num_slices, octave_bins = bins_per_octave, + fstart = freq_start) - freqs = combined_df.index.values - Pxx = combined_df.values - - include_freqs_mask = freqs <= max_freq - Pxx = Pxx[include_freqs_mask] - freqs = freqs[include_freqs_mask] - - if db_scale: - Pxx = 10 * np.log10(Pxx) - - trace = [go.Heatmap(x=bins, y=freqs, z=Pxx, colorscale='Jet')] - layout = go.Layout( - yaxis={'title': 'Frequency (Hz)'}, - xaxis={'title': 'Time (s)'}, - ) - - fig = go.Figure(data=trace, layout=layout) - - if log_scale_y_axis: - fig.update_yaxes(type="log") - - fig.update_traces(showscale=False) - - data_df = pd.DataFrame(Pxx, index=freqs, columns=bins) - - return data_df, fig + return df_psd, spectrum_over_time(df_psd, y_limit=max_freq, log_z = db_scale, log_y = log_scale_y_axis) def octave_psd_bar_plot(df: pd.DataFrame, bins_per_octave: int = 3, f_start: float = 20.0, yaxis_title: str = '', @@ -538,3 +503,250 @@ def animate_quaternion(df, rate=6., scale=1.): fig.layout.updatemenus[0].buttons[0].args[1]['frame']['duration'] = int(1000//rate) return fig + + +def spectrum_over_time( + df: pd.DataFrame, + plot_type: typing.Literal["Heatmap", "Surface", "Waterfall", "Animation"] = "Heatmap", + var_column: str = 'variable', + var_to_process: str = None, + x_column: str = 'timestamp', + y_column: str = 'frequency (Hz)', + z_column: str = 'value', + y_limit: float = None, + log_y: bool = False, + log_z: bool = False, + round_time: bool = True, + waterfall_line_sequence: bool = True, + waterfall_line_color: str = '#EE7F27', +) -> pd.DataFrame: + """ + Generate a 3D Plotly figure from a stacked spectrum to visualize how the frequency content changes over time + :param df: the input dataframe with columns defining the frequency content, timestamps, values, and variables, + see the following functions which provides outputs that would then be passed into this function as an input: + - :py:func:`~endaq.calc.fft.rolling_fft()` + - :py:func:`~endaq.calc.psd.rolling_psd()` + - :py:func:`~endaq.calc.shock.rolling_shock_spectrum()` + :param plot_type: the type of plot to display the spectrum, options are: + - `Heatmap`: a 2D visualization with the color defining the z value (the default) + - `Surface`: similar to `Heatmap` but the z value is also projected "off the page" + - `Waterfall`: distinct lines are plotted per time slice in a 3D view + - `Animation`: a 2D display of the waterfall but the spectrum changes with animation frames + :param var_column: the column name in the dataframe that defines the different variables, default is `"variable"` + :param var_to_process: the variable value in the `var_column` to filter the input df down to, + if none is provided (the default) this function will filter to the first value + :param x_column: the column name in the dataframe that defines the timestamps, default is `"timestamp"` + :param y_column: the column name in the dataframe that defines the frequency, default is `"frequency (Hz)"` + :param z_column: the column name in the dataframe that defines the values, default is `"value"` + :param y_limit: the limit of the y axis (frequency) to include in the figure, + default is None meaning it will display all content + :param log_y: if `True` the y axis (frequency) will be in a log scale, default is `False` + :param log_z: if `True` the z axis (values) will be in a log scale, default is `False` + :param round_time: if `True` (default) the time values will be rounded to the nearest second for datetimes and hundredths of a second for floats + :param waterfall_line_sequence: if `True` the waterfall line colors are defined with a color scale, + if `False` all lines will have the same color, default is `True` + :param waterfall_line_color: the color to use for all lines in the Waterfall plot if `waterfall_line_sequence` is `False` + :return: a Plotly figure visualizing the spectrum over time + + + Here's a few examples from a dataset recorded with an enDAQ sensor on a motorcycle + as it revved the engine which resulted in changing frequency content + .. code:: python + import endaq + endaq.plot.utilities.set_theme() + import pandas as pd + + #Get the Data + df_vibe = pd.read_csv('https://info.endaq.com/hubfs/data/motorcycle-vibration-moving-frequency.csv',index_col=0) + df_vibe = df_vibe - df_vibe.median() + + #Calculate a Rolling FFT + fft = endaq.calc.fft.rolling_fft(df_vibe, num_slices=200, add_resultant=True) + + #Visualize the Rolling FFT as a Heatmap + heatmap = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Heatmap', log_y=False, log_z=False, y_limit=200, var_to_process='Resultant') + heatmap.show() + + #Visualize as a Surface Plot + surface = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Surface', log_y=False, log_z=False, y_limit=200, var_to_process='Resultant') + surface.show() + + #Visualize as a Waterfall + waterfall = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Waterfall', log_y=False, log_z=False, y_limit=200, var_to_process='Resultant') + waterfall.show() + + #Visualize as an Animation + animation = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Animation', log_y=False, log_z=False, y_limit=200, var_to_process='Resultant') + animation.show() + + .. plotly:: + :fig-vars: heatmap, surface, waterfall, animation + + import endaq + endaq.plot.utilities.set_theme() + import pandas as pd + + #Get the Data + df_vibe = pd.read_csv('https://info.endaq.com/hubfs/data/motorcycle-vibration-moving-frequency.csv',index_col=0) + df_vibe = df_vibe - df_vibe.median() + + #Calculate a Rolling FFT + fft = endaq.calc.fft.rolling_fft(df_vibe, num_slices=200, add_resultant=True) + + #Visualize the Rolling FFT as a Heatmap + heatmap = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Heatmap', log_y=False, log_z=False, y_limit=200, var_to_process='Resultant') + heatmap.show() + + #Visualize as a Surface Plot + surface = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Surface', log_y=False, log_z=False, y_limit=200, var_to_process='Resultant') + surface.show() + + #Visualize as a Waterfall + waterfall = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Waterfall', log_y=False, log_z=False, y_limit=200, var_to_process='Resultant') + waterfall.show() + + #Visualize as an Animation + animation = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Animation', log_y=False, log_z=False, y_limit=200, var_to_process='Resultant') + animation.show() + + """ + df = df.copy() + + #Filter to one variable + if var_to_process is None: + var_to_process = df[var_column].unique()[0] + df = df.loc[df[var_column] == var_to_process] + + #Round time + if round_time: + if isinstance(df[x_column].iloc[0], datetime.datetime): + df[x_column] = df[x_column].round('s') + else: + df[x_column] = np.round(df[x_column].to_numpy(),2) + + #Filter frequency + if y_limit is not None: + df = df.loc[df[y_column] < y_limit] + + #Remove 0s + df = df.loc[df[z_column] > 0] + + #Check Length of Dataframe + if len(df) > 50000: + warnings.warn( + "plot data is very large, may be unresponsive", + RuntimeWarning, + ) + + #Heatmap & Surface + if plot_type in ["Heatmap", "Surface"]: + df[y_column] = np.round(df[y_column].to_numpy(),2) + df = df.pivot(columns=x_column,index=y_column,values=z_column) + + #Deal with Dateimes + x_type = float + if isinstance(df.columns[0], datetime.datetime): + if plot_type == "Surface": + df.columns = (df.columns - df.columns[0]).total_seconds() + else: + x_type = str + + #Build Dictionary of Plot Data, Apply Log Scale If Needed + data_dict = { + 'x': df.columns.astype(x_type), + 'y': df.index.astype(float), + 'z': df.to_numpy().astype(float), + 'connectgaps' : True + } + if log_z: + data_dict['z'] = np.log10(data_dict['z']) + + #Generate Figures + if plot_type == "Heatmap": + fig = go.Figure(data = go.Heatmap(data_dict, zsmooth='best')).update_layout( + xaxis_title_text='Timestamp', yaxis_title_text='Frequency (Hz)') + if log_y: + fig.update_layout(yaxis_type='log') + else: + fig = go.Figure(data = go.Surface(data_dict)) + fig.update_traces(showscale=False) + + #Waterfall + elif plot_type == 'Waterfall': + #Define Colors + if waterfall_line_sequence: + color_sequence = colors.sample_colorscale( + [[0.0, '#6914F0'], + [0.2, '#3764FF'], + [0.4, '#2DB473'], + [0.6, '#FAC85F'], + [0.8, '#EE7F27'], + [1.0, '#D72D2D']], len(df[x_column].unique())) + else: + color_sequence = [waterfall_line_color] + + #Deal with Datetime + df['label_column'] = df[x_column] + if isinstance(df[x_column].iloc[0], datetime.datetime): + df[x_column] = (df[x_column] - df[x_column].iloc[0]).dt.total_seconds() + df['label_column'] = df['label_column'].dt.tz_convert(None).astype(str) + + #Generate Figure + fig = px.line_3d( + df, + x = x_column, + y = y_column, + z = z_column, + color = 'label_column', + color_discrete_sequence = color_sequence).update_layout( + legend_orientation = 'v', + legend_y = 1, + legend_title_text = 'Timestamps' + ) + + #Animation + elif plot_type == 'Animation': + #Define Range (If Undefined it Starts with Fits to First Frame) + range = [df['value'].min(), df['value'].max()] + if log_z: + range = np.log10(range) + + #Deal with Datetime + if isinstance(df[x_column].iloc[0], datetime.datetime): + df[x_column] = df[x_column].dt.tz_convert(None).astype(str) + + #Gerate Figure + fig = px.line( + df, + animation_frame = x_column, + x = y_column, + y = z_column, + log_y = log_z, + log_x = log_y + ).update_layout( + showlegend = False, + yaxis_range = range, + yaxis_title_text = '', + xaxis_title_text = 'Frequency (Hz)' + ) + + else: + raise ValueError(f"invalid plot type {plot_type}") + + #Add Labels to 3D plots + if plot_type in ["Surface", "Waterfall"]: + fig.update_scenes( + aspectratio_x=2.0, + aspectratio_y=1.0, + aspectratio_z=0.3, + xaxis_title_text='Timestamp', + yaxis_title_text='Frequency (Hz)', + zaxis_title_text= '' + ) + if log_y: + fig.update_scenes(yaxis_type = 'log') + if log_z: + fig.update_scenes(zaxis_type = 'log') + fig.update_layout(scene_camera = dict(eye=dict(x=-1.5, y=-1.5, z=1))) + + return fig From 7f8957c0595dece43e6d8c9eafc1fa2bfd83283f Mon Sep 17 00:00:00 2001 From: shanlyMIDE <35080650+shanlyMIDE@users.noreply.github.com> Date: Thu, 7 Apr 2022 17:12:02 -0400 Subject: [PATCH 02/24] rolling_fft --- endaq/calc/fft.py | 73 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 73 insertions(+) diff --git a/endaq/calc/fft.py b/endaq/calc/fft.py index b6901d9..43b4f06 100644 --- a/endaq/calc/fft.py +++ b/endaq/calc/fft.py @@ -35,6 +35,79 @@ def aggregate_fft(df, **kwargs): return psd.welch(df, **kwargs) +def rolling_fft( + df: pd.DataFrame, + bin_width: float = 1.0, + num_slices: int = 100, + indexes: np.array = None, + index_values: np.array = None, + slice_width: float = None, + add_resultant: bool = True, + disable_warnings: bool = True, + **kwargs, +) -> pd.DataFrame: + """ + Compute FFTs for defined slices of a time series data set using :py:func:`~endaq.calc.fft.aggregate_fft()` + :param df: the input dataframe with an index defining the time in seconds or datetime + :param bin_width: the bin width or resolution in Hz for the FFT + :param num_slices: the number of slices to split the time series into, default is 100, + this is ignored if `indexes` is defined + :param indexes: the center index locations (not value) of each slice to compute the FFT + :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 FFT column will + also be computed + :param disable_warnings: if `True` (default) it disables the warnings on the PSD length + :param kwargs: Other parameters to pass directly to :py:func:`~endaq.calc.fft.aggregate_fft()` + :return: a dataframe containing all the FFTs, stacked on each other + + See example use cases and syntax at :py:func:`~endaq.plot.plots.spectrum_over_time()` + which visualizes the output of this function in Heatmaps, Waterfall plots, + Surface plots, and Animations + + """ + if disable_warnings: + warnings.filterwarnings('ignore', '.*greater than.*', ) + + length = len(df) + + #Define center index locations of each slice if not provided + if indexes is None: + if index_values is not None: + indexes = np.zeros(len(index_values),int) + for i in range(len(indexes)): + indexes[i] = int((np.abs(df.index - index_values[i])).argmin()) + else: + indexes = np.linspace(0, length, num_slices, endpoint = False, dtype=int) + indexes = indexes + int(indexes[1]/2) + + #Calculate slice step size + spacing = utils.sample_spacing(df) + if slice_width is None: + slice_width = spacing * length / len(indexes) + num = int(slice_width / spacing / 2) + + #Loop through and compute fft + fft = pd.DataFrame() + for i in indexes: + window_start = max(0, i - num) + window_end = min(length, i + num) + slice_fft = aggregate_fft( + df.iloc[window_start:window_end], + bin_width = bin_width, + **kwargs + ) + if add_resultant: + slice_fft['Resultant'] = slice_fft.pow(2).sum(axis=1).pow(0.5) + + slice_fft = slice_fft.reset_index().melt(id_vars=slice_fft.index.name) + slice_fft['timestamp'] = df.index[i] + fft = pd.concat([fft,slice_fft]) + + return fft + + def fft( df: pd.DataFrame, output: typing.Literal[None, "magnitude", "angle", "complex"] = None, From f80c01fef9343c8057220fb472091fdc8115565e Mon Sep 17 00:00:00 2001 From: shanlyMIDE <35080650+shanlyMIDE@users.noreply.github.com> Date: Thu, 7 Apr 2022 17:12:35 -0400 Subject: [PATCH 03/24] rolling_psd --- endaq/calc/psd.py | 86 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 86 insertions(+) diff --git a/endaq/calc/psd.py b/endaq/calc/psd.py index 6304ad8..7007c49 100644 --- a/endaq/calc/psd.py +++ b/endaq/calc/psd.py @@ -335,3 +335,89 @@ def vc_curves( # The PSD must already scale by ∆f -> need only scale by √∆f? # TODO make `density` parameter, scale differently depending on mode return np.sqrt(accel_psd.index[1]) * df_vel_oct.apply(np.sqrt) + + +def rolling_psd( + df: pd.DataFrame, + bin_width: float = 1.0, + octave_bins: float = None, + fstart: float = 1.0, + num_slices: int = 100, + indexes: np.array = None, + index_values: np.array = None, + slice_width: float = None, + add_resultant: bool = True, + disable_warnings: bool = True, + **kwargs, +) -> pd.DataFrame: + """ + Compute PSDs for defined slices of a time series data set using :py:func:`~endaq.calc.psd.welch()` + :param df: the input dataframe with an index defining the time in seconds or datetime + :param bin_width: the bin width or resolution in Hz for the PSD, defaults to 1, + this is ignored if `octave_bins` is defined + :param octave_bins: the number of frequency bins in each octave, defaults to `None` + :param fstart: the lowest frequency for an octave PSD, defaults to 1 + :param num_slices: the number of slices to split the time series into, default is 100, + this is ignored if `indexes` is defined + :param indexes: the center index locations (not value) of each slice to compute the PSD + :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` (default) the root sum of squares of each PSD column will + also be computed + :param disable_warnings: if `True` (default) it disables the warnings on the PSD length + :param kwargs: Other parameters to pass directly to :py:func:`psd.welch` + :return: a dataframe containing all the PSDs, stacked on each other + + See example use cases and syntax at :py:func:`~endaq.plot.plots.spectrum_over_time()` + which visualizes the output of this function in Heatmaps, Waterfall plots, + Surface plots, and Animations + + """ + if disable_warnings: + warnings.filterwarnings('ignore', '.*greater than.*', ) + warnings.filterwarnings('ignore', '.*empty frequency*', ) + + length = len(df) + + #Define center index locations of each slice if not provided + if indexes is None: + if index_values is not None: + indexes = np.zeros(len(index_values),int) + for i in range(len(indexes)): + indexes[i] = int((np.abs(df.index - index_values[i])).argmin()) + else: + indexes = np.linspace(0, length, num_slices, endpoint = False, dtype=int) + indexes = indexes + int(indexes[1]/2) + + #Calculate slice step size + spacing = utils.sample_spacing(df) + if slice_width is None: + slice_width = spacing * length / len(indexes) + num = int(slice_width / spacing / 2) + + #Define bin_width if octave spacing + if octave_bins is not None: + bin_width = 1 / slice_width + + #Loop through and compute PSD + psd = pd.DataFrame() + for i in indexes: + window_start = max(0, i - num) + window_end = min(length, i + num) + slice_psd = welch( + df.iloc[window_start:window_end], + bin_width = bin_width, + **kwargs + ) + if octave_bins is not None: + slice_psd = to_octave(slice_psd, octave_bins=octave_bins, fstart = fstart) + + if add_resultant: + slice_psd['Resultant'] = slice_psd.sum(axis=1) + + slice_psd = slice_psd.reset_index().melt(id_vars=slice_psd.index.name) + slice_psd['timestamp'] = df.index[i] + psd = pd.concat([psd,slice_psd]) + + return psd From dc56ab8f2a89f1d10699a8acbcd670cb060e750d Mon Sep 17 00:00:00 2001 From: shanlyMIDE <35080650+shanlyMIDE@users.noreply.github.com> Date: Thu, 7 Apr 2022 17:14:11 -0400 Subject: [PATCH 04/24] rolling_shock_spectrum --- endaq/calc/shock.py | 115 +++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 109 insertions(+), 6 deletions(-) diff --git a/endaq/calc/shock.py b/endaq/calc/shock.py index a354d0b..b60d31d 100644 --- a/endaq/calc/shock.py +++ b/endaq/calc/shock.py @@ -405,7 +405,9 @@ def relative_displacement_static(accel: pd.DataFrame, omega: float, damp: float def shock_spectrum( accel: pd.DataFrame, - freqs: np.ndarray, + freqs: np.ndarray = None, + init_freq: float = 0.5, + bins_per_octave: float = 12.0, damp: float = 0.0, mode: typing.Literal["srs", "pvss"] = "srs", two_sided: bool = False, @@ -413,13 +415,15 @@ def shock_spectrum( ) -> pd.DataFrame: """ Calculate the shock spectrum of an acceleration signal. - :param accel: the absolute acceleration `y"` - :param freqs: the natural frequencies across which to calculate the spectrum + :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 :param damp: the damping coefficient `ζ`, related to the Q-factor by `ζ = 1/(2Q)`; defaults to 0 :param mode: the type of spectrum to calculate: - - `'srs'` (default) specifies the Shock Response Spectrum (SRS) - `'pvss'` specifies the Pseudo-Velocity Shock Spectrum (PVSS) :param two_sided: whether to return for each frequency: @@ -428,9 +432,7 @@ def shock_spectrum( :param aggregate_axes: whether to calculate the column-wise resultant (`True`) or calculate spectra along each column independently (`False`; default) :return: the shock spectrum - .. seealso:: - - `Pseudo Velocity Shock Spectrum Rules For Analysis Of Mechanical Shock, Howard A. Gaberson `_ - `An Introduction To The Shock Response Spectrum, Tom Irvine, 9 July 2012 `_ - `SciPy transfer functions `_ @@ -440,6 +442,8 @@ def shock_spectrum( Documentation for the biquad function used to implement the transfer function. """ + if freqs is None: + freqs = utils.logfreqs(accel, init_freq = init_freq, bins_per_octave = bins_per_octave) if two_sided and aggregate_axes: raise ValueError("cannot enable both options `two_sided` and `aggregate_axes`") freqs = np.asarray(freqs) @@ -506,6 +510,105 @@ def shock_spectrum( ) +def rolling_shock_spectrum( + df: pd.DataFrame, + damp: float = 0.0, + mode: typing.Literal["srs", "pvss"] = "srs", + add_resultant: bool = True, + init_freq: float = 0.5, + bins_per_octave: float = 12.0, + num_slices: int = 100, + indexes: np.array = None, + index_values: np.array = None, + slice_width: float = None, + multiplier: float = 1.0, + disable_warnings: bool = True, +) -> pd.DataFrame: + """ + Compute Shock Response Spectrums for defined slices of a time series data set using :py:func:`~endaq.calc.shock.shock_spectrum()` + :param df: the input dataframe with an index defining the time in seconds or datetime + :param damp: the damping coefficient `ζ`, related to the Q-factor by + `ζ = 1/(2Q)`; defaults to 0 + :param mode: the type of spectrum to calculate: + - `'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 aggregate_axes: whether to calculate the column-wise resultant (`True`) + or calculate spectra along each column independently (`False`; default) + :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 + :param num_slices: the number of slices to split the time series into, default is 100, + this is ignored if `indexes` is defined + :param indexes: the center index locations (not value) of each slice to compute the 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 multiplier: applies a scale to the output + - 386.09 to convert from g to inches (in) + - 9806.65 to convert from g to millimeters (mm) + :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 + + See example use cases and syntax at :py:func:`~endaq.plot.plots.spectrum_over_time()` + which visualizes the output of this function in Heatmaps, Waterfall plots, + Surface plots, and Animations + + """ + if disable_warnings: + warnings.filterwarnings('ignore', '.*too short*', ) + + length = len(df) + + #Define center index locations of each slice if not provided + if indexes is None: + if index_values is not None: + indexes = np.zeros(len(index_values),int) + for i in range(len(indexes)): + indexes[i] = int((np.abs(df.index - index_values[i])).argmin()) + else: + indexes = np.linspace(0, length, num_slices, endpoint = False, dtype=int) + indexes = indexes + int(indexes[1]/2) + + #Calculate slice step size + spacing = utils.sample_spacing(df) + if slice_width is None: + slice_width = spacing * length / len(indexes) + num = int(slice_width / spacing / 2) + + #Loop through and compute shock spectrum + srs = pd.DataFrame() + for i in indexes: + window_start = max(0, i - num) + window_end = min(length, i + num) + slice_srs = shock_spectrum( + df.iloc[window_start:window_end] * multiplier, + mode = mode, + damp = damp, + init_freq = init_freq, + bins_per_octave = bins_per_octave, + ) + if add_resultant: + slice_srs['Resultant'] = shock_spectrum( + df.iloc[window_start:window_end] * multiplier, + mode = mode, + damp = damp, + init_freq = init_freq, + bins_per_octave = bins_per_octave, + aggregate_axes=True + )['resultant'] + + + slice_srs = slice_srs.reset_index().melt(id_vars=slice_srs.index.name) + slice_srs['timestamp'] = df.index[i] + srs = pd.concat([srs,slice_srs]) + + return srs + + @dataclass class HalfSineWavePulse: """ From e75b7b200c1896875f07dbf79e2356d2ef98109a Mon Sep 17 00:00:00 2001 From: shanlyMIDE <35080650+shanlyMIDE@users.noreply.github.com> Date: Thu, 7 Apr 2022 17:18:38 -0400 Subject: [PATCH 05/24] import warnings --- endaq/calc/fft.py | 1 + 1 file changed, 1 insertion(+) diff --git a/endaq/calc/fft.py b/endaq/calc/fft.py index 43b4f06..afdc17d 100644 --- a/endaq/calc/fft.py +++ b/endaq/calc/fft.py @@ -2,6 +2,7 @@ import typing from typing import Optional +import warnings import numpy as np import pandas as pd From bc2e39c93b4b67bba8c58d672631fa30dd4d21ff Mon Sep 17 00:00:00 2001 From: shanlyMIDE <35080650+shanlyMIDE@users.noreply.github.com> Date: Thu, 7 Apr 2022 17:23:30 -0400 Subject: [PATCH 06/24] octave_spectrogram fix --- endaq/plot/plots.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/endaq/plot/plots.py b/endaq/plot/plots.py index d0d99eb..2b0b93b 100644 --- a/endaq/plot/plots.py +++ b/endaq/plot/plots.py @@ -293,7 +293,7 @@ def octave_spectrogram(df: pd.DataFrame, window: float, bins_per_octave: int = 3 if len(df.columns) != 1: raise ValueError("The parameter 'df' must have only one column of data!") - num_slices = int(len(df)*utils.sample_spacing(df) / window) + num_slices = int(len(df)*sample_spacing(df) / window) df_psd = rolling_psd(df, num_slices=num_slices, octave_bins = bins_per_octave, fstart = freq_start) From d1edba612367c801ab39c21b66b62b83084f090f Mon Sep 17 00:00:00 2001 From: shanlyMIDE <35080650+shanlyMIDE@users.noreply.github.com> Date: Thu, 7 Apr 2022 17:31:08 -0400 Subject: [PATCH 07/24] spectrums_over_time added import --- endaq/plot/plots.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/endaq/plot/plots.py b/endaq/plot/plots.py index 2b0b93b..743f83e 100644 --- a/endaq/plot/plots.py +++ b/endaq/plot/plots.py @@ -13,7 +13,7 @@ import datetime from endaq.calc import sample_spacing -from endaq.calc.psd import to_octave, welch +from endaq.calc.psd import to_octave, welch, rolling_psd from .utilities import determine_plotly_map_zoom, get_center_of_coordinates from .dashboards import rolling_enveloped_dashboard From 5fce50e81c2a2bc62778dd7aa0a37a16edc6ee58 Mon Sep 17 00:00:00 2001 From: shanlyMIDE <35080650+shanlyMIDE@users.noreply.github.com> Date: Fri, 8 Apr 2022 09:42:47 -0400 Subject: [PATCH 08/24] rolling_psd_scaling --- endaq/calc/psd.py | 16 +++++++++++++++- 1 file changed, 15 insertions(+), 1 deletion(-) diff --git a/endaq/calc/psd.py b/endaq/calc/psd.py index 7007c49..695dd0f 100644 --- a/endaq/calc/psd.py +++ b/endaq/calc/psd.py @@ -342,6 +342,11 @@ def rolling_psd( bin_width: float = 1.0, octave_bins: float = None, fstart: float = 1.0, + scaling: typing.Literal[None, "density", "spectrum", "parseval", "unit"] = None, + agg: typing.Union[ + typing.Literal["mean", "sum"], + typing.Callable[[np.ndarray, int], float], + ] = "mean", num_slices: int = 100, indexes: np.array = None, index_values: np.array = None, @@ -357,6 +362,14 @@ def rolling_psd( this is ignored if `octave_bins` is defined :param octave_bins: the number of frequency bins in each octave, defaults to `None` :param fstart: the lowest frequency for an octave PSD, defaults to 1 + :param scaling: the scaling of the output; `"density"` & `"spectrum"` + correspond to the same options in ``scipy.signal.welch``; `"parseval"` + will maintain the "energy" between the input & output, s.t. + ``welch(df, scaling="parseval").sum(axis="rows")`` is roughly equal to + ``df.abs().pow(2).sum(axis="rows")``; + `"unit"` will maintain the units and scale of the input data. + :param agg: the method for aggregating values into bins (only used if converting to octave); `'mean'` preserves + the PSD's area-under-the-curve, `'sum'` preserves the PSD's "energy" :param num_slices: the number of slices to split the time series into, default is 100, this is ignored if `indexes` is defined :param indexes: the center index locations (not value) of each slice to compute the PSD @@ -408,10 +421,11 @@ def rolling_psd( slice_psd = welch( df.iloc[window_start:window_end], bin_width = bin_width, + scaling = scaling, **kwargs ) if octave_bins is not None: - slice_psd = to_octave(slice_psd, octave_bins=octave_bins, fstart = fstart) + slice_psd = to_octave(slice_psd, octave_bins=octave_bins, fstart = fstart, agg=agg) if add_resultant: slice_psd['Resultant'] = slice_psd.sum(axis=1) From 6611d9d52c68eb75fae47a4d0b18897ecb2b4a14 Mon Sep 17 00:00:00 2001 From: shanlyMIDE <35080650+shanlyMIDE@users.noreply.github.com> Date: Fri, 8 Apr 2022 10:33:52 -0400 Subject: [PATCH 09/24] Update fft.py --- endaq/calc/fft.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/endaq/calc/fft.py b/endaq/calc/fft.py index afdc17d..7f87dbf 100644 --- a/endaq/calc/fft.py +++ b/endaq/calc/fft.py @@ -11,7 +11,7 @@ from endaq.calc import psd, utils -__all__ = ["aggregate_fft", "fft", "rfft", "dct", "dst"] +__all__ = ["aggregate_fft", "rolling_fft", "fft", "rfft", "dct", "dst"] def aggregate_fft(df, **kwargs): @@ -49,6 +49,7 @@ def rolling_fft( ) -> pd.DataFrame: """ Compute FFTs for defined slices of a time series data set using :py:func:`~endaq.calc.fft.aggregate_fft()` + :param df: the input dataframe with an index defining the time in seconds or datetime :param bin_width: the bin width or resolution in Hz for the FFT :param num_slices: the number of slices to split the time series into, default is 100, From 60e29ce98d46ce107286ee079c470fe6b84978b9 Mon Sep 17 00:00:00 2001 From: shanlyMIDE <35080650+shanlyMIDE@users.noreply.github.com> Date: Fri, 8 Apr 2022 10:34:20 -0400 Subject: [PATCH 10/24] Update psd.py --- endaq/calc/psd.py | 1 + 1 file changed, 1 insertion(+) diff --git a/endaq/calc/psd.py b/endaq/calc/psd.py index 695dd0f..b2ea14e 100644 --- a/endaq/calc/psd.py +++ b/endaq/calc/psd.py @@ -357,6 +357,7 @@ def rolling_psd( ) -> pd.DataFrame: """ Compute PSDs for defined slices of a time series data set using :py:func:`~endaq.calc.psd.welch()` + :param df: the input dataframe with an index defining the time in seconds or datetime :param bin_width: the bin width or resolution in Hz for the PSD, defaults to 1, this is ignored if `octave_bins` is defined From 78c2b0e5b1b9f9fe2bcf64fb07812c0356102ae5 Mon Sep 17 00:00:00 2001 From: shanlyMIDE <35080650+shanlyMIDE@users.noreply.github.com> Date: Fri, 8 Apr 2022 10:35:14 -0400 Subject: [PATCH 11/24] Update shock.py --- endaq/calc/shock.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/endaq/calc/shock.py b/endaq/calc/shock.py index b60d31d..fe0395d 100644 --- a/endaq/calc/shock.py +++ b/endaq/calc/shock.py @@ -415,6 +415,7 @@ def shock_spectrum( ) -> pd.DataFrame: """ Calculate the shock spectrum of an acceleration signal. + :param accel: the absolute acceleration `y"` :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 @@ -526,6 +527,7 @@ def rolling_shock_spectrum( ) -> pd.DataFrame: """ Compute Shock Response Spectrums for defined slices of a time series data set using :py:func:`~endaq.calc.shock.shock_spectrum()` + :param df: the input dataframe with an index defining the time in seconds or datetime :param damp: the damping coefficient `ζ`, related to the Q-factor by `ζ = 1/(2Q)`; defaults to 0 From ba8c1b5b24c6b0ad2f5fa1790621287097882988 Mon Sep 17 00:00:00 2001 From: shanlyMIDE <35080650+shanlyMIDE@users.noreply.github.com> Date: Fri, 8 Apr 2022 10:36:31 -0400 Subject: [PATCH 12/24] Update plots.py --- endaq/plot/plots.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/endaq/plot/plots.py b/endaq/plot/plots.py index 743f83e..0bf15aa 100644 --- a/endaq/plot/plots.py +++ b/endaq/plot/plots.py @@ -28,6 +28,7 @@ 'rolling_min_max_envelope', 'around_peak', 'animate_quaternion', + 'spectrum_over_time' ] DEFAULT_ATTRIBUTES_TO_PLOT_INDIVIDUALLY = np.array([ @@ -522,6 +523,7 @@ def spectrum_over_time( ) -> pd.DataFrame: """ Generate a 3D Plotly figure from a stacked spectrum to visualize how the frequency content changes over time + :param df: the input dataframe with columns defining the frequency content, timestamps, values, and variables, see the following functions which provides outputs that would then be passed into this function as an input: - :py:func:`~endaq.calc.fft.rolling_fft()` @@ -551,6 +553,7 @@ def spectrum_over_time( Here's a few examples from a dataset recorded with an enDAQ sensor on a motorcycle as it revved the engine which resulted in changing frequency content + .. code:: python import endaq endaq.plot.utilities.set_theme() From b49a7eb49e68ef840038bbbee93e1470be017d0b Mon Sep 17 00:00:00 2001 From: shanlyMIDE <35080650+shanlyMIDE@users.noreply.github.com> Date: Fri, 8 Apr 2022 10:42:25 -0400 Subject: [PATCH 13/24] Update plots.py --- endaq/plot/plots.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/endaq/plot/plots.py b/endaq/plot/plots.py index 0bf15aa..0c0570a 100644 --- a/endaq/plot/plots.py +++ b/endaq/plot/plots.py @@ -279,7 +279,8 @@ def octave_spectrogram(df: pd.DataFrame, window: float, bins_per_octave: int = 3 ) -> go.Figure: """ Produces an octave spectrogram of the given data, this is a wrapper around - :py:func:`~endaq.calc.psd.rolling_psd()` and :py:func:`~endaq.plot.plots.spectrum_over_time()` + :py:func:`~endaq.calc.psd.rolling_psd()` and :py:func:`~endaq.plot.spectrum_over_time()` + :param df: The dataframe of sensor data. This must only have 1 column. :param window: The time window for each of the columns in the spectrogram :param bins_per_octave: The number of frequency bins per octave From 197777bf1518e636889058c7a9f649f1c02c3f82 Mon Sep 17 00:00:00 2001 From: shanlyMIDE <35080650+shanlyMIDE@users.noreply.github.com> Date: Fri, 8 Apr 2022 10:42:56 -0400 Subject: [PATCH 14/24] Update fft.py --- endaq/calc/fft.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/endaq/calc/fft.py b/endaq/calc/fft.py index 7f87dbf..9b0371d 100644 --- a/endaq/calc/fft.py +++ b/endaq/calc/fft.py @@ -64,7 +64,7 @@ def rolling_fft( :param kwargs: Other parameters to pass directly to :py:func:`~endaq.calc.fft.aggregate_fft()` :return: a dataframe containing all the FFTs, stacked on each other - See example use cases and syntax at :py:func:`~endaq.plot.plots.spectrum_over_time()` + See example use cases and syntax at :py:func:`~endaq.plot.spectrum_over_time()` which visualizes the output of this function in Heatmaps, Waterfall plots, Surface plots, and Animations From bee6993b59a17d976f64e10c1f96ee37d565c369 Mon Sep 17 00:00:00 2001 From: shanlyMIDE <35080650+shanlyMIDE@users.noreply.github.com> Date: Fri, 8 Apr 2022 10:44:53 -0400 Subject: [PATCH 15/24] Update psd.py --- endaq/calc/psd.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/endaq/calc/psd.py b/endaq/calc/psd.py index b2ea14e..b39c50b 100644 --- a/endaq/calc/psd.py +++ b/endaq/calc/psd.py @@ -383,7 +383,7 @@ def rolling_psd( :param kwargs: Other parameters to pass directly to :py:func:`psd.welch` :return: a dataframe containing all the PSDs, stacked on each other - See example use cases and syntax at :py:func:`~endaq.plot.plots.spectrum_over_time()` + See example use cases and syntax at :py:func:`~endaq.plot.spectrum_over_time()` which visualizes the output of this function in Heatmaps, Waterfall plots, Surface plots, and Animations From 637d48b5182f5870aee3a2b86b5b64949b779755 Mon Sep 17 00:00:00 2001 From: shanlyMIDE <35080650+shanlyMIDE@users.noreply.github.com> Date: Fri, 8 Apr 2022 10:46:36 -0400 Subject: [PATCH 16/24] Update shock.py --- endaq/calc/shock.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/endaq/calc/shock.py b/endaq/calc/shock.py index fe0395d..c5463a7 100644 --- a/endaq/calc/shock.py +++ b/endaq/calc/shock.py @@ -532,6 +532,7 @@ def rolling_shock_spectrum( :param damp: the damping coefficient `ζ`, related to the Q-factor by `ζ = 1/(2Q)`; defaults to 0 :param mode: the type of spectrum to calculate: + - `'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 @@ -550,12 +551,13 @@ def rolling_shock_spectrum( :param add_resultant: if `True` the root sum of squares of each shock spectrum column will also be computed, default is `False` :param multiplier: applies a scale to the output + - 386.09 to convert from g to inches (in) - 9806.65 to convert from g to millimeters (mm) :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 - See example use cases and syntax at :py:func:`~endaq.plot.plots.spectrum_over_time()` + See example use cases and syntax at :py:func:`~endaq.plot.spectrum_over_time()` which visualizes the output of this function in Heatmaps, Waterfall plots, Surface plots, and Animations From c6001e1a13c03c23f55ed976360a8607a7491ec8 Mon Sep 17 00:00:00 2001 From: shanlyMIDE <35080650+shanlyMIDE@users.noreply.github.com> Date: Fri, 8 Apr 2022 11:02:54 -0400 Subject: [PATCH 17/24] Update plots.py --- endaq/plot/plots.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/endaq/plot/plots.py b/endaq/plot/plots.py index 0c0570a..3dc1e26 100644 --- a/endaq/plot/plots.py +++ b/endaq/plot/plots.py @@ -521,7 +521,7 @@ def spectrum_over_time( round_time: bool = True, waterfall_line_sequence: bool = True, waterfall_line_color: str = '#EE7F27', -) -> pd.DataFrame: +) -> go.Figure: """ Generate a 3D Plotly figure from a stacked spectrum to visualize how the frequency content changes over time From 8e99dbb59d508d1bce841ecec63d261bbdba26e6 Mon Sep 17 00:00:00 2001 From: "MIDE\\shanly" Date: Fri, 8 Apr 2022 17:14:57 -0400 Subject: [PATCH 18/24] made changes Pycharm and Connor suggested (I think)! --- endaq/calc/fft.py | 47 ++++++++++++----------------- endaq/calc/psd.py | 60 +++++++++++++++++-------------------- endaq/calc/shock.py | 73 ++++++++++++++++++++------------------------- endaq/calc/utils.py | 46 ++++++++++++++++++++++++++++ 4 files changed, 125 insertions(+), 101 deletions(-) diff --git a/endaq/calc/fft.py b/endaq/calc/fft.py index 9b0371d..76b24b7 100644 --- a/endaq/calc/fft.py +++ b/endaq/calc/fft.py @@ -69,45 +69,36 @@ def rolling_fft( Surface plots, and Animations """ - if disable_warnings: - warnings.filterwarnings('ignore', '.*greater than.*', ) - length = len(df) - - #Define center index locations of each slice if not provided - if indexes is None: - if index_values is not None: - indexes = np.zeros(len(index_values),int) - for i in range(len(indexes)): - indexes[i] = int((np.abs(df.index - index_values[i])).argmin()) - else: - indexes = np.linspace(0, length, num_slices, endpoint = False, dtype=int) - indexes = indexes + int(indexes[1]/2) - - #Calculate slice step size - spacing = utils.sample_spacing(df) - if slice_width is None: - slice_width = spacing * length / len(indexes) - num = int(slice_width / spacing / 2) + indexes, slice_width, num, length = utils._rolling_slice_definitions( + df, + indexes=indexes, + index_values=index_values, + num_slices=num_slices, + slice_width=slice_width + ) - #Loop through and compute fft - fft = pd.DataFrame() + # Loop through and compute fft + fft_df = pd.DataFrame() for i in indexes: window_start = max(0, i - num) window_end = min(length, i + num) - slice_fft = aggregate_fft( - df.iloc[window_start:window_end], - bin_width = bin_width, - **kwargs - ) + with warnings.catch_warnings(): + if disable_warnings: + warnings.filterwarnings('ignore', '.*greater than.*') + slice_fft = aggregate_fft( + df.iloc[window_start:window_end], + bin_width=bin_width, + **kwargs + ) if add_resultant: slice_fft['Resultant'] = slice_fft.pow(2).sum(axis=1).pow(0.5) slice_fft = slice_fft.reset_index().melt(id_vars=slice_fft.index.name) slice_fft['timestamp'] = df.index[i] - fft = pd.concat([fft,slice_fft]) + fft_df = pd.concat([fft_df, slice_fft]) - return fft + return fft_df def fft( diff --git a/endaq/calc/psd.py b/endaq/calc/psd.py index b39c50b..ea532d6 100644 --- a/endaq/calc/psd.py +++ b/endaq/calc/psd.py @@ -388,51 +388,45 @@ def rolling_psd( Surface plots, and Animations """ - if disable_warnings: - warnings.filterwarnings('ignore', '.*greater than.*', ) - warnings.filterwarnings('ignore', '.*empty frequency*', ) - - length = len(df) - - #Define center index locations of each slice if not provided - if indexes is None: - if index_values is not None: - indexes = np.zeros(len(index_values),int) - for i in range(len(indexes)): - indexes[i] = int((np.abs(df.index - index_values[i])).argmin()) - else: - indexes = np.linspace(0, length, num_slices, endpoint = False, dtype=int) - indexes = indexes + int(indexes[1]/2) - - #Calculate slice step size - spacing = utils.sample_spacing(df) - if slice_width is None: - slice_width = spacing * length / len(indexes) - num = int(slice_width / spacing / 2) - - #Define bin_width if octave spacing + + indexes, slice_width, num, length = utils._rolling_slice_definitions( + df, + indexes=indexes, + index_values=index_values, + num_slices=num_slices, + slice_width=slice_width + ) + + # Define bin_width if octave spacing if octave_bins is not None: bin_width = 1 / slice_width - - #Loop through and compute PSD + + # Loop through and compute PSD psd = pd.DataFrame() for i in indexes: window_start = max(0, i - num) window_end = min(length, i + num) - slice_psd = welch( - df.iloc[window_start:window_end], - bin_width = bin_width, - scaling = scaling, - **kwargs - ) + with warnings.catch_warnings(): + if disable_warnings: + warnings.filterwarnings('ignore', '.*greater than.*', ) + slice_psd = welch( + df.iloc[window_start:window_end], + bin_width=bin_width, + scaling=scaling, + **kwargs + ) + if octave_bins is not None: - slice_psd = to_octave(slice_psd, octave_bins=octave_bins, fstart = fstart, agg=agg) + with warnings.catch_warnings(): + if disable_warnings: + warnings.filterwarnings('ignore', '.*empty frequency.*') + slice_psd = to_octave(slice_psd, octave_bins=octave_bins, fstart=fstart, agg=agg) if add_resultant: slice_psd['Resultant'] = slice_psd.sum(axis=1) slice_psd = slice_psd.reset_index().melt(id_vars=slice_psd.index.name) slice_psd['timestamp'] = df.index[i] - psd = pd.concat([psd,slice_psd]) + psd = pd.concat([psd, slice_psd]) return psd diff --git a/endaq/calc/shock.py b/endaq/calc/shock.py index c5463a7..3409ac2 100644 --- a/endaq/calc/shock.py +++ b/endaq/calc/shock.py @@ -537,7 +537,7 @@ def rolling_shock_spectrum( - `'pvss'` specifies the Pseudo-Velocity Shock Spectrum (PVSS) :param add_resultant: if `True` (default) the column-wise resultant will also be computed - :param aggregate_axes: whether to calculate the column-wise resultant (`True`) + :param add_resultant: whether to calculate the column-wise resultant (`True`) or calculate spectra along each column independently (`False`; default) :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 @@ -562,53 +562,46 @@ def rolling_shock_spectrum( Surface plots, and Animations """ - if disable_warnings: - warnings.filterwarnings('ignore', '.*too short*', ) - - length = len(df) - - #Define center index locations of each slice if not provided - if indexes is None: - if index_values is not None: - indexes = np.zeros(len(index_values),int) - for i in range(len(indexes)): - indexes[i] = int((np.abs(df.index - index_values[i])).argmin()) - else: - indexes = np.linspace(0, length, num_slices, endpoint = False, dtype=int) - indexes = indexes + int(indexes[1]/2) - - #Calculate slice step size - spacing = utils.sample_spacing(df) - if slice_width is None: - slice_width = spacing * length / len(indexes) - num = int(slice_width / spacing / 2) - - #Loop through and compute shock spectrum + + indexes, slice_width, num, length = utils._rolling_slice_definitions( + df, + indexes=indexes, + index_values=index_values, + num_slices=num_slices, + slice_width=slice_width + ) + + # Loop through and compute shock spectrum srs = pd.DataFrame() for i in indexes: window_start = max(0, i - num) window_end = min(length, i + num) - slice_srs = shock_spectrum( - df.iloc[window_start:window_end] * multiplier, - mode = mode, - damp = damp, - init_freq = init_freq, - bins_per_octave = bins_per_octave, - ) + with warnings.catch_warnings(): + if disable_warnings: + warnings.filterwarnings('ignore', '.*too short*', ) + slice_srs = shock_spectrum( + df.iloc[window_start:window_end] * multiplier, + mode=mode, + damp=damp, + init_freq=init_freq, + bins_per_octave=bins_per_octave, + ) if add_resultant: - slice_srs['Resultant'] = shock_spectrum( - df.iloc[window_start:window_end] * multiplier, - mode = mode, - damp = damp, - init_freq = init_freq, - bins_per_octave = bins_per_octave, - aggregate_axes=True - )['resultant'] - + with warnings.catch_warnings(): + if disable_warnings: + warnings.filterwarnings('ignore', '.*too short*', ) + slice_srs['Resultant'] = shock_spectrum( + df.iloc[window_start:window_end] * multiplier, + mode=mode, + damp=damp, + init_freq=init_freq, + bins_per_octave=bins_per_octave, + aggregate_axes=True + )['resultant'] slice_srs = slice_srs.reset_index().melt(id_vars=slice_srs.index.name) slice_srs['timestamp'] = df.index[i] - srs = pd.concat([srs,slice_srs]) + srs = pd.concat([srs, slice_srs]) return srs diff --git a/endaq/calc/utils.py b/endaq/calc/utils.py index 87e3ebc..a3bffa0 100644 --- a/endaq/calc/utils.py +++ b/endaq/calc/utils.py @@ -160,3 +160,49 @@ def resample(df: pd.DataFrame, sample_rate: Optional[float] = None) -> pd.DataFr resampled_df.index.name = df.index.name return resampled_df + + +def _rolling_slice_definitions( + df: pd.DataFrame, + num_slices: int = 100, + indexes: np.array = None, + index_values: np.array = None, + slice_width: float = None, +): + """ + Compute parameters needed to define index locations and slice width for rolling computations + + :param df: the input dataframe with an index defining the time in seconds or datetime + :param num_slices: the number of slices to split the time series into, default is 100, + this is ignored if `indexes` is defined + :param indexes: the center index locations (not value) of each slice to compute the FFT + :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 + :return: a tuple of `indexes`, `slice_width`, `num`, and `length` + + See example use cases and syntax at :py:func:`~endaq.plot.spectrum_over_time()` + which visualizes the output of this function in Heatmaps, Waterfall plots, + Surface plots, and Animations + + """ + + length = len(df) + + # Define center index locations of each slice if not provided + if indexes is None: + if index_values is not None: + indexes = np.zeros(len(index_values), int) + for i in range(len(indexes)): + indexes[i] = int((np.abs(df.index - index_values[i])).argmin()) + else: + indexes = np.linspace(0, length, num_slices, endpoint=False, dtype=int) + indexes = indexes + int(indexes[1] / 2) + + # Calculate slice step size + spacing = sample_spacing(df) + if slice_width is None: + slice_width = spacing * length / len(indexes) + num = int(slice_width / spacing / 2) + + return indexes, slice_width, num, length From 450070c035007c797d06897928f0fa36e958795c Mon Sep 17 00:00:00 2001 From: "MIDE\\shanly" Date: Sun, 10 Apr 2022 15:58:32 -0400 Subject: [PATCH 19/24] reformated the plots per pycharm, changed variable names, and added some more plot options --- endaq/calc/psd.py | 6 +- endaq/calc/shock.py | 5 +- endaq/plot/plots.py | 343 ++++++++++++++++++++++++++++---------------- 3 files changed, 228 insertions(+), 126 deletions(-) diff --git a/endaq/calc/psd.py b/endaq/calc/psd.py index ea532d6..bd65469 100644 --- a/endaq/calc/psd.py +++ b/endaq/calc/psd.py @@ -402,7 +402,7 @@ def rolling_psd( bin_width = 1 / slice_width # Loop through and compute PSD - psd = pd.DataFrame() + psd_df = pd.DataFrame() for i in indexes: window_start = max(0, i - num) window_end = min(length, i + num) @@ -427,6 +427,6 @@ def rolling_psd( slice_psd = slice_psd.reset_index().melt(id_vars=slice_psd.index.name) slice_psd['timestamp'] = df.index[i] - psd = pd.concat([psd, slice_psd]) + psd_df = pd.concat([psd_df, slice_psd]) - return psd + return psd_df diff --git a/endaq/calc/shock.py b/endaq/calc/shock.py index 3409ac2..9142444 100644 --- a/endaq/calc/shock.py +++ b/endaq/calc/shock.py @@ -532,9 +532,8 @@ def rolling_shock_spectrum( :param damp: the damping coefficient `ζ`, related to the Q-factor by `ζ = 1/(2Q)`; defaults to 0 :param mode: the type of spectrum to calculate: - - - `'srs'` (default) specifies the Shock Response Spectrum (SRS) - - `'pvss'` specifies the Pseudo-Velocity Shock Spectrum (PVSS) + * `'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`) diff --git a/endaq/plot/plots.py b/endaq/plot/plots.py index 3dc1e26..d067386 100644 --- a/endaq/plot/plots.py +++ b/endaq/plot/plots.py @@ -1,5 +1,8 @@ from __future__ import annotations +import typing +import warnings + import numpy as np import pandas as pd import plotly.graph_objects as go @@ -276,7 +279,7 @@ def gen_map(df_map: pd.DataFrame, mapbox_access_token: str, filter_points_by_pos def octave_spectrogram(df: pd.DataFrame, window: float, bins_per_octave: int = 3, freq_start: float = 20.0, max_freq: float = float('inf'), db_scale: bool = True, log_scale_y_axis: bool = True - ) -> go.Figure: + ) -> tuple[pd.DataFrame, go.Figure]: """ Produces an octave spectrogram of the given data, this is a wrapper around :py:func:`~endaq.calc.psd.rolling_psd()` and :py:func:`~endaq.plot.spectrum_over_time()` @@ -294,13 +297,13 @@ def octave_spectrogram(df: pd.DataFrame, window: float, bins_per_octave: int = 3 """ if len(df.columns) != 1: raise ValueError("The parameter 'df' must have only one column of data!") - - num_slices = int(len(df)*sample_spacing(df) / window) - df_psd = rolling_psd(df, num_slices=num_slices, octave_bins = bins_per_octave, - fstart = freq_start) - - return df_psd, spectrum_over_time(df_psd, y_limit=max_freq, log_z = db_scale, log_y = log_scale_y_axis) + num_slices = int(len(df) * sample_spacing(df) / window) + + df_psd = rolling_psd(df, num_slices=num_slices, octave_bins=bins_per_octave, + fstart=freq_start) + + return df_psd, spectrum_over_time(df_psd, freq_max=max_freq, log_val=db_scale, log_freq=log_scale_y_axis) def octave_psd_bar_plot(df: pd.DataFrame, bins_per_octave: int = 3, f_start: float = 20.0, yaxis_title: str = '', @@ -509,46 +512,63 @@ def animate_quaternion(df, rate=6., scale=1.): def spectrum_over_time( df: pd.DataFrame, - plot_type: typing.Literal["Heatmap", "Surface", "Waterfall", "Animation"] = "Heatmap", + plot_type: typing.Literal["Heatmap", + "Surface", + "Waterfall", + "Animation", + "Peak", + "Lines"] = "Heatmap", var_column: str = 'variable', var_to_process: str = None, - x_column: str = 'timestamp', - y_column: str = 'frequency (Hz)', - z_column: str = 'value', - y_limit: float = None, - log_y: bool = False, - log_z: bool = False, + time_column: str = 'timestamp', + freq_column: str = 'frequency (Hz)', + val_column: str = 'value', + freq_min: float = None, + freq_max: float = None, + log_freq: bool = False, + log_val: bool = False, round_time: bool = True, waterfall_line_sequence: bool = True, waterfall_line_color: str = '#EE7F27', + min_median_max_line_color: str = '#6914F0', ) -> go.Figure: """ Generate a 3D Plotly figure from a stacked spectrum to visualize how the frequency content changes over time :param df: the input dataframe with columns defining the frequency content, timestamps, values, and variables, see the following functions which provides outputs that would then be passed into this function as an input: - - :py:func:`~endaq.calc.fft.rolling_fft()` - - :py:func:`~endaq.calc.psd.rolling_psd()` - - :py:func:`~endaq.calc.shock.rolling_shock_spectrum()` + * :py:func:`~endaq.calc.fft.rolling_fft()` + * :py:func:`~endaq.calc.psd.rolling_psd()` + * :py:func:`~endaq.calc.shock.rolling_shock_spectrum()` + * :py:func:`~endaq.batch.GetDataBuilder.add_psd()` + * :py:func:`~endaq.batch.GetDataBuilder.add_pvss()` :param plot_type: the type of plot to display the spectrum, options are: - - `Heatmap`: a 2D visualization with the color defining the z value (the default) - - `Surface`: similar to `Heatmap` but the z value is also projected "off the page" - - `Waterfall`: distinct lines are plotted per time slice in a 3D view - - `Animation`: a 2D display of the waterfall but the spectrum changes with animation frames + * `Heatmap`: a 2D visualization with the color defining the z value (the default) + * `Surface`: similar to `Heatmap` but the z value is also projected "off the page" + * `Waterfall`: distinct lines are plotted per time slice in a 3D view + * `Animation`: a 2D display of the waterfall but the spectrum changes with animation frames + * `Peak`: per timestamp the peak frequency is determined and plotted against time + * `Lines`: the value in each frequency bin is plotted against time :param var_column: the column name in the dataframe that defines the different variables, default is `"variable"` :param var_to_process: the variable value in the `var_column` to filter the input df down to, if none is provided (the default) this function will filter to the first value - :param x_column: the column name in the dataframe that defines the timestamps, default is `"timestamp"` - :param y_column: the column name in the dataframe that defines the frequency, default is `"frequency (Hz)"` - :param z_column: the column name in the dataframe that defines the values, default is `"value"` - :param y_limit: the limit of the y axis (frequency) to include in the figure, + :param time_column: the column name in the dataframe that defines the timestamps, default is `"timestamp"` + :param freq_column: the column name in the dataframe that defines the frequency, default is `"frequency (Hz)"` + :param val_column: the column name in the dataframe that defines the values, default is `"value"` + :param freq_min: the minimum of the y axis (frequency) to include in the figure, + default is None meaning it will display all content + :param freq_max: the maximum of the y axis (frequency) to include in the figure, default is None meaning it will display all content - :param log_y: if `True` the y axis (frequency) will be in a log scale, default is `False` - :param log_z: if `True` the z axis (values) will be in a log scale, default is `False` - :param round_time: if `True` (default) the time values will be rounded to the nearest second for datetimes and hundredths of a second for floats + :param log_freq: if `True` the frequency will be in a log scale, default is `False` + :param log_val: if `True` the values will be in a log scale, default is `False` + :param round_time: if `True` (default) the time values will be rounded to the nearest second for datetimes and + hundredths of a second for floats :param waterfall_line_sequence: if `True` the waterfall line colors are defined with a color scale, if `False` all lines will have the same color, default is `True` - :param waterfall_line_color: the color to use for all lines in the Waterfall plot if `waterfall_line_sequence` is `False` + :param waterfall_line_color: the color to use for all lines in the Waterfall plot if `waterfall_line_sequence` is + `False` + :param min_median_max_line_color: the color to use for the min, max, and median lines in the Animation, if set to + `None` these lines won't be added, default is `'#6914F0'` :return: a Plotly figure visualizing the spectrum over time @@ -556,6 +576,7 @@ def spectrum_over_time( as it revved the engine which resulted in changing frequency content .. code:: python + import endaq endaq.plot.utilities.set_theme() import pandas as pd @@ -568,23 +589,35 @@ def spectrum_over_time( fft = endaq.calc.fft.rolling_fft(df_vibe, num_slices=200, add_resultant=True) #Visualize the Rolling FFT as a Heatmap - heatmap = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Heatmap', log_y=False, log_z=False, y_limit=200, var_to_process='Resultant') + heatmap = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Heatmap', freq_max=200, var_to_process='Resultant') heatmap.show() + #Visualize as the peak frequency + peak = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Peak', freq_max=200, var_to_process='Resultant') + peak.show() + #Visualize as a Surface Plot - surface = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Surface', log_y=False, log_z=False, y_limit=200, var_to_process='Resultant') + surface = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Surface', freq_max=200, var_to_process='Resultant') surface.show() #Visualize as a Waterfall - waterfall = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Waterfall', log_y=False, log_z=False, y_limit=200, var_to_process='Resultant') + waterfall = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Waterfall', freq_max=200, var_to_process='Resultant') waterfall.show() #Visualize as an Animation - animation = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Animation', log_y=False, log_z=False, y_limit=200, var_to_process='Resultant') + animation = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Animation', log_val=True, freq_max=200, var_to_process='Resultant') animation.show() + #Calculate a Rolling PSD with Units as RMS**2 + psd = endaq.calc.psd.rolling_psd(df_vibe, num_slices=200, scaling='parseval', add_resultant=True, octave_bins=1, fstart=1, agg='sum') + psd['value'] = psd['value'] ** 0.5 + + #Visualize the energy in each frequency bin + lines = endaq.plot.plots.spectrum_over_time(psd, plot_type = 'Lines', log_val=False, var_to_process='Resultant') + lines.show() + .. plotly:: - :fig-vars: heatmap, surface, waterfall, animation + :fig-vars: heatmap, peak, surface, waterfall, animation, lines import endaq endaq.plot.utilities.set_theme() @@ -598,146 +631,216 @@ def spectrum_over_time( fft = endaq.calc.fft.rolling_fft(df_vibe, num_slices=200, add_resultant=True) #Visualize the Rolling FFT as a Heatmap - heatmap = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Heatmap', log_y=False, log_z=False, y_limit=200, var_to_process='Resultant') + heatmap = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Heatmap', freq_max=200, var_to_process='Resultant') heatmap.show() + #Visualize as the peak frequency + peak = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Peak', freq_max=200, var_to_process='Resultant') + peak.show() + #Visualize as a Surface Plot - surface = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Surface', log_y=False, log_z=False, y_limit=200, var_to_process='Resultant') + surface = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Surface', freq_max=200, var_to_process='Resultant') surface.show() #Visualize as a Waterfall - waterfall = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Waterfall', log_y=False, log_z=False, y_limit=200, var_to_process='Resultant') + waterfall = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Waterfall', freq_max=200, var_to_process='Resultant') waterfall.show() #Visualize as an Animation - animation = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Animation', log_y=False, log_z=False, y_limit=200, var_to_process='Resultant') + animation = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Animation', log_val=True, freq_max=200, var_to_process='Resultant') animation.show() + #Calculate a Rolling PSD with Units as RMS**2 + psd = endaq.calc.psd.rolling_psd(df_vibe, num_slices=200, scaling='parseval', add_resultant=True, octave_bins=1, fstart=1, agg='sum') + psd['value'] = psd['value'] ** 0.5 + + #Visualize the energy in each frequency bin + lines = endaq.plot.plots.spectrum_over_time(psd, plot_type = 'Lines', log_val=False, var_to_process='Resultant') + lines.show() + """ df = df.copy() - #Filter to one variable + # Filter to one variable if var_to_process is None: var_to_process = df[var_column].unique()[0] df = df.loc[df[var_column] == var_to_process] - #Round time + # Round time if round_time: - if isinstance(df[x_column].iloc[0], datetime.datetime): - df[x_column] = df[x_column].round('s') + if isinstance(df[time_column].iloc[0], datetime.datetime): + df[time_column] = df[time_column].round('s') else: - df[x_column] = np.round(df[x_column].to_numpy(),2) + df[time_column] = np.round(df[time_column].to_numpy(), 2) - #Filter frequency - if y_limit is not None: - df = df.loc[df[y_column] < y_limit] + # Filter frequency + df = df.loc[df[freq_column] > 0.0] + if freq_max is not None: + df = df.loc[df[freq_column] < freq_max] + if freq_min is not None: + df = df.loc[df[freq_column] > freq_min] - #Remove 0s - df = df.loc[df[z_column] > 0] + # Remove 0s + df = df.loc[df[val_column] > 0] - #Check Length of Dataframe + # Check Length of Dataframe if len(df) > 50000: warnings.warn( - "plot data is very large, may be unresponsive", + "plot data is very large, may be unresponsive, suggest limiting frequency range and/or using less slices", RuntimeWarning, ) - #Heatmap & Surface - if plot_type in ["Heatmap", "Surface"]: - df[y_column] = np.round(df[y_column].to_numpy(),2) - df = df.pivot(columns=x_column,index=y_column,values=z_column) + # Create pivot table + df_pivot = df.copy() + first_step = df_pivot[freq_column].iloc[1] - df_pivot[freq_column].iloc[0] + second_step = df_pivot[freq_column].iloc[2] - df_pivot[freq_column].iloc[1] + if first_step == second_step: + round_freq = np.round(df_pivot[freq_column].min(), 0) + df_pivot[freq_column] = np.round(df_pivot[freq_column].to_numpy() / round_freq, 0) * round_freq + df_pivot = df_pivot.pivot(columns=time_column, index=freq_column, values=val_column) - #Deal with Dateimes + # Heatmap & Surface + if plot_type in ["Heatmap", "Surface"]: + # Deal with Datetime x_type = float - if isinstance(df.columns[0], datetime.datetime): + if isinstance(df_pivot.columns[0], datetime.datetime): if plot_type == "Surface": - df.columns = (df.columns - df.columns[0]).total_seconds() + df_pivot.columns = (df_pivot.columns - df_pivot.columns[0]).total_seconds() else: x_type = str - #Build Dictionary of Plot Data, Apply Log Scale If Needed + # Build Dictionary of Plot Data, Apply Log Scale If Needed data_dict = { - 'x': df.columns.astype(x_type), - 'y': df.index.astype(float), - 'z': df.to_numpy().astype(float), - 'connectgaps' : True + 'x': df_pivot.columns.astype(x_type), + 'y': df_pivot.index.astype(float), + 'z': df_pivot.to_numpy().astype(float), + 'connectgaps': True } - if log_z: + if log_val: data_dict['z'] = np.log10(data_dict['z']) - #Generate Figures + # Generate Figures if plot_type == "Heatmap": - fig = go.Figure(data = go.Heatmap(data_dict, zsmooth='best')).update_layout( + fig = go.Figure(data=go.Heatmap(data_dict, zsmooth='best')).update_layout( xaxis_title_text='Timestamp', yaxis_title_text='Frequency (Hz)') - if log_y: - fig.update_layout(yaxis_type='log') + if log_freq: + fig.update_layout(yaxis_type='log') else: - fig = go.Figure(data = go.Surface(data_dict)) + fig = go.Figure(data=go.Surface(data_dict)) fig.update_traces(showscale=False) - #Waterfall + # Waterfall elif plot_type == 'Waterfall': - #Define Colors + # Define Colors if waterfall_line_sequence: color_sequence = colors.sample_colorscale( - [[0.0, '#6914F0'], - [0.2, '#3764FF'], - [0.4, '#2DB473'], - [0.6, '#FAC85F'], - [0.8, '#EE7F27'], - [1.0, '#D72D2D']], len(df[x_column].unique())) + [[0.0, '#6914F0'], + [0.2, '#3764FF'], + [0.4, '#2DB473'], + [0.6, '#FAC85F'], + [0.8, '#EE7F27'], + [1.0, '#D72D2D']], len(df[time_column].unique())) else: color_sequence = [waterfall_line_color] - #Deal with Datetime - df['label_column'] = df[x_column] - if isinstance(df[x_column].iloc[0], datetime.datetime): - df[x_column] = (df[x_column] - df[x_column].iloc[0]).dt.total_seconds() - df['label_column'] = df['label_column'].dt.tz_convert(None).astype(str) + # Deal with Datetime + df['label_column'] = df[time_column] + if isinstance(df[time_column].iloc[0], datetime.datetime): + df[time_column] = (df[time_column] - df[time_column].iloc[0]).dt.total_seconds() + df['label_column'] = df['label_column'].dt.tz_convert(None).astype(str) - #Generate Figure + # Generate Figure fig = px.line_3d( df, - x = x_column, - y = y_column, - z = z_column, - color = 'label_column', - color_discrete_sequence = color_sequence).update_layout( - legend_orientation = 'v', - legend_y = 1, - legend_title_text = 'Timestamps' - ) - - #Animation - elif plot_type == 'Animation': - #Define Range (If Undefined it Starts with Fits to First Frame) - range = [df['value'].min(), df['value'].max()] - if log_z: - range = np.log10(range) + x=time_column, + y=freq_column, + z=val_column, + color='label_column', + color_discrete_sequence=color_sequence).update_layout( + legend_orientation='v', + legend_y=1, + legend_title_text='Timestamps' + ) - #Deal with Datetime - if isinstance(df[x_column].iloc[0], datetime.datetime): - df[x_column] = df[x_column].dt.tz_convert(None).astype(str) + # Animation + elif plot_type == 'Animation': + # Deal with Datetime + if isinstance(df[time_column].iloc[0], datetime.datetime): + df[time_column] = df[time_column].dt.tz_convert(None).astype(str) - #Gerate Figure + # Generate Figure fig = px.line( df, - animation_frame = x_column, - x = y_column, - y = z_column, - log_y = log_z, - log_x = log_y + animation_frame=time_column, + x=freq_column, + y=val_column, + log_y=log_val, + log_x=log_freq + ).update_layout( + showlegend=False, + yaxis_title_text='', + xaxis_title_text='Frequency (Hz)' + ) + + # Add min, max, median lines + if min_median_max_line_color is not None: + # Add Max + df_t = df_pivot.max(axis=1).dropna() + fig.add_trace(go.Scattergl( + x=df_t.index, + y=df_t, + mode='lines', + line_color=min_median_max_line_color, + )) + + # Add Min + df_t = df_pivot.min(axis=1).dropna() + fig.add_trace(go.Scattergl( + x=df_t.index, + y=df_t, + mode='lines', + line_color=min_median_max_line_color, + )) + + # Add Median + df_t = df_pivot.median(axis=1).dropna() + fig.add_trace(go.Scattergl( + x=df_t.index, + y=df_t, + mode='lines', + line_color=min_median_max_line_color, + line_dash='dash', + )) + + # Peak Frequencies + elif plot_type == 'Peak': + fig = px.line( + df_pivot.idxmax(), + log_y=log_freq ).update_layout( - showlegend = False, - yaxis_range = range, - yaxis_title_text = '', - xaxis_title_text = 'Frequency (Hz)' + showlegend=False, + yaxis_title_text='Peak Frequency (Hz)', + xaxis_title_text='Timestamp' + ) + + # Lines per Frequency Bin + elif plot_type == 'Lines': + df_pivot.index = np.round(df_pivot.index, 2) + fig = px.line( + df_pivot.T, + log_y=log_val + ).update_layout( + legend_title_text='Frequency Bin (Hz)', + yaxis_title_text='', + xaxis_title_text='Timestamp', + legend_orientation='v', + legend_y=1, ) else: raise ValueError(f"invalid plot type {plot_type}") - - #Add Labels to 3D plots + + # Add Labels to 3D plots if plot_type in ["Surface", "Waterfall"]: fig.update_scenes( aspectratio_x=2.0, @@ -745,12 +848,12 @@ def spectrum_over_time( aspectratio_z=0.3, xaxis_title_text='Timestamp', yaxis_title_text='Frequency (Hz)', - zaxis_title_text= '' - ) - if log_y: - fig.update_scenes(yaxis_type = 'log') - if log_z: - fig.update_scenes(zaxis_type = 'log') - fig.update_layout(scene_camera = dict(eye=dict(x=-1.5, y=-1.5, z=1))) + zaxis_title_text='' + ) + if log_freq: + fig.update_scenes(yaxis_type='log') + if log_val: + fig.update_scenes(zaxis_type='log') + fig.update_layout(scene_camera=dict(eye=dict(x=-1.5, y=-1.5, z=1))) return fig From 5a4ccee55b3e36b5f49166cdcadb9c38b6f7c742 Mon Sep 17 00:00:00 2001 From: "MIDE\\shanly" Date: Sun, 10 Apr 2022 23:41:14 -0400 Subject: [PATCH 20/24] rolling RMS in frequency ranges, improved docs example, changed plot type for peak frequency --- endaq/calc/psd.py | 27 +++++-- endaq/calc/shock.py | 8 +-- endaq/plot/plots.py | 171 +++++++++++++++++++++++++++++--------------- 3 files changed, 137 insertions(+), 69 deletions(-) diff --git a/endaq/calc/psd.py b/endaq/calc/psd.py index bd65469..cf05a7f 100644 --- a/endaq/calc/psd.py +++ b/endaq/calc/psd.py @@ -342,11 +342,12 @@ def rolling_psd( bin_width: float = 1.0, octave_bins: float = None, fstart: float = 1.0, - scaling: typing.Literal[None, "density", "spectrum", "parseval", "unit"] = None, + scaling: typing.Literal[None, "density", "spectrum", "parseval", "unit", "rms"] = None, agg: typing.Union[ typing.Literal["mean", "sum"], typing.Callable[[np.ndarray, int], float], ] = "mean", + freq_splits: np.array = None, num_slices: int = 100, indexes: np.array = None, index_values: np.array = None, @@ -368,9 +369,11 @@ def rolling_psd( will maintain the "energy" between the input & output, s.t. ``welch(df, scaling="parseval").sum(axis="rows")`` is roughly equal to ``df.abs().pow(2).sum(axis="rows")``; - `"unit"` will maintain the units and scale of the input data. - :param agg: the method for aggregating values into bins (only used if converting to octave); `'mean'` preserves - the PSD's area-under-the-curve, `'sum'` preserves the PSD's "energy" + `"unit"` will maintain the units and scale of the input data. `"rms"` will use `"parseval"` for the PSD + calculations and set `agg` to "sum", but then take the square root at the end + :param agg: the method for aggregating values into bins (only used if converting to octave or jagged); `'mean'` preserves + the PSD's area-under-the-curve, `'sum'` preserves the PSD's "energy" + :param freq_splits: the boundaries of the frequency bins to pass to :py:func:`~endaq.calc.psd.to_jagged()` :param num_slices: the number of slices to split the time series into, default is 100, this is ignored if `indexes` is defined :param indexes: the center index locations (not value) of each slice to compute the PSD @@ -397,6 +400,13 @@ def rolling_psd( slice_width=slice_width ) + # Allow for RMS calculations + exp = 1 + if scaling == "rms": + scaling = "parseval" + exp = 0.5 + agg = "sum" + # Define bin_width if octave spacing if octave_bins is not None: bin_width = 1 / slice_width @@ -416,16 +426,19 @@ def rolling_psd( **kwargs ) - if octave_bins is not None: + if octave_bins is not None or freq_splits is not None: with warnings.catch_warnings(): if disable_warnings: warnings.filterwarnings('ignore', '.*empty frequency.*') - slice_psd = to_octave(slice_psd, octave_bins=octave_bins, fstart=fstart, agg=agg) + if freq_splits is None: + slice_psd = to_octave(slice_psd, octave_bins=octave_bins, fstart=fstart, agg=agg) + else: + slice_psd = to_jagged(slice_psd, freq_splits=freq_splits, agg=agg) if add_resultant: slice_psd['Resultant'] = slice_psd.sum(axis=1) - slice_psd = slice_psd.reset_index().melt(id_vars=slice_psd.index.name) + slice_psd = (slice_psd ** exp).reset_index().melt(id_vars=slice_psd.index.name) slice_psd['timestamp'] = df.index[i] psd_df = pd.concat([psd_df, slice_psd]) diff --git a/endaq/calc/shock.py b/endaq/calc/shock.py index 9142444..8c92ac2 100644 --- a/endaq/calc/shock.py +++ b/endaq/calc/shock.py @@ -408,7 +408,7 @@ def shock_spectrum( freqs: np.ndarray = None, init_freq: float = 0.5, bins_per_octave: float = 12.0, - damp: float = 0.0, + damp: float = 0.05, mode: typing.Literal["srs", "pvss"] = "srs", two_sided: bool = False, aggregate_axes: bool = False, @@ -423,7 +423,7 @@ def shock_spectrum( 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 :param damp: the damping coefficient `ζ`, related to the Q-factor by - `ζ = 1/(2Q)`; defaults to 0 + `ζ = 1/(2Q)`; defaults to 0.05 :param mode: the type of spectrum to calculate: - `'srs'` (default) specifies the Shock Response Spectrum (SRS) - `'pvss'` specifies the Pseudo-Velocity Shock Spectrum (PVSS) @@ -513,7 +513,7 @@ def shock_spectrum( def rolling_shock_spectrum( df: pd.DataFrame, - damp: float = 0.0, + damp: float = 0.05, mode: typing.Literal["srs", "pvss"] = "srs", add_resultant: bool = True, init_freq: float = 0.5, @@ -530,7 +530,7 @@ def rolling_shock_spectrum( :param df: the input dataframe with an index defining the time in seconds or datetime :param damp: the damping coefficient `ζ`, related to the Q-factor by - `ζ = 1/(2Q)`; defaults to 0 + `ζ = 1/(2Q)`; defaults to 0.05 :param mode: the type of spectrum to calculate: * `'srs'` (default) specifies the Shock Response Spectrum (SRS) * `'pvss'` specifies the Pseudo-Velocity Shock Spectrum (PVSS) diff --git a/endaq/plot/plots.py b/endaq/plot/plots.py index d067386..02d576e 100644 --- a/endaq/plot/plots.py +++ b/endaq/plot/plots.py @@ -529,6 +529,7 @@ def spectrum_over_time( log_val: bool = False, round_time: bool = True, waterfall_line_sequence: bool = True, + zsmooth: str = "best", waterfall_line_color: str = '#EE7F27', min_median_max_line_color: str = '#6914F0', ) -> go.Figure: @@ -537,18 +538,18 @@ def spectrum_over_time( :param df: the input dataframe with columns defining the frequency content, timestamps, values, and variables, see the following functions which provides outputs that would then be passed into this function as an input: - * :py:func:`~endaq.calc.fft.rolling_fft()` - * :py:func:`~endaq.calc.psd.rolling_psd()` - * :py:func:`~endaq.calc.shock.rolling_shock_spectrum()` - * :py:func:`~endaq.batch.GetDataBuilder.add_psd()` - * :py:func:`~endaq.batch.GetDataBuilder.add_pvss()` + * :py:func:`~endaq.calc.fft.rolling_fft()` + * :py:func:`~endaq.calc.psd.rolling_psd()` + * :py:func:`~endaq.calc.shock.rolling_shock_spectrum()` + * :py:func:`~endaq.batch.GetDataBuilder.add_psd()` + * :py:func:`~endaq.batch.GetDataBuilder.add_pvss()` :param plot_type: the type of plot to display the spectrum, options are: - * `Heatmap`: a 2D visualization with the color defining the z value (the default) - * `Surface`: similar to `Heatmap` but the z value is also projected "off the page" - * `Waterfall`: distinct lines are plotted per time slice in a 3D view - * `Animation`: a 2D display of the waterfall but the spectrum changes with animation frames - * `Peak`: per timestamp the peak frequency is determined and plotted against time - * `Lines`: the value in each frequency bin is plotted against time + * `Heatmap`: a 2D visualization with the color defining the z value (the default) + * `Surface`: similar to `Heatmap` but the z value is also projected "off the page" + * `Waterfall`: distinct lines are plotted per time slice in a 3D view + * `Animation`: a 2D display of the waterfall but the spectrum changes with animation frames + * `Peak`: per timestamp the peak frequency is determined and plotted against time + * `Lines`: the value in each frequency bin is plotted against time :param var_column: the column name in the dataframe that defines the different variables, default is `"variable"` :param var_to_process: the variable value in the `var_column` to filter the input df down to, if none is provided (the default) this function will filter to the first value @@ -565,6 +566,8 @@ def spectrum_over_time( hundredths of a second for floats :param waterfall_line_sequence: if `True` the waterfall line colors are defined with a color scale, if `False` all lines will have the same color, default is `True` + :param zsmooth: the Plotly smooth algorithm to use in the heatmap, default is `"best"` which looks ideal but + `"fast"` will be more responsive, or `False` will attempt no smoothing :param waterfall_line_color: the color to use for all lines in the Waterfall plot if `waterfall_line_sequence` is `False` :param min_median_max_line_color: the color to use for the min, max, and median lines in the Animation, if set to @@ -577,87 +580,139 @@ def spectrum_over_time( .. code:: python + import pandas as pd import endaq + + # Set Theme endaq.plot.utilities.set_theme() - import pandas as pd - #Get the Data - df_vibe = pd.read_csv('https://info.endaq.com/hubfs/data/motorcycle-vibration-moving-frequency.csv',index_col=0) - df_vibe = df_vibe - df_vibe.median() + # Get Vibration Data + df_vibe = pd.read_csv('https://info.endaq.com/hubfs/data/motorcycle-vibration-moving-frequency.csv', index_col=0) - #Calculate a Rolling FFT + # Calculate a Rolling FFT fft = endaq.calc.fft.rolling_fft(df_vibe, num_slices=200, add_resultant=True) - #Visualize the Rolling FFT as a Heatmap - heatmap = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Heatmap', freq_max=200, var_to_process='Resultant') + # Visualize the Rolling FFT as a Heatmap + heatmap = endaq.plot.plots.spectrum_over_time(fft, plot_type='Heatmap', freq_max=200, var_to_process='Resultant') heatmap.show() - #Visualize as the peak frequency - peak = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Peak', freq_max=200, var_to_process='Resultant') + # Plot the Peak Frequency vs Time + peak = endaq.plot.plots.spectrum_over_time(fft, plot_type='Peak', freq_max=200, var_to_process='Resultant') peak.show() - #Visualize as a Surface Plot - surface = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Surface', freq_max=200, var_to_process='Resultant') + # Visualize as a Surface Plot + surface = endaq.plot.plots.spectrum_over_time(fft, plot_type='Surface', freq_max=200, var_to_process='Resultant') surface.show() - #Visualize as a Waterfall - waterfall = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Waterfall', freq_max=200, var_to_process='Resultant') + # Visualize as a Waterfall + waterfall = endaq.plot.plots.spectrum_over_time(fft, plot_type='Waterfall', freq_max=200, var_to_process='Resultant') waterfall.show() - #Visualize as an Animation - animation = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Animation', log_val=True, freq_max=200, var_to_process='Resultant') + # Get a Longer Dataset with DatetimeIndex + engine = endaq.ide.get_primary_sensor_data('https://info.endaq.com/hubfs/data/Commute.ide', measurement_type='accel', + time_mode='datetime') + + # Compute PSD + psd = endaq.calc.psd.rolling_psd(engine, num_slices=500, add_resultant=True, octave_bins=12, fstart=4) + + # Visualize as a Heatmap + heatmap2 = endaq.plot.plots.spectrum_over_time(psd, plot_type='Heatmap', var_to_process='Resultant', zsmooth='best', + log_freq=True, log_val=True) + heatmap2.show() + + # Visualize as an Animation + animation = endaq.plot.plots.spectrum_over_time(psd, plot_type='Animation', var_to_process='Resultant', + log_freq=True, log_val=True) animation.show() - #Calculate a Rolling PSD with Units as RMS**2 - psd = endaq.calc.psd.rolling_psd(df_vibe, num_slices=200, scaling='parseval', add_resultant=True, octave_bins=1, fstart=1, agg='sum') - psd['value'] = psd['value'] ** 0.5 + # Use Rolling PSD to Calculate RMS in Certain Frequency Bins + rms = endaq.calc.psd.rolling_psd(engine, num_slices=500, scaling='rms', add_resultant=True, + freq_splits=[1, 20, 60, 300, 3000]) - #Visualize the energy in each frequency bin - lines = endaq.plot.plots.spectrum_over_time(psd, plot_type = 'Lines', log_val=False, var_to_process='Resultant') + # Plot the RMS per Frequency Bin Over Time + lines = endaq.plot.plots.spectrum_over_time(rms, plot_type='Lines', log_val=False, var_to_process='Resultant') lines.show() + # Compute Pseudo Velocity at Specific Times + pvss = endaq.calc.shock.rolling_shock_spectrum(engine, slice_width=2.0, add_resultant=True, + mode='pvss', init_freq=4, damp=0.05, + index_values=pd.DatetimeIndex(['2016-08-02 12:07:15', + '2016-08-02 12:08:01', + '2016-08-02 12:10:06'], tz='UTC')) + + # Visualize as a Waterfall + waterfall2 = endaq.plot.plots.spectrum_over_time(pvss, plot_type='Waterfall', log_freq=True, log_val=True, + var_to_process='Resultant', waterfall_line_sequence=False) + waterfall2.show() + .. plotly:: - :fig-vars: heatmap, peak, surface, waterfall, animation, lines + :fig-vars: heatmap, peak, surface, waterfall, heatmap2, animation, lines, waterfall2 + import pandas as pd import endaq + + # Set Theme endaq.plot.utilities.set_theme() - import pandas as pd - #Get the Data - df_vibe = pd.read_csv('https://info.endaq.com/hubfs/data/motorcycle-vibration-moving-frequency.csv',index_col=0) - df_vibe = df_vibe - df_vibe.median() + # Get Vibration Data + df_vibe = pd.read_csv('https://info.endaq.com/hubfs/data/motorcycle-vibration-moving-frequency.csv', index_col=0) - #Calculate a Rolling FFT + # Calculate a Rolling FFT fft = endaq.calc.fft.rolling_fft(df_vibe, num_slices=200, add_resultant=True) - #Visualize the Rolling FFT as a Heatmap - heatmap = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Heatmap', freq_max=200, var_to_process='Resultant') + # Visualize the Rolling FFT as a Heatmap + heatmap = endaq.plot.plots.spectrum_over_time(fft, plot_type='Heatmap', freq_max=200, var_to_process='Resultant') heatmap.show() - #Visualize as the peak frequency - peak = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Peak', freq_max=200, var_to_process='Resultant') + # Plot the Peak Frequency vs Time + peak = endaq.plot.plots.spectrum_over_time(fft, plot_type='Peak', freq_max=200, var_to_process='Resultant') peak.show() - #Visualize as a Surface Plot - surface = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Surface', freq_max=200, var_to_process='Resultant') + # Visualize as a Surface Plot + surface = endaq.plot.plots.spectrum_over_time(fft, plot_type='Surface', freq_max=200, var_to_process='Resultant') surface.show() - #Visualize as a Waterfall - waterfall = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Waterfall', freq_max=200, var_to_process='Resultant') + # Visualize as a Waterfall + waterfall = endaq.plot.plots.spectrum_over_time(fft, plot_type='Waterfall', freq_max=200, var_to_process='Resultant') waterfall.show() - #Visualize as an Animation - animation = endaq.plot.plots.spectrum_over_time(fft, plot_type = 'Animation', log_val=True, freq_max=200, var_to_process='Resultant') + # Get a Longer Dataset with DatetimeIndex + engine = endaq.ide.get_primary_sensor_data('https://info.endaq.com/hubfs/data/Commute.ide', measurement_type='accel', + time_mode='datetime') + + # Compute PSD + psd = endaq.calc.psd.rolling_psd(engine, num_slices=500, add_resultant=True, octave_bins=12, fstart=4) + + # Visualize as a Heatmap + heatmap2 = endaq.plot.plots.spectrum_over_time(psd, plot_type='Heatmap', var_to_process='Resultant', zsmooth='best', + log_freq=True, log_val=True) + heatmap2.show() + + # Visualize as an Animation + animation = endaq.plot.plots.spectrum_over_time(psd, plot_type='Animation', var_to_process='Resultant', + log_freq=True, log_val=True) animation.show() - #Calculate a Rolling PSD with Units as RMS**2 - psd = endaq.calc.psd.rolling_psd(df_vibe, num_slices=200, scaling='parseval', add_resultant=True, octave_bins=1, fstart=1, agg='sum') - psd['value'] = psd['value'] ** 0.5 + # Use Rolling PSD to Calculate RMS in Certain Frequency Bins + rms = endaq.calc.psd.rolling_psd(engine, num_slices=500, scaling='rms', add_resultant=True, + freq_splits=[1, 20, 60, 300, 3000]) - #Visualize the energy in each frequency bin - lines = endaq.plot.plots.spectrum_over_time(psd, plot_type = 'Lines', log_val=False, var_to_process='Resultant') + # Plot the RMS per Frequency Bin Over Time + lines = endaq.plot.plots.spectrum_over_time(rms, plot_type='Lines', log_val=False, var_to_process='Resultant') lines.show() + # Compute Pseudo Velocity at Specific Times + pvss = endaq.calc.shock.rolling_shock_spectrum(engine, slice_width=2.0, add_resultant=True, + mode='pvss', init_freq=4, damp=0.05, + index_values=pd.DatetimeIndex(['2016-08-02 12:07:15', + '2016-08-02 12:08:01', + '2016-08-02 12:10:06'], tz='UTC')) + + # Visualize as a Waterfall + waterfall2 = endaq.plot.plots.spectrum_over_time(pvss, plot_type='Waterfall', log_freq=True, log_val=True, + var_to_process='Resultant', waterfall_line_sequence=False) + waterfall2.show() + """ df = df.copy() @@ -684,7 +739,7 @@ def spectrum_over_time( df = df.loc[df[val_column] > 0] # Check Length of Dataframe - if len(df) > 50000: + if len(df) > 100000: warnings.warn( "plot data is very large, may be unresponsive, suggest limiting frequency range and/or using less slices", RuntimeWarning, @@ -721,7 +776,7 @@ def spectrum_over_time( # Generate Figures if plot_type == "Heatmap": - fig = go.Figure(data=go.Heatmap(data_dict, zsmooth='best')).update_layout( + fig = go.Figure(data=go.Heatmap(data_dict, zsmooth=zsmooth)).update_layout( xaxis_title_text='Timestamp', yaxis_title_text='Frequency (Hz)') if log_freq: fig.update_layout(yaxis_type='log') @@ -786,7 +841,7 @@ def spectrum_over_time( if min_median_max_line_color is not None: # Add Max df_t = df_pivot.max(axis=1).dropna() - fig.add_trace(go.Scattergl( + fig.add_trace(go.Scatter( x=df_t.index, y=df_t, mode='lines', @@ -795,7 +850,7 @@ def spectrum_over_time( # Add Min df_t = df_pivot.min(axis=1).dropna() - fig.add_trace(go.Scattergl( + fig.add_trace(go.Scatter( x=df_t.index, y=df_t, mode='lines', @@ -804,7 +859,7 @@ def spectrum_over_time( # Add Median df_t = df_pivot.median(axis=1).dropna() - fig.add_trace(go.Scattergl( + fig.add_trace(go.Scatter( x=df_t.index, y=df_t, mode='lines', @@ -814,7 +869,7 @@ def spectrum_over_time( # Peak Frequencies elif plot_type == 'Peak': - fig = px.line( + fig = px.scatter( df_pivot.idxmax(), log_y=log_freq ).update_layout( From 979fb76aa84f302eaaa759dfcfd0095b1403ca79 Mon Sep 17 00:00:00 2001 From: "MIDE\\shanly" Date: Tue, 12 Apr 2022 16:00:32 -0400 Subject: [PATCH 21/24] rolling spectrum unit tests --- endaq/calc/shock.py | 10 ++-- endaq/plot/plots.py | 75 +++++++++++++++++---------- tests/calc/test_fft.py | 106 ++++++++++++++++++++++++++------------- tests/calc/test_utils.py | 32 ++++++++++++ tests/plot/test_plots.py | 36 +++++++++++++ 5 files changed, 189 insertions(+), 70 deletions(-) create mode 100644 tests/plot/test_plots.py diff --git a/endaq/calc/shock.py b/endaq/calc/shock.py index 8c92ac2..0bede15 100644 --- a/endaq/calc/shock.py +++ b/endaq/calc/shock.py @@ -522,7 +522,6 @@ def rolling_shock_spectrum( indexes: np.array = None, index_values: np.array = None, slice_width: float = None, - multiplier: float = 1.0, disable_warnings: bool = True, ) -> pd.DataFrame: """ @@ -532,6 +531,7 @@ def rolling_shock_spectrum( :param damp: the damping coefficient `ζ`, related to the Q-factor by `ζ = 1/(2Q)`; defaults to 0.05 :param mode: the type of spectrum to calculate: + * `'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 @@ -549,10 +549,6 @@ def rolling_shock_spectrum( 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 multiplier: applies a scale to the output - - - 386.09 to convert from g to inches (in) - - 9806.65 to convert from g to millimeters (mm) :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 @@ -579,7 +575,7 @@ def rolling_shock_spectrum( if disable_warnings: warnings.filterwarnings('ignore', '.*too short*', ) slice_srs = shock_spectrum( - df.iloc[window_start:window_end] * multiplier, + df.iloc[window_start:window_end], mode=mode, damp=damp, init_freq=init_freq, @@ -590,7 +586,7 @@ def rolling_shock_spectrum( if disable_warnings: warnings.filterwarnings('ignore', '.*too short*', ) slice_srs['Resultant'] = shock_spectrum( - df.iloc[window_start:window_end] * multiplier, + df.iloc[window_start:window_end], mode=mode, damp=damp, init_freq=init_freq, diff --git a/endaq/plot/plots.py b/endaq/plot/plots.py index 02d576e..f2b5cd1 100644 --- a/endaq/plot/plots.py +++ b/endaq/plot/plots.py @@ -538,12 +538,14 @@ def spectrum_over_time( :param df: the input dataframe with columns defining the frequency content, timestamps, values, and variables, see the following functions which provides outputs that would then be passed into this function as an input: + * :py:func:`~endaq.calc.fft.rolling_fft()` * :py:func:`~endaq.calc.psd.rolling_psd()` * :py:func:`~endaq.calc.shock.rolling_shock_spectrum()` * :py:func:`~endaq.batch.GetDataBuilder.add_psd()` * :py:func:`~endaq.batch.GetDataBuilder.add_pvss()` :param plot_type: the type of plot to display the spectrum, options are: + * `Heatmap`: a 2D visualization with the color defining the z value (the default) * `Surface`: similar to `Heatmap` but the z value is also projected "off the page" * `Waterfall`: distinct lines are plotted per time slice in a 3D view @@ -608,6 +610,47 @@ def spectrum_over_time( waterfall = endaq.plot.plots.spectrum_over_time(fft, plot_type='Waterfall', freq_max=200, var_to_process='Resultant') waterfall.show() + .. plotly:: + :fig-vars: heatmap, peak, surface, waterfall + + import pandas as pd + import endaq + + # Set Theme + endaq.plot.utilities.set_theme() + + # Get Vibration Data + df_vibe = pd.read_csv('https://info.endaq.com/hubfs/data/motorcycle-vibration-moving-frequency.csv', index_col=0) + + # Calculate a Rolling FFT + fft = endaq.calc.fft.rolling_fft(df_vibe, num_slices=200, add_resultant=True) + + # Visualize the Rolling FFT as a Heatmap + heatmap = endaq.plot.plots.spectrum_over_time(fft, plot_type='Heatmap', freq_max=200, var_to_process='Resultant') + heatmap.show() + + # Plot the Peak Frequency vs Time + peak = endaq.plot.plots.spectrum_over_time(fft, plot_type='Peak', freq_max=200, var_to_process='Resultant') + peak.show() + + # Visualize as a Surface Plot + surface = endaq.plot.plots.spectrum_over_time(fft, plot_type='Surface', freq_max=200, var_to_process='Resultant') + surface.show() + + # Visualize as a Waterfall + waterfall = endaq.plot.plots.spectrum_over_time(fft, plot_type='Waterfall', freq_max=200, var_to_process='Resultant') + waterfall.show() + + Here's another few examples with a longer dataset with DatetimeIndex of a car engine during a morning commute + + .. code:: python + + import pandas as pd + import endaq + + # Set Theme + endaq.plot.utilities.set_theme() + # Get a Longer Dataset with DatetimeIndex engine = endaq.ide.get_primary_sensor_data('https://info.endaq.com/hubfs/data/Commute.ide', measurement_type='accel', time_mode='datetime') @@ -646,7 +689,7 @@ def spectrum_over_time( waterfall2.show() .. plotly:: - :fig-vars: heatmap, peak, surface, waterfall, heatmap2, animation, lines, waterfall2 + :fig-vars: heatmap2, animation, lines, waterfall2 import pandas as pd import endaq @@ -654,28 +697,6 @@ def spectrum_over_time( # Set Theme endaq.plot.utilities.set_theme() - # Get Vibration Data - df_vibe = pd.read_csv('https://info.endaq.com/hubfs/data/motorcycle-vibration-moving-frequency.csv', index_col=0) - - # Calculate a Rolling FFT - fft = endaq.calc.fft.rolling_fft(df_vibe, num_slices=200, add_resultant=True) - - # Visualize the Rolling FFT as a Heatmap - heatmap = endaq.plot.plots.spectrum_over_time(fft, plot_type='Heatmap', freq_max=200, var_to_process='Resultant') - heatmap.show() - - # Plot the Peak Frequency vs Time - peak = endaq.plot.plots.spectrum_over_time(fft, plot_type='Peak', freq_max=200, var_to_process='Resultant') - peak.show() - - # Visualize as a Surface Plot - surface = endaq.plot.plots.spectrum_over_time(fft, plot_type='Surface', freq_max=200, var_to_process='Resultant') - surface.show() - - # Visualize as a Waterfall - waterfall = endaq.plot.plots.spectrum_over_time(fft, plot_type='Waterfall', freq_max=200, var_to_process='Resultant') - waterfall.show() - # Get a Longer Dataset with DatetimeIndex engine = endaq.ide.get_primary_sensor_data('https://info.endaq.com/hubfs/data/Commute.ide', measurement_type='accel', time_mode='datetime') @@ -731,9 +752,9 @@ def spectrum_over_time( # Filter frequency df = df.loc[df[freq_column] > 0.0] if freq_max is not None: - df = df.loc[df[freq_column] < freq_max] + df = df.loc[df[freq_column] <= freq_max] if freq_min is not None: - df = df.loc[df[freq_column] > freq_min] + df = df.loc[df[freq_column] >= freq_min] # Remove 0s df = df.loc[df[val_column] > 0] @@ -749,10 +770,10 @@ def spectrum_over_time( df_pivot = df.copy() first_step = df_pivot[freq_column].iloc[1] - df_pivot[freq_column].iloc[0] second_step = df_pivot[freq_column].iloc[2] - df_pivot[freq_column].iloc[1] - if first_step == second_step: + if np.isclose(first_step, second_step): round_freq = np.round(df_pivot[freq_column].min(), 0) df_pivot[freq_column] = np.round(df_pivot[freq_column].to_numpy() / round_freq, 0) * round_freq - df_pivot = df_pivot.pivot(columns=time_column, index=freq_column, values=val_column) + df_pivot = df_pivot.pivot_table(columns=time_column, index=freq_column, values=val_column) # Heatmap & Surface if plot_type in ["Heatmap", "Surface"]: diff --git a/tests/calc/test_fft.py b/tests/calc/test_fft.py index 7e87768..fd8c445 100644 --- a/tests/calc/test_fft.py +++ b/tests/calc/test_fft.py @@ -19,27 +19,27 @@ class TestFFT: @pytest.mark.parametrize("normalization, output, n", [ - (None, None, None), - ("unit", None, None), - ("forward", None, None), - ("backward", None, None), - ("ortho", None, None), - (None, "magnitude", None), - (None, "angle", None), - (None, "complex", None), + (None, None, None), + ("unit", None, None), + ("forward", None, None), + ("backward", None, None), + ("ortho", None, None), + (None, "magnitude", None), + (None, "angle", None), + (None, "complex", None), ], ) def test_fft(self, normalization, output, n): - data = np.array([0.0]*10 + [-0.5, -1.0, 1.0, 0.5] + [0.0]*10) + data = np.array([0.0] * 10 + [-0.5, -1.0, 1.0, 0.5] + [0.0] * 10) fxx = np.fft.fftshift(np.fft.fft(data)) n = len(data) - scales = {"unit": 2./n, "forward": 1./n, "backward": 1., "ortho": np.sqrt(1./n)} + scales = {"unit": 2. / n, "forward": 1. / n, "backward": 1., "ortho": np.sqrt(1. / n)} - scale = scales.get(normalization, 2./n) + scale = scales.get(normalization, 2. / n) fxx *= scale if output == "angle": @@ -62,12 +62,12 @@ def test_fft(self, normalization, output, n): @pytest.mark.parametrize( "normalization, output, n, raises", [ - ("sideways", None, None, pytest.raises(ValueError)), - (None, "octonian", None, pytest.raises(ValueError)), - (None, None, -1, pytest.raises(ValueError)), - (None, None, 0, pytest.raises(ValueError)), - (None, None, 0., pytest.raises(TypeError)), - (None, None, "five", pytest.raises(TypeError)), + ("sideways", None, None, pytest.raises(ValueError)), + (None, "octonian", None, pytest.raises(ValueError)), + (None, None, -1, pytest.raises(ValueError)), + (None, None, 0, pytest.raises(ValueError)), + (None, None, 0., pytest.raises(TypeError)), + (None, None, "five", pytest.raises(TypeError)), ], ) def test_fft_exceptions(self, normalization, output, n, raises): @@ -75,7 +75,6 @@ def test_fft_exceptions(self, normalization, output, n, raises): df = pd.DataFrame(np.arange(10)) with raises: - fft.fft(df, norm=normalization, output=output, nfft=n) @@ -83,27 +82,27 @@ class TestRFFT: @pytest.mark.parametrize("normalization, output, n", [ - (None, None, None), - ("unit", None, None), - ("forward", None, None), - ("backward", None, None), - ("ortho", None, None), - (None, "magnitude", None), - (None, "angle", None), - (None, "complex", None), + (None, None, None), + ("unit", None, None), + ("forward", None, None), + ("backward", None, None), + ("ortho", None, None), + (None, "magnitude", None), + (None, "angle", None), + (None, "complex", None), ], ) def test_rfft(self, normalization, output, n): - data = np.array([0.0]*10 + [-0.5, -1.0, 1.0, 0.5] + [0.0]*10) + data = np.array([0.0] * 10 + [-0.5, -1.0, 1.0, 0.5] + [0.0] * 10) fxx = np.fft.rfft(data) n = len(data) - scales = {"unit": 2./n, "forward": 1./n, "backward": 1., "ortho": np.sqrt(1./n)} + scales = {"unit": 2. / n, "forward": 1. / n, "backward": 1., "ortho": np.sqrt(1. / n)} - scale = scales.get(normalization, 2./n) + scale = scales.get(normalization, 2. / n) fxx *= scale if output == "angle": @@ -126,12 +125,12 @@ def test_rfft(self, normalization, output, n): @pytest.mark.parametrize( "normalization, output, n, raises", [ - ("sideways", None, None, pytest.raises(ValueError)), - (None, "octonian", None, pytest.raises(ValueError)), - (None, None, -1, pytest.raises(ValueError)), - (None, None, 0, pytest.raises(ValueError)), - (None, None, 0., pytest.raises(TypeError)), - (None, None, "five", pytest.raises(TypeError)), + ("sideways", None, None, pytest.raises(ValueError)), + (None, "octonian", None, pytest.raises(ValueError)), + (None, None, -1, pytest.raises(ValueError)), + (None, None, 0, pytest.raises(ValueError)), + (None, None, 0., pytest.raises(TypeError)), + (None, None, "five", pytest.raises(TypeError)), ], ) def test_fft_exceptions(self, normalization, output, n, raises): @@ -139,5 +138,40 @@ def test_fft_exceptions(self, normalization, output, n, raises): df = pd.DataFrame(np.arange(10)) with raises: - fft.rfft(df, norm=normalization, output=output, nfft=n) + + +def test_rolling_fft(): + # Build Time Array + time = np.arange(1000) / 1000 + + # Build Dataframe with Noise + df_accel = pd.DataFrame({ + 'time': np.concatenate((time, time + 1.0)), + 'A': np.concatenate(( + np.sin(2 * np.pi * 10 * time), + np.sin(2 * np.pi * 13 * time)) + ) + }).set_index('time') + + # Add Random + df_accel.A = df_accel.A + np.random.random(2000) * 2 - 1 + + # Get the Rolling FFT Dataframe + df_rolling_fft = fft.rolling_fft( + df_accel, + num_slices=2, + add_resultant=False + ) + + # Get the Individual FFTs + first_fft = fft.aggregate_fft(df_accel[:1.0]) + second_fft = fft.aggregate_fft(df_accel[1.0:]) + + # Do Assertions + npt.assert_almost_equal( + df_rolling_fft['value'].iloc[:501].to_numpy(), + first_fft['A'].to_numpy()) + npt.assert_almost_equal( + df_rolling_fft['value'].iloc[-501:].to_numpy(), + second_fft['A'].to_numpy()) diff --git a/tests/calc/test_utils.py b/tests/calc/test_utils.py index 60b42af..e7ad432 100644 --- a/tests/calc/test_utils.py +++ b/tests/calc/test_utils.py @@ -67,3 +67,35 @@ def test_uniform_resample_is_uniform_for_datetime64ns(): times = np.arange('2021-02', '2021-03', np.timedelta64(1, 's'), dtype='datetime64[ns]') df = pd.DataFrame(np.arange(len(times)), index=times) assert len(np.unique(np.diff(utils.resample(df).index))) == 1 + + +def test_rolling_slice_definitions(): + # Build dataframe with 1 second of data + df = pd.DataFrame({ + 'time': np.arange(1000) / 1000, + 'A': np.ones(1000) + }).set_index('time') + + indexes, slice_width, num, length = utils._rolling_slice_definitions( + df=df, + num_slices=2 + ) + assert slice_width == 0.5 + + indexes, slice_width, num, length = utils._rolling_slice_definitions( + df=df, + index_values=[0.1, 0.9] + ) + assert indexes[1] == 900 + + df.index = pd.to_datetime(df.index, unit='s') + indexes, slice_width, num, length = utils._rolling_slice_definitions( + df=df, + ) + assert num == 5 + + indexes, slice_width, num, length = utils._rolling_slice_definitions( + df=df, + index_values=pd.DatetimeIndex(['1970-01-01 00:00:00.9951']) + ) + assert indexes[0] == 995 diff --git a/tests/plot/test_plots.py b/tests/plot/test_plots.py new file mode 100644 index 0000000..f9b8e00 --- /dev/null +++ b/tests/plot/test_plots.py @@ -0,0 +1,36 @@ +import pytest + +import numpy as np +import pandas as pd + +from endaq import plot + + +def test_spectrums_over_time(): + freqs = np.arange(10) + times = np.zeros_like(freqs) + + df = pd.DataFrame({ + 'frequency (Hz)': np.concatenate((freqs, freqs, freqs)), + 'timestamp': np.concatenate((times, times + 30, times + 60)), + 'value': np.random.random(30) + }) + df['variable'] = 'X' + + fig = plot.spectrum_over_time(df, 'Heatmap') + assert fig['data'][0]['type'] == 'heatmap' + + fig = plot.spectrum_over_time(df, 'Surface') + assert fig['data'][0]['type'] == 'surface' + + fig = plot.spectrum_over_time(df, 'Waterfall') + assert fig['data'][0]['type'] == 'scatter3d' + + fig = plot.spectrum_over_time(df, 'Animation') + assert len(fig['data']) == 4 + + fig = plot.spectrum_over_time(df, 'Lines') + assert len(fig['data']) == 9 + + fig = plot.spectrum_over_time(df, 'Peak') + assert fig['data'][0]['marker']['symbol'] == 'circle' From 760f41dff3d03aebe232212dbf9b8757cd69bdc8 Mon Sep 17 00:00:00 2001 From: "MIDE\\shanly" Date: Fri, 22 Apr 2022 15:50:01 -0400 Subject: [PATCH 22/24] used scipy.spectrogram for evenly spaced windows --- endaq/calc/fft.py | 52 ++++++----- endaq/calc/psd.py | 213 +++++++++++++++++++++++++++++++++++++++------- 2 files changed, 216 insertions(+), 49 deletions(-) diff --git a/endaq/calc/fft.py b/endaq/calc/fft.py index 76b24b7..eaa649a 100644 --- a/endaq/calc/fft.py +++ b/endaq/calc/fft.py @@ -69,6 +69,9 @@ def rolling_fft( Surface plots, and Animations """ + use_spectrogram = True + if (indexes is not None) | (index_values is not None): + use_spectrogram = False indexes, slice_width, num, length = utils._rolling_slice_definitions( df, @@ -78,25 +81,36 @@ def rolling_fft( slice_width=slice_width ) - # Loop through and compute fft - fft_df = pd.DataFrame() - for i in indexes: - window_start = max(0, i - num) - window_end = min(length, i + num) - with warnings.catch_warnings(): - if disable_warnings: - warnings.filterwarnings('ignore', '.*greater than.*') - slice_fft = aggregate_fft( - df.iloc[window_start:window_end], - bin_width=bin_width, - **kwargs - ) - if add_resultant: - slice_fft['Resultant'] = slice_fft.pow(2).sum(axis=1).pow(0.5) - - slice_fft = slice_fft.reset_index().melt(id_vars=slice_fft.index.name) - slice_fft['timestamp'] = df.index[i] - fft_df = pd.concat([fft_df, slice_fft]) + if use_spectrogram: + fft_df = psd.spectrogram( + df=df, + num_slices=num_slices, + bin_width=bin_width, + scaling='unit', + add_resultant=add_resultant, + disable_warnings=disable_warnings, + **kwargs + ) + else: + # Loop through and compute fft + fft_df = pd.DataFrame() + for i in indexes: + window_start = max(0, i - num) + window_end = min(length, i + num) + with warnings.catch_warnings(): + if disable_warnings: + warnings.filterwarnings('ignore', '.*greater than.*') + slice_fft = aggregate_fft( + df.iloc[window_start:window_end], + bin_width=bin_width, + **kwargs + ) + if add_resultant: + slice_fft['Resultant'] = slice_fft.pow(2).sum(axis=1).pow(0.5) + + slice_fft = slice_fft.reset_index().melt(id_vars=slice_fft.index.name) + slice_fft['timestamp'] = df.index[i] + fft_df = pd.concat([fft_df, slice_fft]) return fft_df diff --git a/endaq/calc/psd.py b/endaq/calc/psd.py index cf05a7f..de4e0f9 100644 --- a/endaq/calc/psd.py +++ b/endaq/calc/psd.py @@ -4,6 +4,7 @@ import typing import warnings +import datetime import numpy as np import pandas as pd @@ -391,6 +392,9 @@ def rolling_psd( Surface plots, and Animations """ + use_spectrogram = True + if (indexes is not None) | (index_values is not None): + use_spectrogram = False indexes, slice_width, num, length = utils._rolling_slice_definitions( df, @@ -400,46 +404,195 @@ def rolling_psd( slice_width=slice_width ) - # Allow for RMS calculations - exp = 1 - if scaling == "rms": - scaling = "parseval" - exp = 0.5 + if use_spectrogram: + psd_df = spectrogram( + df=df, + bin_width=bin_width, + num_slices=num_slices, + scaling=scaling, + octave_bins=octave_bins, + fstart=fstart, + freq_splits=freq_splits, + add_resultant=add_resultant, + disable_warnings=disable_warnings, + **kwargs + ) + else: + # Allow for RMS calculations + exp = 1 + if scaling == "rms": + scaling = "parseval" + exp = 0.5 + agg = "sum" + + # Define bin_width if octave spacing + if octave_bins is not None: + bin_width = 1 / slice_width + + # Loop through and compute PSD + psd_df = pd.DataFrame() + for i in indexes: + window_start = max(0, i - num) + window_end = min(length, i + num) + with warnings.catch_warnings(): + if disable_warnings: + warnings.filterwarnings('ignore', '.*greater than.*', ) + slice_psd = welch( + df.iloc[window_start:window_end], + bin_width=bin_width, + scaling=scaling, + **kwargs + ) + + if octave_bins is not None or freq_splits is not None: + with warnings.catch_warnings(): + if disable_warnings: + warnings.filterwarnings('ignore', '.*empty frequency.*') + if freq_splits is None: + slice_psd = to_octave(slice_psd, octave_bins=octave_bins, fstart=fstart, agg=agg) + else: + slice_psd = to_jagged(slice_psd, freq_splits=freq_splits, agg=agg) + + if add_resultant: + slice_psd['Resultant'] = slice_psd.sum(axis=1) + + slice_psd = (slice_psd ** exp).reset_index().melt(id_vars=slice_psd.index.name) + slice_psd['timestamp'] = df.index[i] + psd_df = pd.concat([psd_df, slice_psd]) + + return psd_df + + +def spectrogram( + df: pd.DataFrame, + num_slices: int = 100, + scaling: typing.Literal[None, "density", "spectrum", "parseval", "unit", "rms"] = None, + bin_width: float = 1.0, + octave_bins: float = None, + fstart: float = 1.0, + agg: typing.Union[ + typing.Literal["mean", "sum"], + typing.Callable[[np.ndarray, int], float], + ] = "mean", + freq_splits: np.array = None, + add_resultant: bool = True, + disable_warnings: bool = True, + **kwargs, +) -> pd.DataFrame: + """ + Compute a spectrogram for a time series data set using :py:func:`scipy.signal.spectrogram()` + + :param df: the input dataframe with an index defining the time in seconds or datetime + :param num_slices: the number of slices to split the time series into, default is 100 + :param scaling: the scaling of the output; `"density"` & `"spectrum"` + correspond to the same options in ``scipy.signal.welch``; `"parseval"` + will maintain the "energy" between the input & output, s.t. + ``welch(df, scaling="parseval").sum(axis="rows")`` is roughly equal to + ``df.abs().pow(2).sum(axis="rows")``; + `"unit"` will maintain the units and scale of the input data. `"rms"` will use `"parseval"` for the PSD + calculations and set `agg` to "sum", but then take the square root at the end + :param bin_width: the bin width or resolution in Hz for the PSD or FFT, defaults to 1 + :param octave_bins: the number of frequency bins in each octave, defaults to `None` + :param fstart: the lowest frequency for an octave spaced output, defaults to 1 + :param agg: the method for aggregating values into bins (only used if converting to octave or jagged); `'mean'` preserves + the PSD's area-under-the-curve, `'sum'` preserves the PSD's "energy" + :param freq_splits: the boundaries of the frequency bins to pass to :py:func:`~endaq.calc.psd.to_jagged()` + :param add_resultant: if `True` (default) the root sum of squares of each PSD column will + also be computed + :param disable_warnings: if `True` (default) it disables the warnings in helper functions + :param kwargs: Other parameters to pass directly to :py:func:`scipy.signal.spectrogram()` + :return: a dataframe containing all the spectrograms "melted" with columns defining the value, frequency, timeslice, + and original column from the input dataframe + + .. seealso:: + + - `SciPy Spectrogram method `_ + Documentation for the function wrapped internally. + """ + + # Apply Scaling + if scaling == "parseval": + kwargs["scaling"] = "density" + elif scaling == "unit": + kwargs["scaling"] = "density" + elif scaling == "rms": + kwargs["scaling"] = "density" agg = "sum" + elif scaling is not None: + kwargs["scaling"] = scaling + + # Windowing, Default to Boxcar + if "window" not in kwargs: + kwargs["window"] = "boxcar" + + # Determine Scipy.Spectrogram Stats + fs = 1 / utils.sample_spacing(df) + nperseg = int(fs / bin_width) + num_f = int(len(df) / fs / bin_width) + noverlap = nperseg - (nperseg * num_f / num_slices) + + # Loop Through and Compute Spectrogram + spec_df = pd.DataFrame() + if add_resultant: + res_df = pd.DataFrame() + for c in df.columns: + f, t, spec = scipy.signal.spectrogram( + x=df[c].to_numpy(), + fs=fs, + nperseg=nperseg, + mode='psd', + noverlap=noverlap, + **kwargs + ) - # Define bin_width if octave spacing - if octave_bins is not None: - bin_width = 1 / slice_width - - # Loop through and compute PSD - psd_df = pd.DataFrame() - for i in indexes: - window_start = max(0, i - num) - window_end = min(length, i + num) - with warnings.catch_warnings(): - if disable_warnings: - warnings.filterwarnings('ignore', '.*greater than.*', ) - slice_psd = welch( - df.iloc[window_start:window_end], - bin_width=bin_width, - scaling=scaling, - **kwargs - ) + # create dataframe + spec_df_col = pd.DataFrame( + data=spec, + columns=t, + index=f + ) + spec_df_col.index.name = 'frequency (Hz)' + # Do octave spacing if octave_bins is not None or freq_splits is not None: with warnings.catch_warnings(): if disable_warnings: warnings.filterwarnings('ignore', '.*empty frequency.*') if freq_splits is None: - slice_psd = to_octave(slice_psd, octave_bins=octave_bins, fstart=fstart, agg=agg) + spec_df_col = to_octave(spec_df_col, octave_bins=octave_bins, fstart=fstart, agg=agg) else: - slice_psd = to_jagged(slice_psd, freq_splits=freq_splits, agg=agg) + spec_df_col = to_jagged(spec_df_col, freq_splits=freq_splits, agg=agg) + # Melt & concatenate dataframe + spec_df_col = spec_df_col.reset_index().melt(id_vars='frequency (Hz)') + spec_df_col.columns = ['frequency (Hz)', 'timestamp', 'value'] + spec_df_col['variable'] = c + spec_df = pd.concat([spec_df, spec_df_col]) + + # add resultant if add_resultant: - slice_psd['Resultant'] = slice_psd.sum(axis=1) + if len(res_df) == 0: + res_df = spec_df_col + res_df['variable'] = "Resultant" + else: + res_df.value += spec_df_col.value + + # Add resultant + if add_resultant: + spec_df = pd.concat([spec_df, res_df]) + + # Update timestamps to be based on input, deal with datetime + if isinstance(df.index[0], datetime.datetime): + spec_df.timestamp = pd.to_timedelta(spec_df.timestamp.astype(float), unit='s') + df.index[0] + else: + spec_df.timestamp = spec_df.timestamp.astype(float) + df.index[0] - slice_psd = (slice_psd ** exp).reset_index().melt(id_vars=slice_psd.index.name) - slice_psd['timestamp'] = df.index[i] - psd_df = pd.concat([psd_df, slice_psd]) + # Apply scales + if scaling == "parseval": + spec_df.value = spec_df.value * f[1] + if scaling == "rms": + spec_df.value = spec_df.value ** 0.5 + elif scaling == "unit": + spec_df.value = (spec_df.value * 2 * f[1]) ** 0.5 - return psd_df + return spec_df From 179fc4b5f81dd6713764d4fffa5423dbaa29a574 Mon Sep 17 00:00:00 2001 From: "MIDE\\shanly" Date: Thu, 28 Apr 2022 16:56:32 -0400 Subject: [PATCH 23/24] added tests for all new functions, updated test_fft for correct formating --- endaq/calc/psd.py | 4 +- endaq/calc/shock.py | 2 +- tests/calc/test_fft.py | 123 ++++++++++++++++++++++----------------- tests/calc/test_psd.py | 49 ++++++++++++++++ tests/calc/test_shock.py | 63 ++++++++++++++++++++ tests/plot/test_plots.py | 53 ++++++++++++----- 6 files changed, 227 insertions(+), 67 deletions(-) diff --git a/endaq/calc/psd.py b/endaq/calc/psd.py index de4e0f9..0508794 100644 --- a/endaq/calc/psd.py +++ b/endaq/calc/psd.py @@ -405,6 +405,8 @@ def rolling_psd( ) if use_spectrogram: + if "window" not in kwargs: + kwargs["window"] = "hann" psd_df = spectrogram( df=df, bin_width=bin_width, @@ -529,7 +531,7 @@ def spectrogram( fs = 1 / utils.sample_spacing(df) nperseg = int(fs / bin_width) num_f = int(len(df) / fs / bin_width) - noverlap = nperseg - (nperseg * num_f / num_slices) + noverlap = nperseg - (nperseg * num_f // num_slices) # Loop Through and Compute Spectrogram spec_df = pd.DataFrame() diff --git a/endaq/calc/shock.py b/endaq/calc/shock.py index 0bede15..8b39965 100644 --- a/endaq/calc/shock.py +++ b/endaq/calc/shock.py @@ -444,7 +444,7 @@ def shock_spectrum( function. """ if freqs is None: - freqs = utils.logfreqs(accel, init_freq = init_freq, bins_per_octave = bins_per_octave) + freqs = utils.logfreqs(accel, init_freq=init_freq, bins_per_octave=bins_per_octave) if two_sided and aggregate_axes: raise ValueError("cannot enable both options `two_sided` and `aggregate_axes`") freqs = np.asarray(freqs) diff --git a/tests/calc/test_fft.py b/tests/calc/test_fft.py index fd8c445..e0a3496 100644 --- a/tests/calc/test_fft.py +++ b/tests/calc/test_fft.py @@ -19,27 +19,27 @@ class TestFFT: @pytest.mark.parametrize("normalization, output, n", [ - (None, None, None), - ("unit", None, None), - ("forward", None, None), - ("backward", None, None), - ("ortho", None, None), - (None, "magnitude", None), - (None, "angle", None), - (None, "complex", None), + (None, None, None), + ("unit", None, None), + ("forward", None, None), + ("backward", None, None), + ("ortho", None, None), + (None, "magnitude", None), + (None, "angle", None), + (None, "complex", None), ], ) def test_fft(self, normalization, output, n): - data = np.array([0.0] * 10 + [-0.5, -1.0, 1.0, 0.5] + [0.0] * 10) + data = np.array([0.0]*10 + [-0.5, -1.0, 1.0, 0.5] + [0.0]*10) fxx = np.fft.fftshift(np.fft.fft(data)) n = len(data) - scales = {"unit": 2. / n, "forward": 1. / n, "backward": 1., "ortho": np.sqrt(1. / n)} + scales = {"unit": 2./n, "forward": 1./n, "backward": 1., "ortho": np.sqrt(1./n)} - scale = scales.get(normalization, 2. / n) + scale = scales.get(normalization, 2./n) fxx *= scale if output == "angle": @@ -62,12 +62,12 @@ def test_fft(self, normalization, output, n): @pytest.mark.parametrize( "normalization, output, n, raises", [ - ("sideways", None, None, pytest.raises(ValueError)), - (None, "octonian", None, pytest.raises(ValueError)), - (None, None, -1, pytest.raises(ValueError)), - (None, None, 0, pytest.raises(ValueError)), - (None, None, 0., pytest.raises(TypeError)), - (None, None, "five", pytest.raises(TypeError)), + ("sideways", None, None, pytest.raises(ValueError)), + (None, "octonian", None, pytest.raises(ValueError)), + (None, None, -1, pytest.raises(ValueError)), + (None, None, 0, pytest.raises(ValueError)), + (None, None, 0., pytest.raises(TypeError)), + (None, None, "five", pytest.raises(TypeError)), ], ) def test_fft_exceptions(self, normalization, output, n, raises): @@ -75,6 +75,7 @@ def test_fft_exceptions(self, normalization, output, n, raises): df = pd.DataFrame(np.arange(10)) with raises: + fft.fft(df, norm=normalization, output=output, nfft=n) @@ -82,27 +83,27 @@ class TestRFFT: @pytest.mark.parametrize("normalization, output, n", [ - (None, None, None), - ("unit", None, None), - ("forward", None, None), - ("backward", None, None), - ("ortho", None, None), - (None, "magnitude", None), - (None, "angle", None), - (None, "complex", None), + (None, None, None), + ("unit", None, None), + ("forward", None, None), + ("backward", None, None), + ("ortho", None, None), + (None, "magnitude", None), + (None, "angle", None), + (None, "complex", None), ], ) def test_rfft(self, normalization, output, n): - data = np.array([0.0] * 10 + [-0.5, -1.0, 1.0, 0.5] + [0.0] * 10) + data = np.array([0.0]*10 + [-0.5, -1.0, 1.0, 0.5] + [0.0]*10) fxx = np.fft.rfft(data) n = len(data) - scales = {"unit": 2. / n, "forward": 1. / n, "backward": 1., "ortho": np.sqrt(1. / n)} + scales = {"unit": 2./n, "forward": 1./n, "backward": 1., "ortho": np.sqrt(1./n)} - scale = scales.get(normalization, 2. / n) + scale = scales.get(normalization, 2./n) fxx *= scale if output == "angle": @@ -125,12 +126,12 @@ def test_rfft(self, normalization, output, n): @pytest.mark.parametrize( "normalization, output, n, raises", [ - ("sideways", None, None, pytest.raises(ValueError)), - (None, "octonian", None, pytest.raises(ValueError)), - (None, None, -1, pytest.raises(ValueError)), - (None, None, 0, pytest.raises(ValueError)), - (None, None, 0., pytest.raises(TypeError)), - (None, None, "five", pytest.raises(TypeError)), + ("sideways", None, None, pytest.raises(ValueError)), + (None, "octonian", None, pytest.raises(ValueError)), + (None, None, -1, pytest.raises(ValueError)), + (None, None, 0, pytest.raises(ValueError)), + (None, None, 0., pytest.raises(TypeError)), + (None, None, "five", pytest.raises(TypeError)), ], ) def test_fft_exceptions(self, normalization, output, n, raises): @@ -141,7 +142,8 @@ def test_fft_exceptions(self, normalization, output, n, raises): fft.rfft(df, norm=normalization, output=output, nfft=n) -def test_rolling_fft(): +@pytest.fixture +def df_test(): # Build Time Array time = np.arange(1000) / 1000 @@ -157,21 +159,38 @@ def test_rolling_fft(): # Add Random df_accel.A = df_accel.A + np.random.random(2000) * 2 - 1 - # Get the Rolling FFT Dataframe - df_rolling_fft = fft.rolling_fft( - df_accel, - num_slices=2, - add_resultant=False - ) + return df_accel + + +class TestRollingFFT: + + def test_using_spectrogram(self, df_test): + df_rolling_fft = fft.rolling_fft( + df_test, + num_slices=2, + add_resultant=False + ) + times = df_rolling_fft.timestamp.unique() + + npt.assert_almost_equal( + df_rolling_fft[df_rolling_fft.timestamp == times[0]].value.to_numpy(), + fft.aggregate_fft(df_test[:1.0])['A'].to_numpy()) + npt.assert_almost_equal( + df_rolling_fft[df_rolling_fft.timestamp == times[1]].value.to_numpy(), + fft.aggregate_fft(df_test[1.0:])['A'].to_numpy()) + + def test_defined_slices(self, df_test): + df_rolling_fft = fft.rolling_fft( + df_test, + index_values=[1.0, 1.5], + slice_width=0.5, + add_resultant=False, + ) - # Get the Individual FFTs - first_fft = fft.aggregate_fft(df_accel[:1.0]) - second_fft = fft.aggregate_fft(df_accel[1.0:]) - - # Do Assertions - npt.assert_almost_equal( - df_rolling_fft['value'].iloc[:501].to_numpy(), - first_fft['A'].to_numpy()) - npt.assert_almost_equal( - df_rolling_fft['value'].iloc[-501:].to_numpy(), - second_fft['A'].to_numpy()) + # Do Assertions + npt.assert_almost_equal( + df_rolling_fft[df_rolling_fft.timestamp == 1.0].value.to_numpy(), + fft.aggregate_fft(df_test.iloc[750:1250])['A'].to_numpy()) + npt.assert_almost_equal( + df_rolling_fft[df_rolling_fft.timestamp == 1.5].value.to_numpy(), + fft.aggregate_fft(df_test.iloc[1250:1750])['A'].to_numpy()) diff --git a/tests/calc/test_psd.py b/tests/calc/test_psd.py index 13c003c..c38207e 100644 --- a/tests/calc/test_psd.py +++ b/tests/calc/test_psd.py @@ -159,3 +159,52 @@ def test_to_octave(psd_df, agg, expt_f, expt_array): calc_df = psd.to_octave(psd_df, fstart=1, octave_bins=1, agg=agg) assert calc_df.index.to_numpy().tolist() == expt_f assert calc_df.to_numpy().flatten().tolist() == expt_array + + +@pytest.fixture +def df_test(): + # Build Time Array + time = np.arange(1000) / 1000 + + # Build Dataframe with Noise + df_accel = pd.DataFrame({ + 'time': np.concatenate((time, time + 1.0)), + 'A': np.random.random(2000) + }).set_index('time') + + return df_accel + + +@pytest.mark.filterwarnings("ignore:empty frequency bins:RuntimeWarning") +class TestRollingPSD: + + def test_using_spectrogram(self, df_test): + df_rolling_psd = psd.rolling_psd( + df_test, + num_slices=2, + add_resultant=False + ) + times = df_rolling_psd.timestamp.unique() + + np.testing.assert_almost_equal( + df_rolling_psd[df_rolling_psd.timestamp == times[0]].value.to_numpy(), + psd.welch(df_test[:1.0])['A'].to_numpy()) + np.testing.assert_almost_equal( + df_rolling_psd[df_rolling_psd.timestamp == times[1]].value.to_numpy(), + psd.welch(df_test[1.0:])['A'].to_numpy()) + + def test_defined_slices(self, df_test): + df_rolling_psd = psd.rolling_psd( + df_test, + index_values=[1.0, 1.5], + slice_width=0.5, + add_resultant=False, + ) + + # Do Assertions + np.testing.assert_almost_equal( + df_rolling_psd[df_rolling_psd.timestamp == 1.0].value.to_numpy(), + psd.welch(df_test.iloc[750:1250])['A'].to_numpy()) + np.testing.assert_almost_equal( + df_rolling_psd[df_rolling_psd.timestamp == 1.5].value.to_numpy(), + psd.welch(df_test.iloc[1250:1750])['A'].to_numpy()) \ No newline at end of file diff --git a/tests/calc/test_shock.py b/tests/calc/test_shock.py index 37cdcf2..00ec0a6 100644 --- a/tests/calc/test_shock.py +++ b/tests/calc/test_shock.py @@ -486,3 +486,66 @@ def test_tuple_like(self): ampl, T = env_half_sine assert ampl == mock.sentinel.amplitude assert T == mock.sentinel.duration + +@pytest.fixture +def df_test(): + half = shock.HalfSineWavePulse( + amplitude=pd.Series([100, 50]), + duration=pd.Series([.1, .05]) + ) + + df_accel = pd.concat([ + half.to_time_series(tstart=0, tstop=1), + half.to_time_series(tstart=1, tstop=2), + ]) + df_accel.columns = ['A', 'B'] + + return df_accel + + +class TestRollingShockSpectrum: + + def test_even_slices(self, df_test): + df_rolling_srs = shock.rolling_shock_spectrum( + df_test, + num_slices=2, + add_resultant=False + ) + + # Get the Individual SRS + first_srs = shock.shock_spectrum(df_test[:1.0]) + second_srs = shock.shock_spectrum(df_test[1.0:]) + + # Do Assertions + npt.assert_almost_equal( + df_rolling_srs[(df_rolling_srs.variable == 'A') & (df_rolling_srs.timestamp == 0.5)]['value'].to_numpy(), + first_srs['A'].to_numpy()) + npt.assert_almost_equal( + df_rolling_srs[(df_rolling_srs.variable == 'B') & (df_rolling_srs.timestamp == 0.5)]['value'].to_numpy(), + first_srs['B'].to_numpy()) + npt.assert_almost_equal( + df_rolling_srs[(df_rolling_srs.variable == 'A') & (df_rolling_srs.timestamp != 0.5)]['value'].to_numpy(), + second_srs['A'].to_numpy()) + npt.assert_almost_equal( + df_rolling_srs[(df_rolling_srs.variable == 'B') & (df_rolling_srs.timestamp != 0.5)]['value'].to_numpy(), + second_srs['B'].to_numpy()) + + def test_defined_slices(self, df_test): + df_rolling_srs = shock.rolling_shock_spectrum( + df_test, + index_values=[1.0], + bins_per_octave=6, + slice_width=0.5, + mode='pvss' + ) + + npt.assert_almost_equal( + df_rolling_srs[df_rolling_srs.variable == 'Resultant']['value'].to_numpy(), + shock.shock_spectrum( + df_test[0.75:1.2499], + bins_per_octave=6, + mode='pvss', + aggregate_axes=True + )['resultant'].to_numpy() + ) + diff --git a/tests/plot/test_plots.py b/tests/plot/test_plots.py index f9b8e00..f87756c 100644 --- a/tests/plot/test_plots.py +++ b/tests/plot/test_plots.py @@ -6,7 +6,8 @@ from endaq import plot -def test_spectrums_over_time(): +@pytest.fixture(params=["float", "datetime"]) +def generate_dataframe(request): freqs = np.arange(10) times = np.zeros_like(freqs) @@ -17,20 +18,46 @@ def test_spectrums_over_time(): }) df['variable'] = 'X' - fig = plot.spectrum_over_time(df, 'Heatmap') - assert fig['data'][0]['type'] == 'heatmap' + if request.param == 'datetime': + df.timestamp = pd.to_datetime(df.timestamp, unit='D', utc=True) + + return df + + +class TestSpectrumOverTime: + + def test_heatmap(self, generate_dataframe): + fig = plot.spectrum_over_time(generate_dataframe, 'Heatmap') + assert fig['data'][0]['type'] == 'heatmap' - fig = plot.spectrum_over_time(df, 'Surface') - assert fig['data'][0]['type'] == 'surface' + def test_surface(self, generate_dataframe): + fig = plot.spectrum_over_time(generate_dataframe, 'Surface') + assert fig['data'][0]['type'] == 'surface' - fig = plot.spectrum_over_time(df, 'Waterfall') - assert fig['data'][0]['type'] == 'scatter3d' + def test_waterfall(self, generate_dataframe): + fig = plot.spectrum_over_time(generate_dataframe, 'Waterfall') + assert fig['data'][0]['type'] == 'scatter3d' - fig = plot.spectrum_over_time(df, 'Animation') - assert len(fig['data']) == 4 + def test_animation(self, generate_dataframe): + fig = plot.spectrum_over_time(generate_dataframe, 'Animation') + assert len(fig['data']) == 4 - fig = plot.spectrum_over_time(df, 'Lines') - assert len(fig['data']) == 9 + def test_lines(self, generate_dataframe): + fig = plot.spectrum_over_time(generate_dataframe, 'Lines') + assert len(fig['data']) == 9 + + def test_peak(self, generate_dataframe): + fig = plot.spectrum_over_time(generate_dataframe, 'Peak') + assert fig['data'][0]['marker']['symbol'] == 'circle' + + +def test_octave_spectrogram(): + time = np.arange(1000) / 1000 + df_accel = pd.DataFrame({ + 'time': np.concatenate((time, time + 1.0, time + 2.0)), + 'A': np.random.random(3000) + }).set_index('time') + + df, fig = plot.octave_spectrogram(df_accel, window=0.1, max_freq=50) + assert fig['data'][0]['type'] == 'heatmap' - fig = plot.spectrum_over_time(df, 'Peak') - assert fig['data'][0]['marker']['symbol'] == 'circle' From d51ec45b0a0570ea3962d4f3e0f67762585d4f54 Mon Sep 17 00:00:00 2001 From: "MIDE\\shanly" Date: Fri, 29 Apr 2022 09:31:47 -0400 Subject: [PATCH 24/24] updated warning filterings --- tests/calc/test_fft.py | 1 + tests/calc/test_psd.py | 1 + tests/calc/test_shock.py | 13 ++++++++----- 3 files changed, 10 insertions(+), 5 deletions(-) diff --git a/tests/calc/test_fft.py b/tests/calc/test_fft.py index e0a3496..45e31af 100644 --- a/tests/calc/test_fft.py +++ b/tests/calc/test_fft.py @@ -179,6 +179,7 @@ def test_using_spectrogram(self, df_test): df_rolling_fft[df_rolling_fft.timestamp == times[1]].value.to_numpy(), fft.aggregate_fft(df_test[1.0:])['A'].to_numpy()) + @pytest.mark.filterwarnings("ignore::UserWarning") def test_defined_slices(self, df_test): df_rolling_fft = fft.rolling_fft( df_test, diff --git a/tests/calc/test_psd.py b/tests/calc/test_psd.py index c38207e..663bb03 100644 --- a/tests/calc/test_psd.py +++ b/tests/calc/test_psd.py @@ -193,6 +193,7 @@ def test_using_spectrogram(self, df_test): df_rolling_psd[df_rolling_psd.timestamp == times[1]].value.to_numpy(), psd.welch(df_test[1.0:])['A'].to_numpy()) + @pytest.mark.filterwarnings("ignore::UserWarning") def test_defined_slices(self, df_test): df_rolling_psd = psd.rolling_psd( df_test, diff --git a/tests/calc/test_shock.py b/tests/calc/test_shock.py index 00ec0a6..cc7e348 100644 --- a/tests/calc/test_shock.py +++ b/tests/calc/test_shock.py @@ -509,12 +509,13 @@ def test_even_slices(self, df_test): df_rolling_srs = shock.rolling_shock_spectrum( df_test, num_slices=2, - add_resultant=False + add_resultant=False, + init_freq=2 ) # Get the Individual SRS - first_srs = shock.shock_spectrum(df_test[:1.0]) - second_srs = shock.shock_spectrum(df_test[1.0:]) + first_srs = shock.shock_spectrum(df_test[:1.0], init_freq=2) + second_srs = shock.shock_spectrum(df_test[1.0:], init_freq=2) # Do Assertions npt.assert_almost_equal( @@ -536,7 +537,8 @@ def test_defined_slices(self, df_test): index_values=[1.0], bins_per_octave=6, slice_width=0.5, - mode='pvss' + mode='pvss', + init_freq=4 ) npt.assert_almost_equal( @@ -545,7 +547,8 @@ def test_defined_slices(self, df_test): df_test[0.75:1.2499], bins_per_octave=6, mode='pvss', - aggregate_axes=True + aggregate_axes=True, + init_freq=4 )['resultant'].to_numpy() )