Butter Filter 和 FFT 的派生结果不会随时间变化

Question

我是信号处理的新手，遇到过不知道对不对的情况。请纠正我，然后我会更新更多细节。

我的数据是here

我从手机（Samsung Galaxy Note 2，采样率 $\approx 99 Hz$）获取加速度计信号。我想分析从 $0.3 Hz$ 到 $5.0 Hz$

的频率

我的程序如下：

组合：假设一个传感器产生 3 个通道 $x$、$y$、$z$。组合是为了产生一个新的通道 $v = \sqrt{(x * x + y * y + z * z)}$
进行中值滤波：使信号平滑
Butter-worth 滤波器：我的截止频率是 $0.3 Hz$ 到 $5.0 Hz$
FFT

下图是我的演示，包含 120 个时间点和 4 个步骤：（可以在我的 video 中探索更多）

我观察到步骤 3 和 4 的结果没有改变，而信号随时间变化

我的问题是是否有什么可以确定这个结果是否正确？提前致谢

下面是我的代码用于应用过滤器

from __future__ import division
import numpy as np
from numpy.fft import rfft, rfftfreq
from numpy import absolute
import matplotlib.pyplot as plt
from scipy.signal import medfilt, hilbert
import pandas as pd

chunk = 120
LOW_CUT = 0.3
HIGH_CUT = 5.0
FS = 99
freqs = rfftfreq(chunk, 1 / FS)
_accel = pd.read_csv('data.csv')
for k, g in _accel.groupby(np.arange(len(_accel)) // chunk):
    _v = g['v'].values
    _v = medfilt(_v, 7)
    _v = butter_bandpass_filter(_v, LOW_CUT, HIGH_CUT, FS, order=4)
    v = 1 / chunk * absolute(rfft(_v))
    plt.stem(freqs, v)

更新1另一个link下载数据https://1drv.ms/u/s!At6qHz_a5mXhgp1KcAYpvsiJeTXsmg

更新 2 更新代码中的采样率 FS = 99

更新 3 将块大小增加到 512，plotted data again. Made a video of result without bandpass

Answer 1

我快速尝试一下这个问题，下面是我的代码片段

data = _accel['v'].tolist()    
Fs = 99
# remove the DC part, to help the plotting later
data = data - np.mean(data)
# Perform FFT for real data, on the whole 6000 samples, 
# using 4096 discrete frequencies, which is dense enough to capture 
# the frequency information within 0.3-5 Hz.    
fdata = rfft(data,4096)

# the frequencies we are looking at in the FFT
freqs = map(lambda x: float(x)*Fs/4096, range(0,4097))

# Plot
plt.plot(freqs[0:2049],fdata)
plt.xlabel('Frequency')
plt.show()

结果图确实包含您感兴趣的波段中的信息。 Plot of frequency magnitude

我猜你的问题是选择的块太小了。频域的分辨率为Fs/N，其中N为进行FFT的点数（通常为时域信号向量的长度）。所以，如果你想捕获 0.3-5Hz 范围内的信息，我假设你需要大约 0.2Hz 的分辨率，这意味着 N 应该至少为 500。你为 window 长度选择 120 显然是不够。

Answer 2

问题在于，每次在循环中处理新的数据块时，都会使用默认状态初始化过滤（如果之前的所有样本均为零，则对应于过滤器的状态）。结果，滤波器几乎没有时间在初始瞬态之后稳定下来（由从 "previous" 零到实际数据样本值的步骤引起），然后对下一个数据块再次执行相同的操作。

解决此问题的一种方法是在使用 FFT 处理数据块之前一次性过滤整个数据集：

_v = _accel['v'].values
_v = medfilt(_v, 7)
_v = butter_bandpass_filter1(_v, LOW_CUT, HIGH_CUT, FS, order=4)
for k in np.arange(1,len(_accel)//chunk):
     v = _v[chunk*k:chunk*(k+1)]
     v = 1 / chunk * absolute(rfft(v))
     plt.stem(freqs, v)

或者，您也可以跟踪过滤器状态（下面的zi）：

def butter_bandpass_filter(data, lowcut, highcut, fs, zi, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    if (zi == None):
        zi = lfilter_zi(b,a)
    y,zf = lfilter(b, a, data, zi=zi)
    return y,zf

zi = None
for k, g in _accel.groupby(np.arange(len(_accel)) // chunk):
    _v = g['v'].values
    _v = medfilt(_v, 7)
    _v,zi = butter_bandpass_filter(_v, LOW_CUT, HIGH_CUT, FS, zi, order=4)
    v = 1 / chunk * absolute(rfft(_v))
    plt.stem(freqs, v)

Butter Filter 和 FFT 的派生结果不会随时间变化

Derived Result from Butter Filter and FFT doesn't change over time

python

signals

fft

sensors