Envelope for Least Square Filtering and Smoothing

Signal Envelope

scipy.signal.hilbert
scipy.interpolate.interp1d

  • 복잡한 신호의 전반적인 형상을 파악하기 위해 envelope을 추출합니다.
  • python에서는 scipy.signal.hilbert를 통해 analytic signal을 추출하는 방법을 씁니다.

    1
    2
    x = np.arange(0, 20.1, 0.1)
    y = abs(np.sin(x)) * np.sin(x*20)


  • 이와 같은 신호에 hilbert를 다음과 같이 적용하면 다음과 같이 숨겨진 sine wave를 찾아냅니다.

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    from scipy.signal import hilbert, chirp

    duration = 1.0
    fs = 400.0
    samples = int(fs*duration)
    t = np.arange(samples) / fs

    analytic_signal = hilbert(y)
    amplitude_envelope = np.abs(analytic_signal)
    instantaneous_phase = np.unwrap(np.angle(analytic_signal))
    instantaneous_frequency = np.diff(instantaneous_phase) / (2.0*np.pi) * fs

    fig, ax = plt.subplots(nrows=2, sharex=True)
    ax[0].plot(x, y, c="b")
    ax[0].plot(x, amplitude_envelope)
    ax[1].plot(x[1:], instantaneous_frequency)
    ax[1].set_xlabel("time in seconds")


  • 그러나, 위 신호 그대로의 모습을 추출하려면 이론식을 찾는 hilbert는 답답한 면이 있습니다.

  • 그럴 때는 아래와 같이 그래프의 요철을 감지해서 interpolation 하는 방법이 직관적입니다.

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    def get_envelope_v1(x, y):
    x_list, y_list = list(x), list(y)
    assert len(x_list) == len(y_list)

    # First data
    ui, ux, uy = [0], [x_list[0]], [y_list[0]]
    li, lx, ly = [0], [x_list[0]], [y_list[0]]

    # Find upper peaks and lower peaks
    for i in range(1, len(x_list)-1):
    if y_list[i] >= y_list[i-1] and y_list[i] >= y_list[i+1]:
    ui.append(i)
    ux.append(x_list[i])
    uy.append(y_list[i])
    if y_list[i] <= y_list[i-1] and y_list[i] <= y_list[i+1]:
    li.append(i)
    lx.append(x_list[i])
    ly.append(y_list[i])

    # Last data
    ui.append(len(x_list)-1)
    ux.append(x_list[-1])
    uy.append(y_list[-1])
    li.append(len(y_list)-1)
    lx.append(x_list[-1])
    ly.append(y_list[-1])

    if len(ux) == 2 or len(lx) == 2:
    return [], []

    else:
    func_ub = interp1d(ux, uy, kind='cubic', bounds_error=False)
    func_lb = interp1d(lx, ly, kind='cubic', bounds_error=False)

    ub, lb = [], []
    for i in x_list:
    ub = func_ub(x_list)
    lb = func_lb(x_list)

    ub = np.array([y, ub]).max(axis=0)
    lb = np.array([y, lb]).min(axis=0)

    return ub, lb
    1
    2
    3
    4
    5
    6
    7
    8
    uy, ly = get_envelope_v1(x, y)

    fig, ax = plt.subplots()
    ax.plot(x, y, c="b")
    ax.plot(x, uy, c="m")
    ax.plot(x, ly, c="c")
    ax.set_xlim(0, 20)
    plt.show()


  • signal envelope은 보통 전체적인 형상을 추출하는데 사용되지만, smoothing 제어에 사용될 수도 있습니다.

Least Squares Smoothing

scipy.signal.savgol_filter
Least Squares Filtering and Smoothing, including Savitzky-Golay
Plateau Detection

  • 신호와 함께 존재하기 마련인 잡음을 제거하기 위해 smoothing을 적용합니다.
  • savinsky-golay filter와 같은 least squares filter는 mean이나 median에 비해 신호를 보존하고 잡음을 선택적으로 잘 제거한다는 장점이 있습니다.
  • 그러나 치명적인 단점이 있는데, signal이 크게 바뀔 경우 overshooting이 발생한다는 것입니다.


  • envelope을 이용해 overshooting의 범위를 제한해 보겠습니다.

  • 그러나 데이터가 이렇게 생기면 envelope도 깔끔하진 않습니다.
  • 같은 코드를 적용해도, 아래와 같이 envelope의 undershooting이 커서 문제가 됩니다.


  • envelop과 원래 데이터를 다시 비교해서 envelope의 범위를 제어합니다.

  • get_envelope()에 아래 부분을 추가합니다.

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    # upper through check
    ux_chk, uy_chk = [x_list[0]], [y_list[0]]
    for ii in range(1, len(ui)):
    slope = (uy[ii]-uy[ii-1])/(ux[ii]-ux[ii-1])
    intercept = -slope*ux[ii] + uy[ii]

    for xi in range(ui[ii-1], ui[ii]):
    y_pred = slope*x_list[xi] + intercept
    if y_pred < y_list[xi] and (x_list[xi] not in ux_chk) and (y_list[xi] not in uy_chk):
    ux_chk.append(x_list[xi])
    uy_chk.append(y_list[xi])
    ux_chk.append(x_list[ui[ii]])
    uy_chk.append(y_list[ui[ii]])

    # lower through check
    lx_chk, ly_chk = [x_list[0]], [y_list[0]]
    for ii in range(1, len(li)):
    slope = (ly[ii]-ly[ii-1])/(lx[ii]-lx[ii-1])
    intercept = -slope*lx[ii] + ly[ii]

    for xi in range(li[ii-1], li[ii]):
    y_pred = slope*x_list[xi] + intercept
    if y_pred > y_list[xi] and (x_list[xi] not in lx_chk) and (y_list[xi] not in ly_chk):
    lx_chk.append(x_list[xi])
    ly_chk.append(y_list[xi])
    lx_chk.append(x_list[li[ii]])
    ly_chk.append(y_list[li[ii]])


  • envelope의 범위가 제어되었습니다.

  • 이제 envelope을 벗어난 signal 대신 envelope을 적용합니다.

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    def savgol(x, y, box_size, recur=1):
    ub, lb = get_envelope(x, y)

    for _ in range(recur):
    y = sg(y, box_size, 2, mode='nearest')
    if len(ub) > 0:
    y = np.where(y > ub, ub, np.where(y < lb, lb, y))

    return y

    y_savgol = savgol(X, Y, 15)


  • Savgol filter의 overshooting 문제가 해결되었습니다.

  • 강제로 집어넣은 부분이 다른 부분에 비해 다소 날카롭습니다.
  • 그러나 overshooting에 비해서는 완화된 것을 확인할 수 있습니다.
  • 전체 코드는 여기에서 확인 가능합니다.


도움이 되셨나요?

카페인을 투입하시면 다음 포스팅으로 변환됩니다

Share