目的

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

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

优化函数

$\min_g \sum_i^N (y_i - g_i)^2 + \lambda \sum_i^{N-2} (g_i - 2g_{t+1} + g_{i+2})^2 \tag{1}$
即
$\min_\mathbf{g} ||\mathbf{y} - \mathbf{g}||^2 + \lambda' ||\nabla^2\mathbf{g}||^2 \tag{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}$
对上式求梯度得：
$\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 DN = len(ts)D1 = D_matrix(N)D2 = D_matrix(N-1)D = D2 @ D1g = np.linalg.inv((np.eye(N)+lamb*D.T@D))@ tsreturn g

测试

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inlineN = 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()