我是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
我遵循此代码来计算功率谱。
编辑2
这没有回答我的问题,因为它计算整个数组/列表的均值,但我想对部分数组求平均值。
编辑3
通过使用面膜的解决方案可以减少时间。实际上,我有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
以下对选定区域进行平均:
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
。
这样做的原因很繁琐,但更难以纠正,但基本上:DC分量在DFT中仅表示一次,而大多数其他频率分量则在两次半振幅时分别表示两次(正频率和负频率) 。更加令人讨厌的例外是奈奎斯特频率,如果信号长度为偶数,则以全振幅表示一次,但如果信号长度为奇数,则以半振幅表示两次。如果要对振幅求平均,那么所有这些都不会影响您:在线性系统中,两次代表会补偿一半的振幅。但是您要平均功率,即在平均之前对值进行平方,因此这种补偿无法实现。
我已经粘贴了所有代码。这段代码还显示了如何使用一个numpy数组堆叠多个信号,从而解决了有关避免在多通道情况下出现循环的后续问题。请记住,正确的电源axis
参数都来numpy.fft.fft()
和我fft2ap()
。
本文收集自互联网,转载请注明来源。
如有侵权,请联系[email protected] 删除。
我来说两句