6

我制作了一个 python 代码来使用 Weierstrass 变换来平滑给定的信号,这基本上是归一化高斯与信号的卷积。

代码如下:


#Importing relevant libraries  
from __future__ import division  
from scipy.signal import fftconvolve   
import numpy as np 

def smooth_func(sig, x, t= 0.002):   
    N = len(x)   
    x1 = x[-1]   
    x0 = x[0]    


# defining a new array y which is symmetric around zero, to make the gaussian symmetric.  
    y = np.linspace(-(x1-x0)/2, (x1-x0)/2, N)   
    #gaussian centered around zero.  
    gaus = np.exp(-y**(2)/t)     

#using fftconvolve to speed up the convolution; gaus.sum() is the normalization constant.  
    return fftconvolve(sig, gaus/gaus.sum(), mode='same')

如果我为一个阶跃函数运行此代码,它会平滑拐角,但在边界处它会解释另一个拐角并使其平滑,因此会在边界处产生不必要的行为。我用下面链接中显示的图来解释这一点。
边界效应

如果我们直接积分求卷积就不会出现这个问题。因此问题不在于 Weierstras 变换,因此问题在于 scipy 的 fftconvolve 函数。

要理解为什么会出现这个问题,我们首先需要了解 fftconvolve 在 scipy 中的工作原理。
fftconvolve 函数基本上使用卷积定理来加速计算。
简而言之:
卷积(int1,int2)=ifft(fft(int1)*fft(int2))
如果我们直接应用这个定理,我们不会得到想要的结果。为了得到想要的结果,我们需要对一个大小为 max(int1,int2) 的数组进行 fft。但这会导致不希望的边界效应。这是因为在 fft 代码中,如果 size(int) 大于大小(在其上采用 fft),它会将输入填充为零,然后采用 fft。这种零填充正是造成不希望的边界效应的原因。

你能建议一种消除这种边界效应的方法吗?

我试图通过一个简单的技巧将其删除。平滑函数后,我将平滑信号的值与边界附近的原始信号进行比较,如果它们不匹配,我将平滑函数的值替换为该点的输入信号。
如下:


i = 0 
eps=1e-3
while abs(smooth[i]-sig[i])> eps: #compairing the signals on the left boundary
    smooth[i] = sig[i]
    i = i + 1
j = -1

while abs(smooth[j]-sig[j])> eps: # compairing on the right boundary.
    smooth[j] = sig[j]
    j = j - 1

这个方法有个问题,因为使用了epsilon,smoothened函数有小的跳转,如下图:

上面的方法可以做些改变来解决这个边界问题吗?

4

2 回答 2

7

最好的方法可能是使用mode = 'valid'

The output consists only of those elements that do not rely on the zero-padding.

除非您可以包装信号,或者正在处理的信号是较大信号的摘录(在这种情况下:处理完整信号然后裁剪感兴趣的区域),否则在进行卷积时总是会产生边缘效应。你必须选择你想如何处理它们。使用mode = valid只是将它们裁剪掉,这是一个很好的解决方案。如果您知道信号始终是“阶梯状”的,则可以适当地扩展已处理信号的前端和末端。

于 2012-04-04T09:16:02.263 回答
3

对称滤波器内核在末端产生的内容取决于您假设数据超出末端的内容。

如果您不喜欢当前结果的外观(假设两端均为零),请尝试使用另一个假设来扩展数据,例如数据的反映或多项式回归延续。将两端的数据至少扩展过滤器内核长度的一半(除非您的扩展是零,这与非循环卷积所需的现有零填充一起免费提供)。然后在过滤后删除添加的结尾扩展,看看您是否喜欢假设的外观。如果不是,请尝试另一个假设。或者更好的是,如果有的话,使用超出末端的实际数据。

于 2012-04-06T00:32:27.243 回答