获得范围内频率平均值的最快方法[关闭]

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 没有回答我的问题，因为它计算整个数组/列表的平均值，但我想要部分数组的平均值。

编辑 3

jez 使用面具的解决方案减少了时间。实际上，我有超过 10 个 1D 信号通道，我想以相同的方式处理它们，即分别在每个通道范围内的平均频率。我认为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。]

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

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

Answer 2

如果您确实有 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 和 plt.psd 上的文档。 Welch 算法是一种有效计算信号功率谱密度的方法（有一些折衷）。

Answer 3

如果您的数组是 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) 或 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)

我希望这会有所帮助

问候