目标
我们希望去除信号中的趋势,使信号看起来是直的
原理
python实现
代码
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
def hp(y, lamb=10):
def D_matrix(N):
#(N-1,N) 元素全为0
D = np.zeros((N-1,N))
#后(N-1,N-1)对角线元素置1
D[:,1:] = np.eye(N-1)
#前(N-1,N-1)对角线元素置-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)
#@是矩阵乘法 D2.shape = (N-2,N-1) D1.shape = (N-1,N) 所以 D.shape = (N-2,N)
D = D2 @ D1
# G = [i+λF]*-1Y
#np.linalg.inv为矩阵求逆 套用公式 D.T@D = F F.shape=(N,N)
g = np.linalg.inv((np.eye(N)+lamb*D.T@D))@ ts
return g
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,10))
#尝试不同的λ
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()
结果
可以看到,λ越大,曲线越平滑