有效计算信号的 50Hz 内容
Efficiently compute 50Hz content of signal
问题陈述
我有一个较长的信号（454912 个样本），我想计算其中 50Hz 的估计值。速度在这里比精度更重要。预计 50Hz 的数量会随时间波动。该值需要代表整个记录，例如平均值。
上下文
信号是从 EEG 电极记录的。当脑电电极与头皮接触不良时，信号中会出现大量50Hz的电力线噪声。我想丢弃来自 EEG 电极的所有数据，这些电极的 50Hz 噪声比其他电极多得多。
尝试过的解决方案
解决问题并不难。从 FFT 到 Welch 的任何方法都可以用来估计功率谱：
import numpy as np
from scipy.signal import welch
# generate an example signal
sfreq = 128.
nsamples = 454912
time = np.arange(nsamples) / sfreq
x = np.sin(2 * np.pi * 50 * time) + np.random.randn(nsamples)
# apply Welch' method to the problem
fs, ps = welch(x, sfreq)
print 'Amount of 50Hz:', ps[np.searchsorted(fs, 50)]
然而，在这里计算所有频率的功率似乎是不必要的，我觉得有一个更有效的解决方案。类似于计算单个 DFFT 步骤的东西？与一些正弦小波卷积？
Welch's method 只是计算信号的多个重叠段的周期图，然后取各段的平均值。这有效地以分辨率换取频域中的降噪。
但是，为每个小段执行大量单独的 FFT 比为较大段计算较少的 FFT 成本更高。根据您的需要，您可能会使用 Welch 方法，但将信号分成更大的部分，and/or 它们之间的重叠较少（两者都会减少 PSD 的方差）。
from matplotlib import pyplot as plt
# default parameters
fs1, ps1 = welch(x, sfreq, nperseg=256, noverlap=128)
# 8x the segment size, keeping the proportional overlap the same
fs2, ps2 = welch(x, sfreq, nperseg=2048, noverlap=1024)
# no overlap between the segments
fs3, ps3 = welch(x, sfreq, nperseg=2048, noverlap=0)
fig, ax1 = plt.subplots(1, 1)
ax1.hold(True)
ax1.loglog(fs1, ps1, label='Welch, defaults')
ax1.loglog(fs2, ps2, label='length=2048, overlap=1024')
ax1.loglog(fs3, ps3, label='length=2048, overlap=0')
ax1.legend(loc=2, fancybox=True)
增加段大小并减少重叠量可以显着提高性能：
In [1]: %timeit welch(x, sfreq)
1 loops, best of 3: 262 ms per loop
In [2]: %timeit welch(x, sfreq, nperseg=2048, noverlap=1024)
10 loops, best of 3: 46.4 ms per loop
In [3]: %timeit welch(x, sfreq, nperseg=2048, noverlap=0)
10 loops, best of 3: 23.2 ms per loop
请注意，window 大小使用 2 的幂是个好主意，因为对长度为 2 的幂的信号进行 FFT 会更快。
更新
您可能会考虑尝试的另一个简单方法是使用以 50Hz 为中心的陷波滤波器对信号进行带通滤波。滤波后信号的包络可以衡量您的信号随时间包含多少 50Hz 功率。
from scipy.signal import filter_design, filtfilt
# a signal whose power at 50Hz varies over time
sfreq = 128.
nsamples = 454912
time = np.arange(nsamples) / sfreq
sinusoid = np.sin(2 * np.pi * 50 * time)
pow50hz = np.zeros(nsamples)
pow50hz[nsamples / 4: 3 * nsamples / 4] = 1
x = pow50hz * sinusoid + np.random.randn(nsamples)
# Chebyshev notch filter centered on 50Hz
nyquist = sfreq / 2.
b, a = filter_design.iirfilter(3, (49. / nyquist, 51. / nyquist), rs=10,
ftype='cheby2')
# filter the signal
xfilt = filtfilt(b, a, x)
fig, ax2 = plt.subplots(1, 1)
ax2.hold(True)
ax2.plot(time[::10], x[::10], label='Raw signal')
ax2.plot(time[::10], xfilt[::10], label='50Hz bandpass-filtered')
ax2.set_xlim(time[0], time[-1])
ax2.set_xlabel('Time')
ax2.legend(fancybox=True)
更新 2
看到@hotpaw2 的回答后，我决定尝试实施 Goertzel algorithm，只是为了好玩。不幸的是它是一个递归算法（因此不能随时间向量化），所以我决定自己写一个 Cython 版本：
# cython: boundscheck=False
# cython: wraparound=False
# cython: cdivision=True
from libc.math cimport cos, M_PI
cpdef double goertzel(double[:] x, double ft, double fs=1.):
"""
The Goertzel algorithm is an efficient method for evaluating single terms
in the Discrete Fourier Transform (DFT) of a signal. It is particularly
useful for measuring the power of individual tones.
Arguments
----------
x double array [nt,]; the signal to be decomposed
ft double scalar; the target frequency at which to evaluate the DFT
fs double scalar; the sample rate of x (same units as ft, default=1)
Returns
----------
p double scalar; the DFT coefficient corresponding to ft
See: <http://en.wikipedia.org/wiki/Goertzel_algorithm>
"""
cdef:
double s
double s_prev = 0
double s_prev2 = 0
double coeff = 2 * cos(2 * M_PI * (ft / fs))
Py_ssize_t N = x.shape[0]
Py_ssize_t ii
for ii in range(N):
s = x[ii] + (coeff * s_prev) - s_prev2
s_prev2 = s_prev
s_prev = s
return s_prev2 * s_prev2 + s_prev * s_prev - coeff * s_prev * s_prev2
这是它的作用：
freqs = np.linspace(49, 51, 1000)
pows = np.array([goertzel(x, ff, sfreq) for ff in freqs])
fig, ax = plt.subplots(1, 1)
ax.plot(freqs, pows, label='DFT coefficients')
ax.set_xlabel('Frequency (Hz)')
ax.legend(loc=1)
这太他妈快了：
In [1]: %timeit goertzel(x, 50, sfreq)
1000 loops, best of 3: 1.98 ms per loop
显然，只有当您只对单个频率而不是频率范围感兴趣时，这种方法才有意义。
对于单个正弦频率，您可以使用 Goertzel 算法或 Goertzel 滤波器，这是一种计算 DFT 或 FFT 结果的单个 bin 幅度的高效计算方法。
您可以 运行 在整个波形上使用此滤波器，或将其与 Welch 方法结合使用，方法是将一系列 Goertzel 滤波器的幅度输出相加，并选择滤波器长度，以便滤波器带宽不太窄（以涵盖 50 Hz 相对于您的采样率可能出现的轻微频率变化）。
Goertzel 滤波器通常与功率估计器结合使用，以确定所选频率的 SNR 是否有效。
问题陈述
我有一个较长的信号（454912 个样本），我想计算其中 50Hz 的估计值。速度在这里比精度更重要。预计 50Hz 的数量会随时间波动。该值需要代表整个记录，例如平均值。
上下文
信号是从 EEG 电极记录的。当脑电电极与头皮接触不良时，信号中会出现大量50Hz的电力线噪声。我想丢弃来自 EEG 电极的所有数据，这些电极的 50Hz 噪声比其他电极多得多。
尝试过的解决方案
解决问题并不难。从 FFT 到 Welch 的任何方法都可以用来估计功率谱：
import numpy as np
from scipy.signal import welch
# generate an example signal
sfreq = 128.
nsamples = 454912
time = np.arange(nsamples) / sfreq
x = np.sin(2 * np.pi * 50 * time) + np.random.randn(nsamples)
# apply Welch' method to the problem
fs, ps = welch(x, sfreq)
print 'Amount of 50Hz:', ps[np.searchsorted(fs, 50)]
然而，在这里计算所有频率的功率似乎是不必要的，我觉得有一个更有效的解决方案。类似于计算单个 DFFT 步骤的东西？与一些正弦小波卷积？
Welch's method 只是计算信号的多个重叠段的周期图，然后取各段的平均值。这有效地以分辨率换取频域中的降噪。
但是，为每个小段执行大量单独的 FFT 比为较大段计算较少的 FFT 成本更高。根据您的需要，您可能会使用 Welch 方法，但将信号分成更大的部分，and/or 它们之间的重叠较少（两者都会减少 PSD 的方差）。
from matplotlib import pyplot as plt
# default parameters
fs1, ps1 = welch(x, sfreq, nperseg=256, noverlap=128)
# 8x the segment size, keeping the proportional overlap the same
fs2, ps2 = welch(x, sfreq, nperseg=2048, noverlap=1024)
# no overlap between the segments
fs3, ps3 = welch(x, sfreq, nperseg=2048, noverlap=0)
fig, ax1 = plt.subplots(1, 1)
ax1.hold(True)
ax1.loglog(fs1, ps1, label='Welch, defaults')
ax1.loglog(fs2, ps2, label='length=2048, overlap=1024')
ax1.loglog(fs3, ps3, label='length=2048, overlap=0')
ax1.legend(loc=2, fancybox=True)
增加段大小并减少重叠量可以显着提高性能：
In [1]: %timeit welch(x, sfreq)
1 loops, best of 3: 262 ms per loop
In [2]: %timeit welch(x, sfreq, nperseg=2048, noverlap=1024)
10 loops, best of 3: 46.4 ms per loop
In [3]: %timeit welch(x, sfreq, nperseg=2048, noverlap=0)
10 loops, best of 3: 23.2 ms per loop
请注意，window 大小使用 2 的幂是个好主意，因为对长度为 2 的幂的信号进行 FFT 会更快。
更新
您可能会考虑尝试的另一个简单方法是使用以 50Hz 为中心的陷波滤波器对信号进行带通滤波。滤波后信号的包络可以衡量您的信号随时间包含多少 50Hz 功率。
from scipy.signal import filter_design, filtfilt
# a signal whose power at 50Hz varies over time
sfreq = 128.
nsamples = 454912
time = np.arange(nsamples) / sfreq
sinusoid = np.sin(2 * np.pi * 50 * time)
pow50hz = np.zeros(nsamples)
pow50hz[nsamples / 4: 3 * nsamples / 4] = 1
x = pow50hz * sinusoid + np.random.randn(nsamples)
# Chebyshev notch filter centered on 50Hz
nyquist = sfreq / 2.
b, a = filter_design.iirfilter(3, (49. / nyquist, 51. / nyquist), rs=10,
ftype='cheby2')
# filter the signal
xfilt = filtfilt(b, a, x)
fig, ax2 = plt.subplots(1, 1)
ax2.hold(True)
ax2.plot(time[::10], x[::10], label='Raw signal')
ax2.plot(time[::10], xfilt[::10], label='50Hz bandpass-filtered')
ax2.set_xlim(time[0], time[-1])
ax2.set_xlabel('Time')
ax2.legend(fancybox=True)
更新 2
看到@hotpaw2 的回答后，我决定尝试实施 Goertzel algorithm，只是为了好玩。不幸的是它是一个递归算法（因此不能随时间向量化），所以我决定自己写一个 Cython 版本：
# cython: boundscheck=False
# cython: wraparound=False
# cython: cdivision=True
from libc.math cimport cos, M_PI
cpdef double goertzel(double[:] x, double ft, double fs=1.):
"""
The Goertzel algorithm is an efficient method for evaluating single terms
in the Discrete Fourier Transform (DFT) of a signal. It is particularly
useful for measuring the power of individual tones.
Arguments
----------
x double array [nt,]; the signal to be decomposed
ft double scalar; the target frequency at which to evaluate the DFT
fs double scalar; the sample rate of x (same units as ft, default=1)
Returns
----------
p double scalar; the DFT coefficient corresponding to ft
See: <http://en.wikipedia.org/wiki/Goertzel_algorithm>
"""
cdef:
double s
double s_prev = 0
double s_prev2 = 0
double coeff = 2 * cos(2 * M_PI * (ft / fs))
Py_ssize_t N = x.shape[0]
Py_ssize_t ii
for ii in range(N):
s = x[ii] + (coeff * s_prev) - s_prev2
s_prev2 = s_prev
s_prev = s
return s_prev2 * s_prev2 + s_prev * s_prev - coeff * s_prev * s_prev2
这是它的作用：
freqs = np.linspace(49, 51, 1000)
pows = np.array([goertzel(x, ff, sfreq) for ff in freqs])
fig, ax = plt.subplots(1, 1)
ax.plot(freqs, pows, label='DFT coefficients')
ax.set_xlabel('Frequency (Hz)')
ax.legend(loc=1)
这太他妈快了：
In [1]: %timeit goertzel(x, 50, sfreq)
1000 loops, best of 3: 1.98 ms per loop
显然，只有当您只对单个频率而不是频率范围感兴趣时，这种方法才有意义。
对于单个正弦频率，您可以使用 Goertzel 算法或 Goertzel 滤波器，这是一种计算 DFT 或 FFT 结果的单个 bin 幅度的高效计算方法。
您可以 运行 在整个波形上使用此滤波器，或将其与 Welch 方法结合使用，方法是将一系列 Goertzel 滤波器的幅度输出相加，并选择滤波器长度，以便滤波器带宽不太窄（以涵盖 50 Hz 相对于您的采样率可能出现的轻微频率变化）。
Goertzel 滤波器通常与功率估计器结合使用，以确定所选频率的 SNR 是否有效。