在python中平滑曲线,边界没有错误？

首先,我们来说明问题：

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter1d

x,y = np.loadtxt('spectrum.dat').T

# Smooth with a guassian filter
smooth = gaussian_filter1d(y,10)

fig,ax = plt.subplots()
ax.loglog(x,y,color='black')
ax.loglog(x,smooth,color='red')
plt.show()

哎哟!边缘效应在数据的末端(右手大小)特别糟糕,因为这是斜率最陡的地方.如果你在开始时有一个更陡峭的坡度,那么你也会看到更强的边缘效应.

好消息是,有很多方法可以解决这个问题. @ChristianK.的答案显示了如何使用平滑样条来有效地执行低通滤波器.我将演示如何使用其他一些信号处理方法来完成同样的事情.哪个“最好”都取决于您的需求.平滑样条线是直截了当的.使用“发烧友”信号处理方法可以精确控制滤除的频率.

您的数据在对数日志空间中看起来像抛物线,所以让我们在对数日志空间中使用二阶多项式对其进行去趋势,然后应用过滤器.

作为一个简单的例子：

import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter1d

x,y = np.loadtxt('spectrum.dat').T

# Let's detrend by fitting a second-order polynomial in log space
# (Note that your data looks like a parabola in log-log space.)
logx,logy = np.log(x),np.log(y)
model = np.polyfit(logx,logy,2)
trend = np.polyval(model,logx)

# Smooth with a guassian filter
smooth = gaussian_filter1d(logy - trend,10)

# Add the trend back in and convert back to linear space
smooth = np.exp(smooth + trend)

fig,color='red')
plt.show()

请注意,我们仍然有一些边缘效应.这是因为我使用的高斯滤波器引起相移.如果我们真的想要花哨,我们可以解决问题,然后使用零相位滤波器来进一步减少边缘效应.

import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as signal

def main():
    x,y = np.loadtxt('spectrum.dat').T

    logx,np.log(y)
    smooth_log = detrend_zero_phase(logx,logy)
    smooth = np.exp(smooth_log)

    fig,ax = plt.subplots()
    ax.loglog(x,'k-')
    ax.loglog(x,'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,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()