时间序列中的降噪保持锐利的边缘

在来自功率计的时间序列中,过程和传感器都会产生噪声。为了识别步骤,我想在不牺牲边缘陡峭度的情况下过滤噪声。

想法是做一个 rolling(window).mean() => 杀死边缘或 rolling(window).median() => 但如果窗口大小需要小,这会产生谐波噪声问题。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# create a reference signal
xrng = 50
sgn = np.zeros(xrng)
sgn[10:xrng//2] = 1
sgn[xrng//2:xrng-10]=0.5

fig = plt.figure(figsize=(10,6))
plt.plot(sgn, label='raw') 

T=3    # period of the sine like noise (random phase shifts not modeled)
noise1 = (np.random.rand(xrng)-0.5)*0.2          # sensor noise
noise2 = np.sin(np.arange(xrng)*2*np.pi/T)*0.1   # harmonic noise 
sg_n = sgn + noise1 + noise2                     # noised signal 
plt.plot(sg_n, label='noised')

# 1. filter mean (good for hamonic)
mnfltr = np.ones(7)/7             
sg_mn = np.convolve(mnfltr,sg_n, 'same')
plt.plot(sg_mn, label='roll_mean')

# 2. filter median (good for edges)  
median = pd.Series(sg_n).rolling(9).median().shift(-4) 
plt.plot(median, label='roll_med') 

plt.legend()
plt.show()

输出如下: applying mean and median filters to noised signal

有没有办法结合两个过滤器来获得两种好处或任何其他方法?

stack overflow Noise reduction in time series keeping sharp edges
原文答案
author avatar

接受的答案

如果噪声的幅度不会掩盖步长,则可以使用完全不同的方法重建步进信号。

您的设置:

import numpy as np
import matplotlib.pyplot as plt

xrng = 50
sgn = np.zeros(xrng)
sgn[10:xrng//2] = 1
sgn[xrng//2:xrng-10]=0.5

fig = plt.figure(figsize=(10,6))
plt.plot(sgn, label='raw') 

T=3    # period of the sine like noise (random phase shifts not modeled)
noise1 = (np.random.rand(xrng)-0.5)*0.2          # sensor noise
noise2 = np.sin(np.arange(xrng)*2*np.pi/T)*0.1   # harmonic noise 
sg_n = sgn + noise1 + noise2                     # noised signal 
plt.plot(sg_n, label='noised')

噪声信号可以被数字化

bins = np.arange(-.25, 2, .5)
plt.plot((np.digitize(sg_n, bins)-1)/2, '.', markersize=8, label='reconstructed from noiced')
plt.legend();

结果:

reconstructed signal


答案:

作者头像

如果您不知道阶跃何时发生,也不知道它们的幅度,您可以使用 Rupture 包来确定信号稳定的窗口。然后对于每个稳定窗口,您计算均值并将其归因。

# building rupture to find the rupture points
algo = rpt.Pelt(model="rbf").fit(sg_n)
result = algo.predict(pen=2)

# Plot ruptures results
rpt.display(sg_n, result)

# smoothing signal piecewise
start_pt = 0
smoothed = []
for i in result:
    smoothed.append(np.tile(np.mean(sg_n[start_pt:i]), i-start_pt))
    start_pt = i
smoothed = np.hstack(smoothed)

# adding smoothed signal to ruptures plot
plt.plot(smoothed, '.', markersize=8, label='reconstructed from noised')
plt.legend()

这是结果:

Raw, noisy and smoothed data with the ruptures slicing results represented by the background color