高效计算 50Hz 信号内容

Question

问题陈述

我有一个较长的信号（454912 个样本），我想计算其中的 50Hz 量的估计值。在这里，速度比精度更重要。预计存在的 50Hz 的数量会随着时间的推移而波动。该值需要代表整个记录，例如平均值。

上下文

信号是从 EEG 电极记录的。当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 步骤的东西？与一些正弦小波的卷积？

Answer 1

Welch's method 只计算信号的多个重叠段的周期图，然后取跨段的平均值。这有效地以分辨率换取频域中的降噪。

但是，为每个小段执行大量单独的 FFT 将比为较大的段计算更少的 FFT 更昂贵。根据您的需要，您可以不使用 Welch 的方法，但将信号分成更大的片段，和/或减少它们之间的重叠（这两种方法都会减少 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

请注意，最好使用 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

显然，这种方法只有在您只对单个频率而不是频率范围感兴趣时才有意义。

Answer 2

对于单个正弦频率，您可以使用 Goertzel 算法或 Goertzel 滤波器，这是一种计算 DFT 或 FFT 结果的单个 bin 幅度的高效计算方法。

您可以在整个波形上运行此滤波器，也可以将其与 Welch 方法结合使用，方法是将一系列 Goertzel 滤波器的幅度输出相加，并选择滤波器长度以使滤波器带宽不会太窄（以涵盖 50 Hz 与您的采样率的可能的轻微频率变化）。

Goertzel 滤波器通常与功率估算器结合使用，以确定所选频率的 SNR 是否有效。