黄油过滤器和 FFT 的派生结果不会随时间变化

Question

我是信号处理方面的新手，遇到过一种情况，我不确定它是否正确。请纠正我，然后我会更新更多细节。

我的数据是here

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

我的程序如下：

1. 组合：假设一个传感器产生 3 个通道 $x$、$y$、$z$。组合是产生一个新的通道 $v = \sqrt{(x * x + y * y + z * z)}$

2. 执行中值滤波：使信号平滑

3. 黄油价值过滤器：我的截止值是从 $0.3 Hz$ 到 $5.0 Hz$

4. 快速傅里叶变换

   下图是我的 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另一个下载数据的链接https://1drv.ms/u/s!At6qHz_a5mXhgp1KcAYpvsiJeTXsmg

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

更新 3 再次将块大小增加到 512，plotted 数据。制作了结果视频without bandpass

Answer 1

我快速尝试一下这个问题，下面是我的sn-p

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。您选择 120 作为窗口长度显然是不够的。

Answer 2

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

解决此问题的一种方法是在使用 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)