<p><a href="http://en.wikipedia.org/wiki/Welch%27s_method" rel="nofollow noreferrer">Welch's method</a>只计算信号的多个重叠段的周期图,然后取各段的平均值。这有效地交换了频域中的降噪分辨率。在</p>
<p>然而,为每个小片段执行大量单独的fft要比为较大片段计算较少的fft要昂贵得多。根据您的需要,您可以不使用韦尔奇方法,但将您的信号分成更大的段,和/或它们之间的重叠较少(这两种方法都会减少PSD的方差)。在</p>
<pre><code>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)
</code></pre>
<p><img src="https://i.stack.imgur.com/pURiW.png" alt="enter image description here"/></p>
<p>增加段大小和减少重叠量可以显著提高性能:</p>
^{pr2}$
<p>注意,对于窗口大小使用2的幂是一个好主意,因为对长度为2的幂的信号进行FFT比较快。在</p>
<hr/>
<h2>更新</h2>
<p>另一个简单的事情,你可以考虑尝试只是带通滤波你的信号与一个陷波滤波器为中心的50赫兹。你的信号在50赫兹范围内的功率是多少。在</p>
<pre><code>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)
</code></pre>
<p><img src="https://i.stack.imgur.com/qkKPa.png" alt="enter image description here"/></p>
<hr/>
<h2>更新2</h2>
<p>看过@hotpaw2的答案后,我决定尝试实现<a href="http://en.wikipedia.org/wiki/Goertzel_algorithm" rel="nofollow noreferrer">Goertzel algorithm</a>,只是为了好玩。不幸的是,这是一个递归算法(因此不能随着时间的推移进行矢量化),所以我决定给自己写一个Cython版本:</p>
<pre><code># 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
</code></pre>
<p>它的作用是:</p>
<pre><code>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)
</code></pre>
<p><img src="https://i.stack.imgur.com/Y0pEX.png" alt="enter image description here"/></p>
<p>它非常快:</p>
<pre><code>In [1]: %timeit goertzel(x, 50, sfreq)
1000 loops, best of 3: 1.98 ms per loop
</code></pre>
<p>显然,只有当你只对一个频率感兴趣,而不是对一系列频率感兴趣时,这种方法才有意义。在</p>