获得范围内频率平均值的最快方法

Question

我是 python 以及信号处理方面的新手。我正在尝试计算信号的某些频率范围内的 mean 值。

我正在尝试做的事情如下：

import numpy as np
data = <my 1d signal>
lF = <lower frequency>
uF = <upper frequency>
ps = np.abs(np.fft.fft(data)) ** 2 #array of power spectrum

time_step = 1.0 / 2000.0

freqs = np.fft.fftfreq(data.size, time_step) # array of frequencies
idx = np.argsort(freqs) # sorting frequencies

sum = 0
c =0
for i in idx:
    if (freqs[i] >= lF) and (freqs[i] <= uF) :
        sum += ps[i]
        c +=1
avgValue = sum/c
print 'mean value is=',avgValue

我觉得计算还可以，但是像超过15GB的数据需要很多时间，处理时间呈指数级增长。有没有最快的方法可以使我能够以最快的方式在某个频率范围内获得功率谱的平均值。提前致谢。

编辑 1

我跟着this code计算了功率谱。

编辑 2

This 没有回答我的问题，因为它计算整个 array/list 的平均值，但我想要数组的一部分的平均值。

编辑 3

jez 使用掩码的解决方案减少了时间。实际上我有超过 10 个一维信号通道，我想以相同的方式处理它们，即分别在每个通道的范围内平均频率。我认为 python 循环很慢。有什么替代方案吗？像这样：

for i in xrange(0,15):
        data = signals[:, i]
        ps = np.abs(np.fft.fft(data)) ** 2
        freqs = np.fft.fftfreq(data.size, time_step)
        mask = np.logical_and(freqs >= lF, freqs <= uF )
        avgValue = ps[mask].mean()
        print 'mean value is=',avgValue

Answer 1

以下对所选区域进行平均：

mask = numpy.logical_and( freqs >= lF, freqs <= uF )
avgValue = ps[ mask ].mean()

为了正确缩放已计算为 abs(fft 系数)**2 的功率值，您需要乘以 (2.0 / len(data))**2 (Parseval's theorem)

请注意，如果您的频率范围包括奈奎斯特频率，它会变得有些复杂——为了获得精确的结果，该单个频率分量的处理将需要取决于 data.size 是偶数还是奇数）。所以为了简单起见，确保 uF 严格小于 max(freqs)。 [出于类似的原因，您应该确保 lF > 0。]

其原因解释起来很乏味，更正起来更乏味，但基本上：DC 分量在 DFT 中表示一次，而大多数其他频率分量是每次以半振幅表示两次（正频率和负频率）。更烦人的例外是奈奎斯特频率，如果信号长度为偶数，它会以全振幅表示一次，但如果信号长度为奇数，则以半振幅表示两次。如果您对振幅取平均值，所有这些都不会影响您：在线性系统中，被表示两次可以补偿半振幅。但是你正在平均 power，即在平均之前对值进行平方，所以这个补偿不起作用。

I've pasted my code 了解所有这些。此代码还展示了如何处理堆叠在一个 numpy 数组中的多个信号，这解决了关于在多通道情况下避免循环的后续问题。请记住向 numpy.fft.fft() 和我的 fft2ap().

提供正确的 axis 参数 both

Answer 2

如果您的数组是 2 的幂，则使用 FFT 会容易得多。当您执行 fft 时，频率范围从 -pi/timestep 到 pi/timestep（假设频率定义为 w = 2*pi/t，如果使用 f =1/t 表示，请相应地更改值）。您的频谱排列为 0 到 minfreqq--maxfreq 到零。您现在可以使用 fftshift 函数交换频率，您的频谱看起来像 minfreq -- DC -- maxfreq。现在您可以轻松确定所需的频率范围，因为它已经排序。

频率步长 dw=2*pi/(time span) or max-frequency/(N/2) 其中 N 是数组大小。 N/2点为直流或0频率。第 N 个位置是最大频率，现在您可以轻松确定您的范围

Lower_freq_indx=N/2+N/2*Lower_freq/max_freq
Higher_freq_index=N/2+N/2*Higher_freq/Max_freq
avg=sum(ps[lower_freq_indx:Higher_freq_index]/(Higher_freq_index-Lower_freq_index)

希望对您有所帮助

问候

Answer 3

如果您确实有 15 GB 大小的信号，您将无法在可接受的时间内计算 FFT。如果您可以接受通过带通滤波器来近似您的频率范围，则可以避免使用 FFT。理由是 Poisson summation formula，它表示平方和不会被 FFT 改变（或者：保留功率）。停留在时域会让处理时间与信号长度成比例上升。

以下代码设计了巴特沃斯带通滤波器，绘制了滤波器响应并过滤了样本信号：

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

dd = np.random.randn(10**4)  # generate sample data
T = 1./2e3  # sampling interval
n, f_s = len(dd), 1./T  # number of points and sampling frequency

# design band path filter:
f_l, f_u = 50, 500  # Band from 50 Hz to 500 Hz
wp = np.array([f_l, f_u])*2/f_s  # normalized pass band frequnecies
ws = np.array([0.8*f_l, 1.2*f_u])*2/f_s  # normalized stop band frequencies
b, a = signal.iirdesign(wp, ws, gpass=60, gstop=80, ftype="butter",
                        analog=False)

# plot filter response:
w, h = signal.freqz(b, a, whole=False)
ff_w = w*f_s/(2*np.pi)
fg, ax = plt.subplots()
ax.set_title('Butterworth filter amplitude response')
ax.plot(ff_w, np.abs(h))
ax.set_ylabel('relative Amplitude')
ax.grid(True)
ax.set_xlabel('Frequency in Hertz')
fg.canvas.draw()


# do the filtering:
zi = signal.lfilter_zi(b, a)*dd[0]
dd1, _ = signal.lfilter(b, a, dd, zi=zi)

# calculate the avarage:
avg = np.mean(dd1**2)
print("RMS values is %g" % avg)
plt.show()

阅读 Scipy's Filter design 的文档以了解如何修改过滤器的参数。

如果您想继续使用 FFT，请阅读 signal.welch and plt.psd 上的文档。 Welch 算法是一种有效计算信号功率谱密度的方法（有一些权衡）。

获得范围内频率平均值的最快方法

Fastest way to get average value of frequencies within range

python

performance

signal-processing

numpy