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

我是python和信号处理的新手.我试图计算mean信号的某些频率范围内的值.

我想做的是如下:

import numpy as np
data = 
lF = 
uF = 
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

我按照此代码计算功率谱.

编辑2

这不能回答我的问题,因为它计算整个数组/列表的平均值,但我希望平均数组的一部分.

编辑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

jez.. 7

以下内容对选定区域执行平均值:

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

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

请注意,如果您的频率范围包括奈奎斯特频率,它会略微繁琐 - 为了获得精确的结果,那么单个频率分量的处理将需要取决于data.size是偶数还是奇数.因此,为简单起见,请确保uF严格小于max(freqs).[出于类似的原因,你应该确保lF > 0.]

造成这种情况的原因很难解释,甚至更难以纠正,但基本上:直流分量在DFT中表示一次,而大多数其他频率分量在每次半幅度时表示两次(正频率和负频率) .更令人讨厌的例外是奈奎斯特频率,如果信号长度是偶数,则以全幅度表示一次,如果信号长度是奇数,则以半幅度表示两次.如果你是平均幅度,所有这些都不会影响你:在线性系统中,被表示两次补偿为半幅度.但是你的平均功率,即平均值之前的平方值,所以这种补偿不起作用.

我已经粘贴了我的代码来解决所有这些问题.此代码还说明了如何处理堆叠在一个numpy数组中的多个信号,这些信号解决了有关在多通道情况下避免循环的后续问题.请记住,正确的电源axis参数都来numpy.fft.fft()和我fft2ap().

以下内容对选定区域执行平均值:

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

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

请注意,如果您的频率范围包括奈奎斯特频率,它会略微繁琐 - 为了获得精确的结果,那么单个频率分量的处理将需要取决于data.size是偶数还是奇数.因此,为简单起见,请确保uF严格小于max(freqs).[出于类似的原因,你应该确保lF > 0.]

造成这种情况的原因很难解释,甚至更难以纠正,但基本上:直流分量在DFT中表示一次,而大多数其他频率分量在每次半幅度时表示两次(正频率和负频率) .更令人讨厌的例外是奈奎斯特频率,如果信号长度是偶数,则以全幅度表示一次,如果信号长度是奇数,则以半幅度表示两次.如果你是平均幅度,所有这些都不会影响你:在线性系统中,被表示两次补偿为半幅度.但是你的平均功率,即平均值之前的平方值,所以这种补偿不起作用.

我已经粘贴了我的代码来解决所有这些问题.此代码还说明了如何处理堆叠在一个numpy数组中的多个信号,这些信号解决了有关在多通道情况下避免循环的后续问题.请记住,正确的电源axis参数都来numpy.fft.fft()和我fft2ap().

如果您确实有15 GB大小的信号,则无法在可接受的时间内计算FFT.如果可以通过带通滤波器估算频率范围,则可以避免使用FFT.理由是泊松求和公式,其表明平方和不会被FFT改变(或者:保持功率).保持时域将使处理时间与信号长度成比例地增加.

以下代码设计了Butterworth带通滤波器,绘制滤波器响应并过滤采样信号:

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过滤器设计的文档,了解如何修改过滤器的参数.

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