我正在尝试使用 scipy 编写一个简单的低通滤波器,但我需要帮助来定义参数。
我的时序数据中有350万条记录需要过滤,数据以1000Hz采样。
我正在使用 scipy 库中的 signal.firwin 和 signal.lfilter。
我在下面的代码中选择的参数根本不会过滤我的数据。取而代之的是,下面的代码只是简单地生成了一些在图形上看起来像完全相同的数据的东西,除了时间相位失真,它使图形向右移动了略小于 1000 个数据点(1 秒)。
在另一个软件程序中,通过图形用户界面命令运行低通冷杉滤波器产生的输出对于每 10 秒(10,000 个数据点)段具有相似的平均值,但其标准偏差显着降低,因此我们基本上失去了这个特定的噪声数据文件并将其替换为保留平均值同时显示不受更高频率噪声污染的长期趋势的东西。另一个软件的参数对话框包含一个复选框,允许您选择系数的数量,以便“根据样本大小和采样频率进行优化”。(我的是以 1000 赫兹收集的 350 万个样本,但我想要一个将这些输入用作变量的函数。)
*谁能告诉我如何调整下面的代码,以消除所有高于 0.05 赫兹的频率?* 我希望在图表中看到平滑的波浪,而不仅仅是我现在从下面的代码中得到的相同图表的时间失真。
class FilterTheZ0():
def __init__(self,ZSmoothedPylab):
#------------------------------------------------------
# Set the order and cutoff of the filter
#------------------------------------------------------
self.n = 1000
self.ZSmoothedPylab=ZSmoothedPylab
self.l = len(ZSmoothedPylab)
self.x = arange(0,self.l)
self.cutoffFreq = 0.05
#------------------------------------------------------
# Run the filter
#------------------------------------------------------
self.RunLowPassFIR_Filter(self.ZSmoothedPylab, self.n, self.l
, self.x, self.cutoffFreq)
def RunLowPassFIR_Filter(self,data, order, l, x, cutoffFreq):
#------------------------------------------------------
# Set a to be the denominator coefficient vector
#------------------------------------------------------
a = 1
#----------------------------------------------------
# Create the low pass FIR filter
#----------------------------------------------------
b = signal.firwin(self.n, cutoff = self.cutoffFreq, window = "hamming")
#---------------------------------------------------
# Run the same data set through each of the various
# filters that were created above.
#---------------------------------------------------
response = signal.lfilter(b,a,data)
responsePylab=p.array(response)
#--------------------------------------------------
# Plot the input and the various outputs that are
# produced by running each of the various filters
# on the same inputs.
#--------------------------------------------------
plot(x[10000:20000],data[10000:20000])
plot(x[10000:20000],responsePylab[10000:20000])
show()
return