有效计算符号的50Hz含量

Question

问题陈述

我有一个很长的信号（454912个样本），我想计算一个50赫兹信号的估计值。在这里，速度比精度更重要。50赫兹的电流预计会随时间而波动。该值需要代表整个记录，例如平均值。在

上下文

信号是从脑电图电极记录下来的。当脑电图电极与头皮接触不良时，信号中会产生大量的50Hz电力线噪声。我想丢弃所有来自EEG电极的数据，这些电极比其他电极有过多的50赫兹噪声。在

尝试的解决方案

解决问题并不难。从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步骤一样的东西？用一些sinoid小波进行卷积？在

Answer 1

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

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

增加段大小和减少重叠量可以显著提高性能：

^{pr2}$

注意，对于窗口大小使用2的幂是一个好主意，因为对长度为2的幂的信号进行FFT比较快。在

---

更新

另一个简单的事情，你可以考虑尝试只是带通滤波你的信号与一个陷波滤波器为中心的50赫兹。你的信号在50赫兹范围内的功率是多少。在

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结果的单个单元幅值的有效方法。在

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

通常，Goertzel滤波器与功率估计器一起使用，以确定所选频率的SNR是否有效。在