whiten

function of dascore.proc.whiten source

whiten(
    patch: Patch ,
    smooth_size: float | None[float, None] = None,
    tukey_alpha: float = 0.1,
    idft: bool = True,
    **kwargs ,
)-> ‘PatchType’

Band-limited signal whitening. The whitened signal is returned in either frequency or time domains.

Parameters

Parameter Description
patch Patch to transform. Has to have dimensions of time and distance.
smooth_size Size (in Hz) of moving average window, used to compute the spectrum
before whitening. If no value is inputted, smoothing is over
the entire spectrum.
tukey_alpha Alpha parameter for Tukey window applied as windowing to the
smoothed spectrum within the required frequency range. By default
its value is 0.1.
See more details at https://docs.scipy.org/doc/scipy/reference
/generated/scipy.signal.windows.tukey.html
idft If False, returns the whitened result in the frequency domain without
converting it back to the time domain. Defaults to True.
**kwargs Used to specify the dimension and frequency, wavelength, or equivalent
limits. If no input is provided, whitening is also the last axis
with frequency band of [0,Nyquist]
Note

  1. The FFT result is divided by the smoothed spectrum before inverting back to time-domain signal. The phase is not changed.

  2. A tukey window (fixed) is applied to window the smoothed spectrum within the frequency range of interest. Be aware of its effect and consider enlarging the frequency range according to the tukey_alpha parameter.

  3. Amplitude is NOT preserved

  4. If idft = False, since for the purely real input data the negative frequency terms are just the complex conjugates of the corresponding positive-frequency terms, the output does not include the negative frequency terms, and therefore the length of the transformed axis of the output is n//2 + 1. Refer to the dft patch function and its real flag.

Example

import matplotlib.pyplot as plt
import numpy as np
import numpy.fft as fft

import dascore as dc
from dascore.units import Hz

def plot_spectrum(x, T, ax, phase=False):
    fftphase = np.angle(fft.fft(x))
    fftsig = np.abs(fft.fft(x))
    fftlen = fftsig.size
    fftsig = fftsig[0 : int(fftlen / 2) + 1]
    fftphase = fftphase[0 : int(fftlen / 2) + 1]
    freqvec = np.linspace(0, 0.5 / T, fftsig.size)
    if not phase:
        ax.plot(freqvec, fftsig)
        ax.set_xlabel("frequency [Hz]")
        ax.set_ylabel("Amplitude (|H(w)|)")
    else:
        ax.plot(freqvec, fftphase)
        ax.set_xlabel("frequency [Hz]")
        ax.set_ylabel("Phase (radians)")


patch = dc.get_example_patch("dispersion_event")
patch = patch.resample(time=(200 * Hz))

white_patch = patch.whiten(smooth_size=3, time = (10,50))

fig, ((ax1, ax2), (ax3, ax4), (ax5, ax6)) = plt.subplots(3, 2, figsize=(10, 7))
ax1.plot(patch.data[50, :])
ax1.set_title("Original data, distance = 50 m")
ax2.plot(white_patch.data[50, :])
ax2.set_title("Whitened data, distance = 50 m")

plot_spectrum(patch.data[50, :], 1 / 200, ax3)
ax3.set_title("Original data, distance = 50 m")
plot_spectrum(white_patch.data[50, :], 1 / 200, ax4)
ax4.set_title("Whitened data, distance = 50 m")

plot_spectrum(patch.data[50, :], 1 / 200, ax5, phase=True)
ax5.set_title("Original data, distance = 50 m")
plot_spectrum(white_patch.data[50, :], 1 / 200, ax6, phase=True)
ax6.set_title("Whitened data, distance = 50 m")
plt.tight_layout()
plt.show()

```