import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as signal
def main():
x, y = np.loadtxt('spectrum.dat').T
logx, logy = np.log(x), np.log(y)
smooth_log = detrend_zero_phase(logx, logy)
smooth = np.exp(smooth_log)
fig, ax = plt.subplots()
ax.loglog(x, y, 'k-')
ax.loglog(x, smooth, 'r-')
plt.show()
def zero_phase(y):
# Low-pass filter...
b, a = signal.butter(3, 0.05)
# Filtfilt applies the filter twice to avoid phase shifts.
return signal.filtfilt(b, a, y)
def detrend_zero_phase(x, y):
# Fit a second order polynomial (Can't just use scipy.signal.detrend here,
# because we need to know what the trend is to add it back in.)
model = np.polyfit(x, y, 2)
trend = np.polyval(model, x)
# Apply a zero-phase filter to the detrended curve.
smooth = zero_phase(y - trend)
# Add the trend back in
return smooth + trend
main()
最简单的方法是在滤波之前对信号进行去趋势化处理。你所看到的边缘效应主要是由于信号不是静止的(即有一个斜率)。在
首先,让我们演示一下这个问题:
哎哟!边缘效应在数据的末尾(右边的大小)尤其糟糕,因为这是斜率最陡的地方。如果你在开始时有一个更陡的坡度,你会看到一个更强的边缘效应。在
好消息是有很多方法可以纠正这种情况。@ChristianK的答案显示了如何使用平滑样条来有效地实现低通滤波器。我将演示使用其他一些信号处理方法来完成相同的任务。哪一个“最好”取决于你的需要。平滑样条线是直接向前的。使用“更富想象力”的信号处理方法,可以精确控制过滤掉的频率。在
您的数据在log-log空间中看起来像一条抛物线,所以让我们用log-log空间中的二阶多项式来反趋势,然后应用过滤器。在
举个简单的例子:
^{pr2}$请注意,我们仍然有一些边缘效应。这是因为我使用的高斯滤波器会引起相移。如果我们真的想变得花哨,我们可以去趋势化,然后使用零相位滤波器进一步最小化边缘效应。在
如果您的所有数据在日志空间中的修改速度一样慢,我将执行以下操作:
例如:
这就产生了:
相关问题 更多 >
编程相关推荐