如何滤波非周期函数

Question

我刚开始学习Python编程，想知道有没有办法为一个周期性函数创建一个高通滤波器，像这样：

import numpy as np
from scipy.signal import lfilter, firwin, butter
from pylab import figure, plot, show
sample_rate = .0167
nsamples = 480

F_1Hz = 1.38e-4
A_1Hz = 1.0

F_15Hz = .0011
A_15Hz = .5

t = np.arange(nsamples) / sample_rate
signal = A_1Hz * np.sin(2*np.pi*F_1Hz*t) + A_15Hz*np.sin(2*np.pi*F_15Hz*t)
signal[::120] = 2
figure(1)
plot(t,signal,'b')
show()

我想保留较高频率（0.0011赫兹）以及在特定位置的2的尖峰，但0.0011赫兹的幅度需要保持在0.5，而尖峰的幅度需要保持在2，所以不能进行归一化。此外，如果我让这个函数在非周期性间隔（比如只在素数位置出现尖峰）上有2的尖峰，我还能否正确地过滤它，并保持正确的幅度呢？

Answer 1

答案很可能是否定的。

这个直接的回答背后的原因是，你的尖峰（值为2）是在信号的顶部。如果你去掉任何东西，你的信号在尖峰处的强度可能会发生变化。

如果你能把这个：

signal[::120] = 2

变成

signal[::120] += 2

那么就可以构建这样的过滤器。你想去掉什么呢？是低于0.0011赫兹的东西吗？

Answer 2

一种可能的方法是使用自定义的高通滤波器。制作高通滤波器的简单方法是先从低通滤波器开始：

def lp_win_sinc(tw, fc, n):
    m = int(np.ceil( 2./tw) * 2) 
    samps = np.arange(m+1)

    shift = samps - m/2
    shift[m/2] = 1
    h = np.sin(2 * np.pi * fc * shift)/shift
    h[m/2] = 2 * np.pi * fc
    h = h * np.blackman(m+1)
    h = h / h.sum()
    s = np.zeros(n)
    s[:len(h)] = h
    return np.roll(s, -m/2)

然后构建一个简单的高通滤波器：

def hp_win_sinc(tw, fc, n):
    hp = -lp_win_sinc(tw, fc, n)
    hp[0] = hp[0] + 1
    return hp

（这些想法可以在 http://www.dspguide.com/pdfbook.htm 的窗口化 sinc 滤波器章节中找到。）

注意：这些是各自滤波器的脉冲响应。要将它们应用到你的数据上，你可以将脉冲与数据进行卷积，或者对你的数据和脉冲响应进行快速傅里叶变换（fft），然后对它们的乘积进行逆傅里叶变换（inverse fft）。在你的例子中，比如：

hp = hp_win_sinc(0.2, 0.001, len(signal))
f_hp = np.fft.rfft(hp)
f_d = np.fft.rfft(signal)
filt_sig = np.fft.irfft( f_hp * f_d)

绘制这个快速结果会得到：

过滤后的数据

根据你的具体应用，你可能只需要调整增益就能恢复到 2.0 和 0.5 的幅度。希望这对你有帮助。祝你好运！