滤波信号：如何限制滤波器最后一个输出点必须等于最后一个输入点

Question

请帮助我对信号处理的了解不足。

我想平滑一些数据。这是我的代码：

import numpy as np
from scipy.signal import butter, filtfilt

def testButterworth(nyf, x, y):
    b, a = butter(4, 1.5/nyf)
    fl = filtfilt(b, a, y)
    return fl

if __name__ == '__main__':
    positions_recorded = np.loadtxt('original_positions.txt', delimiter='\n')
    number_of_points = len(positions_recorded)
    end = 10
    dt = end/float(number_of_points)
    nyf = 0.5/dt

    x = np.linspace(0, end, number_of_points)
    y = positions_recorded

    fl = testButterworth(nyf, x, y)

我对结果很满意，除了一点： 返回值中的起点和终点等于输入的起点和终点，这对我来说绝对至关重要。如何引入此限制？

UPD 14 年 12 月 15 日 12:04：

我的原始数据是这样的 应用过滤器并放大图表的最后一部分会得到以下结果：

所以，目前我只关心必须等于原始点的最后一点。我尝试以这种方式将数据副本附加到原始列表的末尾：

结果如预期更糟。

然后我尝试以这种方式附加数据：

一个周期结束和下一个周期开始的切片看起来像这样：

Answer 1

问题是如何过滤数据，并要求过滤结果的左端点与数据的左端点匹配，右端点也一样。 （也就是说，一般来说，过滤后的结果应该接近大部分数据点，但不一定完全匹配其中任何一个，但是如果你需要在两个端点都匹配呢？）

为了使过滤后的结果与曲线的端点完全匹配，可以在曲线的任一端添加点填充并调整此填充的 y 位置，以便过滤器有效部分的端点完全匹配匹配原始数据的端点（没有填充）。

一般而言，这可以通过迭代解决方案、调整填充 y 位置直到末端对齐，或者通过计算一些值然后插值以确定所需的 y 位置来完成匹配的端点。我会做第二种方法。

这是我使用的代码，我将数据模拟为正弦波，两侧各有两个扁平部分（请注意，这些扁平部分 不是填充，但我只是在尝试使数据看起来有点像 OP）。

import numpy as np
from scipy.signal import butter, filtfilt
import matplotlib.pyplot as plt

####  op's code

def testButterworth(nyf, x, y):
    #b, a = butter(4, 1.5/nyf)
    b, a = butter(4, 1.5/nyf)
    fl = filtfilt(b, a, y)
    return fl

def do_fit(data):
    positions_recorded = data
    #positions_recorded = np.loadtxt('original_positions.txt', delimiter='\n')
    number_of_points = len(positions_recorded)
    end = 10
    dt = end/float(number_of_points)
    nyf = 0.5/dt

    x = np.linspace(0, end, number_of_points)
    y = positions_recorded

    fx = testButterworth(nyf, x, y)    
    return fx

###  simulate some data (op should have done this too!)

def sim_data():
    t = np.linspace(.1*np.pi, (2.-.1)*np.pi, 100)
    y = np.sin(t)
    c = np.ones(10, dtype=np.float)
    z = np.concatenate((c*y[0], y, c*y[-1]))
    return z

### code to find the required offset padding

def fit_with_pads(v, data, n=1):
    c = np.ones(n, dtype=np.float)
    z = np.concatenate((c*v[0], data, c*v[1]))
    fx = do_fit(z)
    return fx

def get_errors(data, fx):
    n = (len(fx)-len(data))//2
    return np.array((fx[n]-data[0], fx[-n]-data[-1]))

def vary_padding(data, span=.005, n=100):
    errors = np.zeros((4, n)) # Lpad, Rpad, Lerror, Rerror
    offsets = np.linspace(-span, span, n)
    for i in range(n):
        vL, vR = data[0]+offsets[i], data[-1]+offsets[i]
        fx = fit_with_pads((vL, vR), data, n=1)
        errs = get_errors(data, fx)
        errors[:,i] = np.array((vL, vR, errs[0], errs[1]))
    return errors

if __name__ == '__main__':
    data = sim_data()

    fx = do_fit(data)
    errors = vary_padding(data)

    plt.plot(errors[0], errors[2], 'x-')
    plt.plot(errors[1], errors[3], 'o-')

    oR = -0.30958
    oL =  0.30887
    fp = fit_with_pads((oL, oR), data, n=1)[1:-1]

    plt.figure()
    plt.plot(data, 'b')
    plt.plot(fx, 'g')
    plt.plot(fp, 'r')

    plt.show()

在这里，对于填充，我只在两边使用了一个点 (n=1)。然后我计算从第一个和最后一个数据点上下移动填充的一系列值的误差。

对于地块：
首先，我绘制偏移与误差（在拟合和所需数据值之间）。为了找到要使用的偏移量，我只是放大了两条线以找到 y 零交叉的 x 值，但要更准确地做到这一点，可以从这些数据中计算零交叉：

这是原始“数据”、拟合（绿色）和调整后拟合（红色）的图：

并放大 RHS：

这里的重点是红色（调整后的拟合）和蓝色（原始数据）端点匹配，即使纯拟合不匹配。

这是一种有效的方法吗？在各种选项中，这似乎是最合理的，因为通常不会对未显示的数据提出任何声明，而且显示区域也有一个准确应用的过滤器。例如，FFT 通常假设数据为零或超出边界是周期性的。当然，确切地说，应该解释所做的事情。

Answer 2

要做到这一点，您总是会以某种方式作弊，因为应用于真实数据的真实过滤器的行为方式并不符合您的要求。

欺骗您的数据的最佳方法之一是假设它是周期性的。这具有以下优点：1）它与您实际拥有的数据一致，并且您所做的所有更改都是将数据附加到您不知道的区域（因此假设它是周期性的，与其他任何东西一样合理 - 尽管可能违反一些未说明或隐含的假设）； 2）结果将与您的过滤器一致。

您通常可以通过将数据副本附加到真实数据的开头和结尾来解决这个问题，或者只是一小部分，具体取决于您的过滤器。

由于 FFT 假定数据无论如何都是周期性的，这通常是一种快速简便的方法，并且完全准确（而连接数据是对无限周期波形的估计）。以下是阶跃滤波器的 FFT 方法示例。

import numpy as np
import matplotlib.pyplot as plt

x = np.arange(0, 128)
y = (np.sin(.22*(x+10))>0).astype(np.float)

# filter
y2 = np.fft.fft(y)
f0 = np.fft.fftfreq(len(x))
y2[(f0<-.25) | (f0>.25)] = 0
y3 = abs(np.fft.ifft(y2))

plt.plot(x, y)
plt.plot(x, y3)
plt.xlim(-10, 140)
plt.ylim(-.1, 1.1)
plt.show()

请注意端点如何在两端相互弯曲，即使这与波形的周期性不一致（因为两端的线段都被截断了）。这也可以通过调整波形使末端相同来看到（这里我使用x+30而不是x+10，并且这里的末端不需要弯曲来匹配所以它们与末端保持水平数据。

另外请注意，要使端点实际上完全相等，您必须将此图扩展一个点（在任一端），因为它的周期与原始波形的波长完全相同。这样做不是临时的，结果将与您的分析完全一致，但只是代表一直假定无限重复的一个额外点。

最后，这个 FFT 技巧最适用于长度为 2ⁿ 的波形。在 FFT 中，其他长度可以是零填充的。在这种情况下，就像我一开始提到的那样，只对任一端进行连接可能是最好的方法。