HP 滤波 (Hodrick Prescott Filter)

目的

HP 滤波是一种分解信号的低频趋势与高频部分的滤波方法：
$$y_t=g_t+c_t$$

其中 $y$ 代表原始信号，$g_t$ 代表信号的低频趋势，$c_t$为更高频的周期或噪声。

优化函数

$$\underset{g}{\min }\sum _i^N(y_i-g_i{)}^2+\lambda \sum _i^{N-2}(g_i-2g_{t+1}+g_{i+2}{)}^2\text{\hspace{1em}(1)}$$
即
$$\underset{\mathbf{g}}{\min }||\mathbf{y}-\mathbf{g}|{|}^2+{\lambda }^{\prime }||{\nabla }^2\mathbf{g}|{|}^2\text{\hspace{1em}(2)}$$
一方面是拟合：减小滤波后信号 $\mathbf{g}$ 与原信号的误差；另一方面是平滑，限制平滑信号的二阶差分大小。

求解

从(1)(2)式中可以看出，目标函数是未知量的二次函数，所以可以用最小二乘法求解。

首先，用矩阵来表示信号的差分：
$$\nabla \mathbf{g}={\left[\begin{array}{lllll}-1& 1& & & \\ & -1& 1& & \\ & & \cdots & & \\ & & & -1& 1\end{array}\right]}_{(N-1)\times N}\mathbf{g}$$
$${\nabla }^2\mathbf{g}={\left[\begin{array}{lllll}-1& 1& & & \\ & -1& 1& & \\ & & \cdots & & \\ & & & -1& 1\end{array}\right]}_{(N-2)\times (N-1)}\nabla \mathbf{g}$$
所以简单记为
$${\nabla }^2\mathbf{g}=D\mathbf{g}$$
问题(2)转换成矩阵形式
$$\begin{array}{ll}{\min }_{\mathbf{g}}& ||\mathbf{y}-\mathbf{g}|{|}^2+\lambda ||D\mathbf{g}|{|}^2\\ & \\ & =(\mathbf{y}-\mathbf{g}{)}^{\top }(\mathbf{y}-\mathbf{g})+\lambda (D\mathbf{g}{)}^{\top }(D\mathbf{g})\\ & \\ & ={\mathbf{g}}^{\top }(I+\lambda D^{\top }D)\mathbf{g}-2{\mathbf{y}}^{\top }\mathbf{g}+{\mathbf{y}}^{\top }\mathbf{y}\end{array}$$
对上式求梯度得：
$$(I+\lambda D^{\top }D)\mathbf{g}=\mathbf{y}$$
所以解得：
$$\mathbf{g}=(I+\lambda D^{\top }D{)}^{-1}\mathbf{y}$$

python 实现

def hp(y, lamb=10):
    def D_matrix(N):
        D = np.zeros((N-1,N))
        D[:,1:] = np.eye(N-1)
        D[:,:-1] -= np.eye(N-1)
        """D1
        [[-1.  1.  0. ...  0.  0.  0.]
         [ 0. -1.  1. ...  0.  0.  0.]
         [ 0.  0. -1. ...  0.  0.  0.]
         ...
         [ 0.  0.  0. ...  1.  0.  0.]
         [ 0.  0.  0. ... -1.  1.  0.]
         [ 0.  0.  0. ...  0. -1.  1.]]
        """
        return D
    N = len(ts)
    D1 = D_matrix(N)
    D2 = D_matrix(N-1)
    D = D2 @ D1
    g = np.linalg.inv((np.eye(N)+lamb*D.T@D))@ ts
    return g

测试

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

N = 100
t = np.linspace(1,10,N)
ts = np.sin(t) + np.cos(20*t) + np.random.randn(N)*0.1

plt.figure(figsize=(10,12))
for i,l in enumerate([0.1,1,10,100,1000, 10000]):
    plt.subplot(3,2,i+1)
    g = hp(ts,l)
    plt.plot(ts, label='original')
    plt.plot(g, label='filtered')
    plt.legend()
    plt.title('$\lambda$='+str(l))
plt.show()