DataAnalysis2021/05-Spectrogram/nonlinear_thresholding.ipynb

Nonlinear Thresholding

When noise and signal share a common frequency band, spectral filtering would lead to a loss in signal. For these cases, we need different techniques to reduce the disturibing noise. We start with a reocorde signal x(t), contains some signal s(t) and additive noise n(t):

x(t) = s(t) + n(t)

We can rewrite the equation above in time-frequency domain as

X(t, f) = S(t, f) + N(t, f) ~.

Assuming that some time-frequnecy coefficients can be associated by noise, we set all coefficients, which a below this threshold, to zero:

\tilde{X}(t, f) = \left\{
    \begin{array}{@{}ll@{}}
         X(t, f) & \mathrm{if}~ |X(t, f)| \geq \beta(f) \\
          0          & \mathrm{otherwise}
    \end{array}\right.
    ~,

where \beta(f) denotes the threshold function. The threshold function is defined as

\beta(f) = \mathrm{ECDF}_f^{-1} (P = 0.99) ~,

where \mathrm{ECDF}_f^{-1} denotes the inverse cumulative distribution function or quantile function, e.g. before the first arrival of earthquake waves.

In [ ]:

# Load all required packages
import os
import numpy as np
from scipy.signal import stft, istft
import matplotlib.pyplot as plt

In [ ]:

# Load our earthquake data and plot
d = np.load(os.path.join(os.path.expanduser('~'),'work', 'data', 'events', 'bug2019mgoh_Z.npz'))
#d = np.load(os.path.join(os.path.expanduser('~'),'work', 'data', 'events', 'bug2019gbbo_Z.npz'))
#d = np.load(os.path.join(os.path.expanduser('~'),'work', 'data', 'events', 'bug2019ibsd_Z.npz'))
plt.figure(figsize=(15, 8))
plt.plot(d['data']);

In [ ]:

# Compute spectrogram of data
f, t, X = stft(d['data'], fs=1/0.01, nfft=198, nperseg=99)   # Be careful with choice of nfft and nperseg, 
                                                             # because istft does not work for all pairs
print(X.shape)
plt.figure(figsize=(15, 8))
plt.pcolormesh(t, f, np.abs(X))
plt.xlabel("Frequency (Hz)")
plt.ylabel("Time (s)");

In [ ]:

# Functions for thresholding
def threshold(X, quantile=0.99):
    """
    X: time-frequency coefficients
    q: quantile [0-1]
    """
    # Loop over frequencies to build threshold function
    beta = np.zeros(X.shape[0])
    for i in range(X.shape[0]):
        beta[i] = np.quantile(np.abs(X[i, :]), q=quantile)
    
    return beta

def modify_spectrogram(X, beta):
    # Loop over all items in X and apply threshold
    # Task: modify X!
    
    return X

In [ ]:

# Compute STFT of data before first arrival of P
_, _, X_thres = stft(d['data'][:2500], fs=1/0.01, nfft=198, nperseg=99)
# Estimate threshold fucntion
beta = threshold(X_thres)

# Modifiy original spectrogram
X_mod = modify_spectrogram(X, beta)

# Plot modified spectrogram
plt.figure(figsize=(15, 8))
plt.pcolormesh(t, f, np.abs(X_mod))
plt.xlabel("Frequency (Hz)")
plt.ylabel("Time (s)");

In [ ]:

# Inverse STFT of modified spectrogram
t, x_mod = istft(X_mod, fs=1/0.01, nfft=198, nperseg=99)

# Plot corrected seismogram
plt.figure(figsize=(15, 8))
plt.plot(t, x_mod);

In [ ]: