"""FFT plotting and digital lowpass filtering utilities."""
from __future__ import annotations
import numpy as np
import pandas as pd
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from scipy.fft import rfft, rfftfreq
from scipy.signal import butter, sosfiltfilt
from ._utils import detect_time_divisor, get_time_array
# ---------------------------------------------------------------------------
# Helpers
# ---------------------------------------------------------------------------
def _ensure_list_df(dfs):
if isinstance(dfs, pd.DataFrame):
return [dfs], True
return list(dfs), False
def _ensure_list_str(cols):
if isinstance(cols, str):
return [cols]
return list(cols)
_NEVER_FILTER = frozenset(
{
"pcm.vnav.yawPitchRoll.yaw",
}
)
def _check_blocked(col_list: list[str]) -> list[str]:
"""Remove blocked columns with a warning, return the rest."""
safe = []
for c in col_list:
if c in _NEVER_FILTER:
print(
f" [filter] Skipping {c} — blocked "
f"(wrapping discontinuities). "
f"Filter after unwrap_degrees() if needed."
)
else:
safe.append(c)
return safe
# ---------------------------------------------------------------------------
# Digital lowpass filter
# ---------------------------------------------------------------------------
def _save_original(df: pd.DataFrame, col: str) -> None:
"""Save ``col_original`` if it doesn't already exist."""
orig = col + "_original"
if orig not in df.columns:
df[orig] = df[col].copy()
print(f" [filter] Backed up {col} -> {orig} ({len(df)} pts)")
[docs]
def lowpass_filter(
dfs: pd.DataFrame | list[pd.DataFrame],
columns: str | list[str],
cutoff_hz: float,
time_col: str = "time_s",
order: int = 4,
) -> pd.DataFrame | list[pd.DataFrame]:
"""Apply a Butterworth lowpass filter to one or more columns.
On first call for a given column the current data is saved to
``<col>_pre_lp``. Subsequent calls always filter from
``<col>_pre_lp`` so that re-filtering at a different cutoff does
not stack, while still preserving any upstream filtering (e.g.
``zscore_filter``) that was applied before this function.
A ``<col>_original`` backup is also created on the very first
filter call (zscore or lowpass) to allow full reset via
:func:`undo_filters`.
Parameters
----------
dfs : DataFrame or list[DataFrame]
Accepts a single DataFrame or a list; applies to all.
columns : str or list[str]
Column name(s) to filter.
cutoff_hz : float
Cutoff frequency in Hz.
time_col : str
Time column to derive sampling frequency from.
order : int
Butterworth filter order (applied forward-backward, so
effective order is 2x).
Returns
-------
Same type as *dfs*, with filtered columns overwritten and
``_pre_lp`` backups created on first filter.
"""
print(
"""
DEPRECATION WARNING: lowpass_filter is deprecated and will be removed in a future release.
Please reinstall the latest PERDA version. PERDA now features perda.utils.filtering.lowpass_filter, which
works natively with DataInstances.
"""
)
df_list, single = _ensure_list_df(dfs)
col_list = _check_blocked(_ensure_list_str(columns))
if not col_list:
return dfs
for df in df_list:
t = get_time_array(df, time_col)
dt_median = np.median(np.diff(t))
if dt_median <= 0:
raise ValueError(
"Non-positive median time step; cannot determine sample rate."
)
divisor = detect_time_divisor(np.asarray([t[0], t[-1]], dtype=np.float64))
fs = divisor / dt_median
nyq = fs / 2.0
if cutoff_hz >= nyq:
print(
f"[filter] Warning: cutoff ({cutoff_hz} Hz) >= Nyquist "
f"({nyq:.1f} Hz). Skipping filter."
)
return df_list[0] if single else df_list
sos = butter(order, cutoff_hz / nyq, btype="low", output="sos")
for col in col_list:
_save_original(df, col)
pre_lp = col + "_pre_lp"
if pre_lp not in df.columns:
df[pre_lp] = df[col].copy()
print(f" [filter] Backed up {col} -> {pre_lp} ({len(df)} pts)")
signal = df[pre_lp].values.astype(np.float64)
valid = ~np.isnan(signal)
if valid.sum() < 3 * (2 * order + 1):
print(
f" [filter] {col}: too few valid points ({valid.sum()}), skipping"
)
continue
filtered = np.full_like(signal, np.nan)
filtered[valid] = sosfiltfilt(sos, signal[valid])
df[col] = filtered
print(
f" [filter] {col}: lowpass @ {cutoff_hz} Hz ({len(df)} pts, fs={fs:.1f} Hz)"
)
return df_list[0] if single else df_list
[docs]
def lowpass_filter_by_distance(
dfs: pd.DataFrame | list[pd.DataFrame],
columns: str | list[str],
cutoff_freq_per_meter: float,
distance_col: str = "distance",
order: int = 4,
) -> pd.DataFrame | list[pd.DataFrame]:
"""Apply a Butterworth lowpass in the spatial domain (1/m).
Same semantics as :func:`lowpass_filter` but the sample rate is
derived from a distance column instead of time.
Parameters
----------
dfs : DataFrame or list[DataFrame]
columns : str or list[str]
cutoff_freq_per_meter : float
Cutoff in cycles per meter (1/m).
distance_col : str
order : int
Returns
-------
Same type as *dfs*.
"""
print(
"""
DEPRECATION WARNING: lowpass_filter_by_distance is deprecated and will be removed in a future release.
Please reinstall the latest PERDA version. PERDA now features perda.utils.filtering.lowpass_filter_by_distance,
which works natively with DataInstances.
"""
)
df_list, single = _ensure_list_df(dfs)
col_list = _check_blocked(_ensure_list_str(columns))
if not col_list:
return dfs
for df in df_list:
dist = df[distance_col].values.astype(np.float64)
dx_median = np.median(np.diff(dist))
if dx_median <= 0:
raise ValueError("Non-positive median distance step.")
fs = 1.0 / dx_median
nyq = fs / 2.0
if cutoff_freq_per_meter >= nyq:
print(
f"[filter] Warning: cutoff ({cutoff_freq_per_meter} 1/m) >= "
f"Nyquist ({nyq:.3f} 1/m). Skipping."
)
return df_list[0] if single else df_list
sos = butter(order, cutoff_freq_per_meter / nyq, btype="low", output="sos")
for col in col_list:
_save_original(df, col)
pre_lp = col + "_pre_lp"
if pre_lp not in df.columns:
df[pre_lp] = df[col].copy()
print(f" [filter] Backed up {col} -> {pre_lp} ({len(df)} pts)")
signal = df[pre_lp].values.astype(np.float64)
valid = ~np.isnan(signal)
if valid.sum() < 3 * (2 * order + 1):
print(f" [filter] {col}: too few valid points, skipping")
continue
filtered = np.full_like(signal, np.nan)
filtered[valid] = sosfiltfilt(sos, signal[valid])
df[col] = filtered
print(
f" [filter] {col}: lowpass @ {cutoff_freq_per_meter} 1/m "
f"({len(df)} pts, fs={fs:.3f} 1/m)"
)
return df_list[0] if single else df_list
# ---------------------------------------------------------------------------
# Rolling z-score outlier filter
# ---------------------------------------------------------------------------
[docs]
def zscore_filter(
dfs: pd.DataFrame | list[pd.DataFrame],
columns: str | list[str],
window_s: float = 2.0,
threshold: float = 4.8,
time_col: str = "time_s",
) -> pd.DataFrame | list[pd.DataFrame]:
"""Mask outliers using a rolling-window z-score.
For each sample, a z-score is computed relative to the local rolling
mean and standard deviation. Samples exceeding *threshold* are
replaced with ``NaN`` and linearly interpolated.
Always re-filters from ``<col>_original`` (the true intake data)
so that re-running with different parameters does not stack. If
a ``<col>_pre_lp`` backup exists (lowpass was already applied),
it is updated to the new zscore output so subsequent lowpass calls
chain correctly.
Parameters
----------
dfs : DataFrame or list[DataFrame]
columns : str or list[str]
Column name(s) to filter.
window_s : float
Rolling window size in seconds.
threshold : float
Z-score threshold -- samples with ``|z| > threshold`` are masked.
time_col : str
Time column used to convert *window_s* to a sample count.
Returns
-------
Same type as *dfs*, with outliers masked to NaN and
``_original`` backups created on first call.
"""
print(
"""
DEPRECATION WARNING: zscore_filter is deprecated and will be removed in a future release.
Please reinstall the latest PERDA version. PERDA now features perda.utils.filtering.zscore_filter, which
works natively with DataInstances.
"""
)
df_list, single = _ensure_list_df(dfs)
col_list = _check_blocked(_ensure_list_str(columns))
if not col_list:
return dfs
for df in df_list:
t = get_time_array(df, time_col)
dt_median = np.median(np.diff(t))
if dt_median <= 0:
raise ValueError("Non-positive median time step.")
divisor = detect_time_divisor(np.asarray([t[0], t[-1]], dtype=np.float64))
fs = divisor / dt_median
win_samples = max(3, int(round(window_s * fs)))
for col in col_list:
_save_original(df, col)
orig_col = col + "_original"
signal = df[orig_col].values.astype(np.float64)
series = pd.Series(signal)
roll_mean = series.rolling(win_samples, center=True, min_periods=1).mean()
roll_std = series.rolling(win_samples, center=True, min_periods=1).std()
roll_std = roll_std.replace(0, np.nan)
z = ((series - roll_mean) / roll_std).abs()
outliers = z > threshold
n_masked = int(outliers.sum())
filtered = signal.copy()
filtered[outliers.values] = np.nan
if n_masked > 0:
filtered = (
pd.Series(filtered)
.interpolate(method="linear", limit_direction="both")
.values
)
df[col] = filtered
# Keep _pre_lp in sync so subsequent lowpass re-runs
# chain on top of this zscore result.
pre_lp = col + "_pre_lp"
if pre_lp in df.columns:
df[pre_lp] = filtered.copy()
print(f" [zscore] Updated {pre_lp} to new zscore output")
print(
f" [zscore] {col}: window={window_s}s "
f"({win_samples} pts), threshold={threshold}, "
f"masked {n_masked}/{len(df)} pts (interpolated)"
)
return df_list[0] if single else df_list
# ---------------------------------------------------------------------------
# Undo filters
# ---------------------------------------------------------------------------
[docs]
def undo_filters(
df: pd.DataFrame,
columns: str | list[str],
) -> pd.DataFrame:
"""Restore columns to their original (pre-filter) state.
Copies data back from ``<col>_original`` and removes all backup
columns (``_original``, ``_pre_lp``, ``_raw``).
Parameters
----------
df : pd.DataFrame
columns : str or list[str]
Column name(s) to restore.
Returns
-------
pd.DataFrame
The same DataFrame, modified in-place and returned.
"""
print(
"""
DEPRECATION WARNING: undo_filters is deprecated and will be removed in a future release.
Please reinstall the latest PERDA version. The PERDA filtering functions do not modify the DataInstances
in place, so undo_filters is no longer necessary.
"""
)
col_list = _ensure_list_str(columns)
for col in col_list:
orig = col + "_original"
if orig not in df.columns:
print(f" [filter] {col}: no _original backup " f"found, skipping")
continue
df[col] = df[orig].copy()
for suffix in ("_original", "_pre_lp", "_raw"):
backup = col + suffix
if backup in df.columns:
df.drop(columns=backup, inplace=True)
print(f" [filter] {col}: restored from _original, " f"backups removed")
return df
# ---------------------------------------------------------------------------
# FFT spectrum plotting
# ---------------------------------------------------------------------------
[docs]
def plot_fft(
df: pd.DataFrame,
columns: str | list[str],
domain: str = "time",
time_col: str = "time_s",
distance_col: str = "distance",
stacked: bool = True,
cutoff: float | None = None,
log_x: bool = True,
log_y: bool = False,
xlim: tuple[float, float] | None = None,
) -> go.Figure:
"""Plot FFT magnitude spectrum of one or more columns.
Parameters
----------
df : pd.DataFrame
columns : str or list[str]
domain : str
``"time"`` or ``"distance"``.
time_col, distance_col : str
stacked : bool
``True`` -> one subplot per column; ``False`` -> all on one plot.
cutoff : float or None
If given, draw a vertical dashed line at this frequency.
log_x : bool
Logarithmic x-axis.
log_y : bool
Logarithmic y-axis.
xlim : (start, end) or None
Restrict the FFT to a subdomain of the x-axis (time in
seconds or distance in meters, depending on *domain*).
``None`` uses the full range.
Returns
-------
go.Figure
Note
----
No ``max_display_hz`` downsampling is applied here — the FFT
output is already in the frequency domain and thinning would
corrupt the spectrum shape.
"""
print(
"""
DEPRECATION WARNING: plot_fft is deprecated and will be removed in a future release.
Please reinstall the latest PERDA version. PERDA now features perda.plotting.fft_plotter.plot_fft_spectrum,
which replaces this functionality.
"""
)
col_list = _ensure_list_str(columns)
# Slice to subdomain if requested
if xlim is not None:
if domain == "distance":
x_arr = df[distance_col].values.astype(np.float64)
else:
x_arr = get_time_array(df, time_col)
mask = (x_arr >= xlim[0]) & (x_arr <= xlim[1])
df = df.loc[mask]
print(
f"[fft] Restricted to {domain} range "
f"[{xlim[0]}, {xlim[1]}]: {int(mask.sum())} pts"
)
# Determine sample spacing
if domain == "distance":
d = df[distance_col].values.astype(np.float64)
dx = np.median(np.diff(d))
freq_unit = "1/m"
else:
t = get_time_array(df, time_col)
dt = np.median(np.diff(t))
divisor = detect_time_divisor(np.asarray([t[0], t[-1]], dtype=np.float64))
dx = dt / divisor # sample spacing in seconds
freq_unit = "Hz"
if stacked:
n = len(col_list)
fig = make_subplots(
rows=n,
cols=1,
shared_xaxes=True,
subplot_titles=[f"FFT: {c}" for c in col_list],
vertical_spacing=0.04,
)
for i, col in enumerate(col_list, 1):
signal = df[col].dropna().values.astype(np.float64)
xf = rfftfreq(len(signal), d=dx)
yf = np.abs(rfft(signal - np.mean(signal)))
fig.add_trace(
go.Scattergl(x=xf, y=yf, mode="lines", name=col),
row=i,
col=1,
)
if cutoff is not None:
fig.add_vline(
x=cutoff,
line_dash="dash",
line_color="red",
annotation_text=f"{cutoff} {freq_unit}",
row=i,
col=1,
)
fig.update_yaxes(title_text="Magnitude", row=i, col=1)
fig.update_xaxes(title_text=f"Frequency ({freq_unit})", row=n, col=1)
if log_x:
fig.update_xaxes(type="log")
if log_y:
fig.update_yaxes(type="log")
fig.update_layout(height=300 * n, showlegend=False)
else:
fig = go.Figure()
for col in col_list:
signal = df[col].dropna().values.astype(np.float64)
xf = rfftfreq(len(signal), d=dx)
yf = np.abs(rfft(signal - np.mean(signal)))
fig.add_trace(go.Scattergl(x=xf, y=yf, mode="lines", name=col))
if cutoff is not None:
fig.add_vline(
x=cutoff,
line_dash="dash",
line_color="red",
annotation_text=f"{cutoff} {freq_unit}",
)
fig.update_layout(
xaxis_title=f"Frequency ({freq_unit})",
yaxis_title="Magnitude",
title="FFT Spectrum",
height=500,
)
if log_x:
fig.update_xaxes(type="log")
if log_y:
fig.update_yaxes(type="log")
return fig