deploy: b27d7a2

.buildinfo

@@ -1,4 +1,4 @@
 # Sphinx build info version 1
 # This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done.
-config: a6ac25651f02e78abf3a4ac84ba4923e
+config: 713ec541478fa44eebb19d03a9a7ca96
 tags: 645f666f9bcd5a90fca523b33c5a78b7

_downloads/07f4b1b69497c93530bbe772d7eb535a/smooth.py

@@ -12,11 +12,6 @@
 from scipy.signal import welch
 
 from movement import sample_data
-from movement.filtering import (
-    interpolate_over_time,
-    median_filter,
-    savgol_filter,
-)
 
 # %%
 # Load a sample dataset
@@ -30,9 +25,10 @@
 print(ds_wasp)
 
 # %%
-# We see that the dataset contains a single individual (a wasp) with two
-# keypoints tracked in 2D space. The video was recorded at 40 fps and lasts for
-# ~27 seconds.
+# We see that the dataset contains the 2D pose tracks and confidence scores
+# for a single wasp, generated with DeepLabCut. The wasp is tracked at two
+# keypoints: "head" and "stinger" in a video that was recorded at 40 fps and
+# lasts for approximately 27 seconds.
 
 # %%
 # Define a plotting function
@@ -67,7 +63,7 @@ def plot_raw_and_smooth_timeseries_and_psd(
     fig, ax = plt.subplots(2, 1, figsize=(10, 6))
 
     for ds, color, label in zip(
-        [ds_raw, ds_smooth], ["k", "r"], ["raw", "smooth"]
+        [ds_raw, ds_smooth], ["k", "r"], ["raw", "smooth"], strict=False
     ):
         # plot position time series
         pos = ds.position.sel(**selection)
@@ -80,9 +76,10 @@ def plot_raw_and_smooth_timeseries_and_psd(
             label=f"{label} {space}",
         )
 
-        # generate interpolated dataset to avoid NaNs in the PSD calculation
-        ds_interp = interpolate_over_time(ds, max_gap=None, print_report=False)
-        pos_interp = ds_interp.position.sel(**selection)
+        # interpolate data to remove NaNs in the PSD calculation
+        pos_interp = ds.sel(**selection).move.interpolate_over_time(
+            print_report=False
+        )
         # compute and plot the PSD
         freq, psd = welch(pos_interp, fs=ds.fps, nperseg=256)
         ax[1].semilogy(
@@ -111,20 +108,44 @@ def plot_raw_and_smooth_timeseries_and_psd(
 # %%
 # Smoothing with a median filter
 # ------------------------------
-# Here we use the :py:func:`movement.filtering.median_filter` function to
-# apply a rolling window median filter to the wasp dataset.
-# The ``window_length`` parameter is defined in seconds (according to the
-# ``time_unit`` dataset attribute).
+# Using the
+# :py:meth:`median_filter()\
+# <movement.move_accessor.MovementDataset.filtering_wrapper>`
+# method of the ``move`` accessor,
+# we apply a rolling window median filter over a 0.1-second window
+# (4 frames) to the wasp dataset.
+# As the ``window`` parameter is defined in *number of observations*,
+# we can simply multiply the desired time window by the frame rate
+# of the video. We will also create a copy of the dataset to avoid
+# modifying the original data.
+
+window = int(0.1 * ds_wasp.fps)
+ds_wasp_smooth = ds_wasp.copy()
+ds_wasp_smooth.update({"position": ds_wasp_smooth.move.median_filter(window)})
 
-ds_wasp_medfilt = median_filter(ds_wasp, window_length=0.1)
+# %%
+# .. note::
+#    The ``move`` accessor :py:meth:`median_filter()\
+#    <movement.move_accessor.MovementDataset.filtering_wrapper>`
+#    method is a convenience method that applies
+#    :py:func:`movement.filtering.median_filter`
+#    to the ``position`` data variable.
+#    The equivalent function call using the
+#    :py:mod:`movement.filtering` module would be:
+#
+#    .. code-block:: python
+#
+#       from movement.filtering import median_filter
+#
+#       ds_wasp_smooth.update({"position": median_filter(position, window)})
 
 # %%
 # We see from the printed report that the dataset has no missing values
 # neither before nor after smoothing. Let's visualise the effects of the
 # median filter in the time and frequency domains.
 
 plot_raw_and_smooth_timeseries_and_psd(
-    ds_wasp, ds_wasp_medfilt, keypoint="stinger"
+    ds_wasp, ds_wasp_smooth, keypoint="stinger"
 )
 
 # %%
@@ -142,8 +163,8 @@ def plot_raw_and_smooth_timeseries_and_psd(
 # %%
 # Choosing parameters for the median filter
 # ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-# You can control the behaviour of :py:func:`movement.filtering.median_filter`
-# via two parameters: ``window_length`` and ``min_periods``.
+# We can control the behaviour of the median filter
+# via two parameters: ``window`` and ``min_periods``.
 # To better understand the effect of these parameters, let's use a
 # dataset that contains missing values.
 
@@ -154,10 +175,15 @@ def plot_raw_and_smooth_timeseries_and_psd(
 # The dataset contains a single mouse with six keypoints tracked in
 # 2D space. The video was recorded at 30 fps and lasts for ~616 seconds. We can
 # see that there are some missing values, indicated as "nan" in the
-# printed dataset. Let's apply the median filter to this dataset, with
-# the ``window_length`` set to 0.1 seconds.
-
-ds_mouse_medfilt = median_filter(ds_mouse, window_length=0.1)
+# printed dataset.
+# Let's apply the median filter over a 0.1-second window (3 frames)
+# to the dataset.
+
+window = int(0.1 * ds_mouse.fps)
+ds_mouse_smooth = ds_mouse.copy()
+ds_mouse_smooth.update(
+    {"position": ds_mouse_smooth.move.median_filter(window)}
+)
 
 # %%
 # The report informs us that the raw data contains NaN values, most of which
@@ -172,7 +198,9 @@ def plot_raw_and_smooth_timeseries_and_psd(
 # For example, setting ``min_periods=2`` means that two non-NaN values in the
 # window are sufficient for the median to be calculated. Let's try this.
 
-ds_mouse_medfilt = median_filter(ds_mouse, window_length=0.1, min_periods=2)
+ds_mouse_smooth.update(
+    {"position": ds_mouse.move.median_filter(window, min_periods=2)}
+)
 
 # %%
 # We see that this time the number of NaN values has decreased
@@ -183,51 +211,64 @@ def plot_raw_and_smooth_timeseries_and_psd(
 # parts of the data.
 
 plot_raw_and_smooth_timeseries_and_psd(
-    ds_mouse, ds_mouse_medfilt, keypoint="snout", time_range=slice(0, 80)
+    ds_mouse, ds_mouse_smooth, keypoint="snout", time_range=slice(0, 80)
 )
 
 # %%
 # The smoothing once again reduces the power of high-frequency components, but
 # the resulting time series stays quite close to the raw data.
 #
-# What happens if we increase the ``window_length`` to 2 seconds?
+# What happens if we increase the ``window`` to 2 seconds (60 frames)?
 
-ds_mouse_medfilt = median_filter(ds_mouse, window_length=2, min_periods=2)
+window = int(2 * ds_mouse.fps)
+ds_mouse_smooth.update(
+    {"position": ds_mouse.move.median_filter(window, min_periods=2)}
+)
 
 # %%
-# The number of NaN values has decreased even further. That's because the
-# chance of finding at least 2 valid values within a 2 second window is
-# quite high. Let's plot the results for the same keypoint and time range
+# The number of NaN values has decreased even further.
+# That's because the chance of finding at least 2 valid values within
+# a 2-second window (i.e. 60 frames) is quite high.
+# Let's plot the results for the same keypoint and time range
 # as before.
 
 plot_raw_and_smooth_timeseries_and_psd(
-    ds_mouse, ds_mouse_medfilt, keypoint="snout", time_range=slice(0, 80)
+    ds_mouse, ds_mouse_smooth, keypoint="snout", time_range=slice(0, 80)
 )
 # %%
 # We see that the filtered time series is much smoother and it has even
 # "bridged" over some small gaps. That said, it often deviates from the raw
 # data, in ways that may not be desirable, depending on the application.
-# That means that our choice of ``window_length`` may be too large.
-# In general, you should choose a ``window_length`` that is small enough to
+# Here, our choice of ``window`` may be too large.
+# In general, you should choose a ``window`` that is small enough to
 # preserve the original data structure, but large enough to remove
 # "spikes" and high-frequency noise. Always inspect the results to ensure
 # that the filter is not removing important features.
 
 # %%
 # Smoothing with a Savitzky-Golay filter
 # --------------------------------------
-# Here we use the :py:func:`movement.filtering.savgol_filter` function,
-# which is a wrapper around :py:func:`scipy.signal.savgol_filter`.
+# Here we use the
+# :py:meth:`savgol_filter()\
+# <movement.move_accessor.MovementDataset.filtering_wrapper>`
+# method of the ``move`` accessor, which is a convenience method that applies
+# :py:func:`movement.filtering.savgol_filter`
+# (a wrapper around :py:func:`scipy.signal.savgol_filter`),
+# to the ``position`` data variable.
 # The Savitzky-Golay filter is a polynomial smoothing filter that can be
-# applied to time series data on a rolling window basis. A polynomial of
-# degree ``polyorder`` is fitted to the data in each window of length
-# ``window_length``, and the value of the polynomial at the center of the
-# window is used as the output value.
+# applied to time series data on a rolling window basis.
+# A polynomial with a degree specified by ``polyorder`` is applied to each
+# data segment defined by the size ``window``.
+# The value of the polynomial at the midpoint of each ``window`` is then
+# used as the output value.
 #
-# Let's try it on the mouse dataset.
-
-ds_mouse_savgol = savgol_filter(ds_mouse, window_length=0.2, polyorder=2)
+# Let's try it on the mouse dataset, this time using a 0.2-second
+# window (i.e. 6 frames) and the default ``polyorder=2`` for smoothing.
+# As before, we first compute the corresponding number of observations
+# to be used as the ``window`` size.
 
+window = int(0.2 * ds_mouse.fps)
+ds_mouse_smooth.update({"position": ds_mouse.move.savgol_filter(window)})
 
 # %%
 # We see that the number of NaN values has increased after filtering. This is
@@ -238,64 +279,71 @@ def plot_raw_and_smooth_timeseries_and_psd(
 # effects in the time and frequency domains.
 
 plot_raw_and_smooth_timeseries_and_psd(
-    ds_mouse, ds_mouse_savgol, keypoint="snout", time_range=slice(0, 80)
+    ds_mouse, ds_mouse_smooth, keypoint="snout", time_range=slice(0, 80)
 )
 # %%
 # Once again, the power of high-frequency components has been reduced, but more
 # missing values have been introduced.
 
 # %%
-# Now let's take a look at the wasp dataset.
+# Now let's apply the same Savitzky-Golay filter to the wasp dataset.
 
-ds_wasp_savgol = savgol_filter(ds_wasp, window_length=0.2, polyorder=2)
+window = int(0.2 * ds_wasp.fps)
+ds_wasp_smooth.update({"position": ds_wasp.move.savgol_filter(window)})
 
 # %%
 plot_raw_and_smooth_timeseries_and_psd(
-    ds_wasp,
-    ds_wasp_savgol,
-    keypoint="stinger",
+    ds_wasp, ds_wasp_smooth, keypoint="stinger"
 )
 # %%
 # This example shows two important limitations of the Savitzky-Golay filter.
 # First, the filter can introduce artefacts around sharp boundaries. For
 # example, focus on what happens around the sudden drop in position
 # during the final second. Second, the PSD appears to have large periodic
 # drops at certain frequencies. Both of these effects vary with the
-# choice of ``window_length`` and ``polyorder``. You can read more about these
+# choice of ``window`` and ``polyorder``. You can read more about these
 # and other limitations of the Savitzky-Golay filter in
 # `this paper <https://pubs.acs.org/doi/10.1021/acsmeasuresciau.1c00054>`_.
 
 
 # %%
 # Combining multiple smoothing filters
 # ------------------------------------
-# You can also combine multiple smoothing filters by applying them
+# We can also combine multiple smoothing filters by applying them
 # sequentially. For example, we can first apply the median filter with a small
-# ``window_length`` to remove "spikes" and then apply the Savitzky-Golay filter
-# with a larger ``window_length`` to further smooth the data.
+# ``window`` to remove "spikes" and then apply the Savitzky-Golay filter
+# with a larger ``window`` to further smooth the data.
 # Between the two filters, we can interpolate over small gaps to avoid the
 # excessive proliferation of NaN values. Let's try this on the mouse dataset.
-# First, let's apply the median filter.
+# First, we will apply the median filter.
 
-ds_mouse_medfilt = median_filter(ds_mouse, window_length=0.1, min_periods=2)
+window = int(0.1 * ds_mouse.fps)
+ds_mouse_smooth.update(
+    {"position": ds_mouse.move.median_filter(window, min_periods=2)}
+)
 
 # %%
-# Next, let's linearly interpolate over gaps smaller than 1 second.
+# Next, let's linearly interpolate over gaps smaller than 1 second (30 frames).
 
-ds_mouse_medfilt_interp = interpolate_over_time(ds_mouse_medfilt, max_gap=1)
+ds_mouse_smooth.update(
+    {"position": ds_mouse_smooth.move.interpolate_over_time(max_gap=30)}
+)
 
 # %%
-# Finally, let's apply the Savitzky-Golay filter.
+# Finally, let's apply the Savitzky-Golay filter over a 0.4-second window
+# (12 frames).
 
-ds_mouse_medfilt_interp_savgol = savgol_filter(
-    ds_mouse_medfilt_interp, window_length=0.4, polyorder=2
+window = int(0.4 * ds_mouse.fps)
+ds_mouse_smooth.update(
+    {"position": ds_mouse_smooth.move.savgol_filter(window)}
 )
 
 # %%
-# A record of all applied operations is stored in the dataset's ``log``
-# attribute. Let's inspect it to summarise what we've done.
+# A record of all applied operations is stored in the ``log`` attribute of the
+# ``ds_mouse_smooth.position`` data array. Let's inspect it to summarise
+# what we've done.
 
-for entry in ds_mouse_medfilt_interp_savgol.log:
+for entry in ds_mouse_smooth.position.log:
     print(entry)
 
 # %%
@@ -304,11 +352,17 @@ def plot_raw_and_smooth_timeseries_and_psd(
 
 plot_raw_and_smooth_timeseries_and_psd(
     ds_mouse,
-    ds_mouse_medfilt_interp_savgol,
+    ds_mouse_smooth,
     keypoint="snout",
     time_range=slice(0, 80),
 )
 
 # %%
 # Feel free to play around with the parameters of the applied filters and to
 # also look at other keypoints and time ranges.
+
+# %%
+# .. seealso::
+#   :ref:`examples/filter_and_interpolate:Filtering multiple data variables`
+#   in the
+#   :ref:`sphx_glr_examples_filter_and_interpolate.py` example.