在 Python 偵測極值

摘要

偵測訊號的極值。

訊號的產生

舉例產生一個偽造的訊號。

• 3 [Hz] 的訊號 + 0.01 [Hz] 的訊號 + 雜訊

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from scipy import signal

t = np.arange(30 * 100) / 100
x = np.sin(2 * np.pi * t / 3) + np.random.randn(len(t)) * 0.2 + np.sin(2 * np.pi * t / 100)
df = pd.DataFrame({'Data': x}, index=pd.to_datetime(t, unit='s').time)

fig, ax = plt.subplots()
df.plot(ax=ax, xlim=['00:00:00', '00:00:10'])
plt.show()

檢查訊號的特性

def acorr(df: pd.DataFrame, ra: int = 3, fs=100):
    x = np.correlate(df.Data.values, df.Data.values, mode='full')
    t = (np.arange(len(x)) - len(x) / 2 + 0.5) / fs
    x /= x.max()
    fig, ax = plt.subplots()
    ax.set_xlim(-ra, ra)
    ax.set_ylim(0, 1)
    ax.plot(t, x)
    ax.set_title('Autocorrelation')
    ax.grid(axis='x')
    ax.set_xticks(np.linspace(-ra, ra, 2 * ra + 1))
    plt.show()

acorr(df)

可以看出主要的訊號是 1/3 [Hz]。因此，對 1/3 [Hz] 施加低通濾波器。

套用低通濾波器並偵測極值

套用低通濾波器以便偵測極值。

如果在 1/3 [Hz] 施加低通濾波器仍無法正確偵測，請逐步降低頻率。

fs = 10
lpf = 0.1
b, a = signal.butter(5, lpf / fs * 2, 'low')
df['Filter'] = signal.filtfilt(b, a, df.Data.values)
max_idx = signal.argrelmax(df['Filter'].values)
min_idx = signal.argrelmin(df['Filter'].values)

df['max_index'] = False
df.iloc[max_idx[0], 2] = True
df['min_index'] = False
df.iloc[min_idx[0], 3] = True

fig, ax = plt.subplots()
df.plot(ax=ax)
ax.set_xlim(['00:00:00', '00:00:30'])
ax.scatter(
    df.loc[df['max_index'], ['Filter']].index, 
    df.loc[df['max_index'], ['Filter']], 
    color='tab:orange', 
    zorder=3)
ax.scatter(
    df.loc[df['min_index'], ['Filter']].index, 
    df.loc[df['min_index'], ['Filter']], 
    color='tab:green', 
    zorder=3)
plt.show()