Skip to content

Commit

Permalink
add docs, fix test coverage
Browse files Browse the repository at this point in the history
  • Loading branch information
d-chambers committed Oct 31, 2024
1 parent 301cd42 commit b7c2f39
Show file tree
Hide file tree
Showing 3 changed files with 147 additions and 3 deletions.
24 changes: 21 additions & 3 deletions dascore/proc/whiten.py
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,11 @@ def _get_dim_freq_range_from_kwargs(patch, kwargs):
msg = f"passed dim of {dim} to whiten but it is not in patch dimensions."
raise ParameterError(msg)
else: # Something when wrong.
msg = "Whiten kwargs must specify a single patch dimension."
msg = (
"Whiten kwargs must specify a single patch dimension. "
"Double check the supported input parameters as these may have "
f"changed. You passed {kwargs}."
)
raise ParameterError(msg)

return dim, freq_range
Expand Down Expand Up @@ -85,7 +89,8 @@ def whiten(
Spectral whitening of a signal.
The whitened signal is returned in the same domain (eq frequency or
time domain) as the input signal.
time domain) as the input signal. See also the
[Whiten Processing Section](`docs/tutorial/processing.qmd`#whiten).
Parameters
----------
Expand Down Expand Up @@ -116,8 +121,21 @@ def whiten(
Example
-------
>>> import dascore as dc
>>>
>>> patch = dc.get_example_patch()
>>>
>>> # Whiten along time dimension
>>> white_patch = patch.whiten(time=None)
>>>
>>> # Band limited whitening
>>> white_patch = patch.whiten(time=(20, 40))
>>>
>>> # Band limited with taper ends
>>> white_patch = patch.whiten(time=(10, 20, 40, 60))
>>>
>>> # Whitening along distance with amplitude smoothing (0.1/m))
>>> white_patch = patch.whiten(smooth_size=0.1, distance=None)
"""
dim, freq_range = _get_dim_freq_range_from_kwargs(patch, kwargs)
fft_dim = FourierTransformatter().rename_dims(dim)[0]
Expand Down
112 changes: 112 additions & 0 deletions docs/tutorial/processing.qmd
Original file line number Diff line number Diff line change
Expand Up @@ -169,3 +169,115 @@ smoothed.viz.waterfall(scale=.1);
```

Notice the nan values at the start of the time axis. These can be trimmed with [`Patch`.dropna](`dascore.Patch.dropna`).


# Whiten

The [`Patch.whiten`](`dascore.Patch.whiten`) function performs spectral whitening by balancing the amplitude spectra of the patch while leaving the phase (largely) unchanged. Spectral whitening is often a pre-processing step in ambient noise correlation workflows.

To demonstrate, we create some plotting code and an example patch.

```{python}
import matplotlib.pyplot as plt
import numpy as np
import dascore as dc
from dascore.utils.time import to_float
rng = np.random.default_rng()
def plot_time_and_frequency(patch, channel=0):
""" Make plots of time and frequency of patch with single channel."""
sq_patch = patch.select(distance=channel, samples=True).squeeze()
time_array = to_float(patch.get_array("time"))
time = time_array - np.min(time_array)
fig, (td_ax, fd_ax, phase_ax) = plt.subplots(1, 3, figsize=(9, 2.5))
# Plot in time domain
td_ax.plot(time, sq_patch.data, color="tab:blue")
td_ax.set_title("Time Domain")
td_ax.set_xlabel("time (s)")
# Plot freq amplitdue
ft_patch = sq_patch.dft("time", real=True)
freq = ft_patch.get_array("ft_time")
fd_ax.plot(freq, ft_patch.abs().data, color="tab:red")
fd_ax.set_xlabel("Frequency (Hz)")
fd_ax.set_title("Amplitude Spectra")
# plot freq phase
phase_ax.plot(freq, np.angle(ft_patch.data), color="tab:cyan")
phase_ax.set_xlabel("Frequency (Hz)")
phase_ax.set_title("Phase Angle")
# fd_ax.set_xlim(0, 1000)
return fig
def make_noisy_sine_patch():
"""Make a noisy sine wave patch."""
patch = dc.get_example_patch(
"sin_wav",
frequency=[1, 10, 20, 54, 66],
amplitude=[1, 2, 3, 4, 9],
duration=1,
sample_rate=200,
)
rand_noise = (rng.random(patch.data.shape) - 0.5) * 10
patch = patch.new(data = patch.data + rand_noise)
return patch
```

```{python}
patch = make_noisy_sine_patch()
plot_time_and_frequency(patch);
```

The default whitening makes all spectral amplitudes equal.

```{python}
white_patch = patch.whiten()
plot_time_and_frequency(white_patch);
```

The whitening can be restricted to certain frequency bands by specifying the dimension and a frequency range.

```{python}
white_patch = patch.whiten(time=(20, 80))
plot_time_and_frequency(white_patch);
```

Four values can be used to control the start/end of the taper.

```{python}
white_patch = patch.whiten(time=(20, 40, 60, 80))
plot_time_and_frequency(white_patch);
```

Whiten also supports using a smoothed amplitude for normalization which causes less drastic changes to the amplitude spectrum.

```{python}
white_patch = patch.whiten(smooth_size=1)
_ = plot_time_and_frequency(white_patch)
```


```{python}
white_patch = patch.whiten(smooth_size=2, time=(10, 20, 80, 90))
_ = plot_time_and_frequency(white_patch)
```

`Whiten` also accepts patches in the frequency domain, in which case a frequency patches are returned. This might be useful if whiten is only a part of a frequency domain workflow.

```{python}
fft_patch = patch.dft("time", real=True)
dft_white = fft_patch.whiten(smooth_size=1, time=(20, 40, 60, 80))
td_patch = dft_white.idft()
plot_time_and_frequency(td_patch);
```

The `water_level` parameter can be useful for stabilizing frequencies that may have near-zero values.
14 changes: 14 additions & 0 deletions tests/test_proc/test_whiten.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@
import numpy as np
import pytest

import dascore as dc
from dascore import get_example_patch
from dascore.exceptions import ParameterError
from dascore.units import Hz
Expand Down Expand Up @@ -204,6 +205,15 @@ def test_whiten_dft_input(self, test_patch):

assert "ft_time" in dft.coords.coord_map

def test_whiten_df_all_parameters(self, test_patch):
"""Ensure whiten accepts all args in dft form."""
fft_patch = test_patch.dft("time", real=True)

whitened_patch_freq_domain = fft_patch.whiten(
smooth_size=5, time=(30, 60), water_level=.01
)
assert isinstance(whitened_patch_freq_domain, dc.Patch)

def test_no_time_no_kwargs_raises(self, random_patch):
"""Ensure if no kwargs and patch doesn't have time an error is raised."""
patch = random_patch.rename_coords(time="money")
Expand All @@ -222,3 +232,7 @@ def test_multiple_kwargs_raises(self, random_patch):
msg = "must specify a single patch dimension"
with pytest.raises(ParameterError, match=msg):
random_patch.whiten(time=None, distance=None)

def test_helpful_error_message(self, random_patch):
"""Ensure helpful error message is used when a bad kwarg is passed."""

0 comments on commit b7c2f39

Please sign in to comment.