Самый быстрый способ получить среднее значение частот в пределах диапазона

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 для расчета спектра мощности.

EDIT 2

This не отвечает на мой вопрос, как он вычисляет среднее значение по всему массиву/список, но я хочу, имею в виду по части массива.

РЕДАКТИРОВАТЬ 3

Решение от Jez использования маски сокращает время. На самом деле у меня более 10 каналов 1D-сигнала, и я хочу обрабатывать их таким же образом, то есть средние частоты в диапазоне каждого канала отдельно. Я думаю, что петли питона медленны. Есть ли альтернатива для этого? Как это:

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)

Обратите внимание, что это немного затруднительно, если ваш частотный диапазон включает частоту Найквиста - для получения точных результатов, handli ng этой одночастотной составляющей тогда будет зависеть от того, является ли data.size четным или нечетным). Поэтому для простоты убедитесь, что uF строго меньше max(freqs). [По тем же причинам, вы должны обеспечить lF > 0.]

Причины для этого нудно объяснять и еще более утомительным для коррекции, но в основном: постоянная составляющая представлена ​​раз в ДПФ, в то время как большинство других частотных составляющих представлены дважды (положительная частота и отрицательная частота) с половинной амплитудой каждый раз. Еще более раздражающим исключением является частота Найквиста, которая представлена ​​один раз на полной амплитуде, если длина сигнала четная, но вдвое меньше половины амплитуды, если длина сигнала нечетна. Все это не повлияло бы на вас, если бы вы усредняли амплитуда: в линейной системе, представленной дважды, компенсируется наличие половины амплитуды. Но вы усредняете мощность, т. Е. Возведение в квадрат значений перед усреднением, поэтому эта компенсация не работает.

I've pasted my code для grokking все это.Этот код также показывает, как вы можете работать с несколькими сигналами, уложенными в один массив numpy, в котором рассматривается ваш следующий вопрос об исключении циклов в многоканальном случае. Не забудьте указать правильный аргумент axisкак, так и numpy.fft.fft() и моим fft2ap().

Answer 2

Гораздо проще работать с FFT, если ваши массивы имеют мощность 2. Когда вы выполняете fft, частоты варьируются от -pi/timestep до pi/timestep (если предположить, что частота определяется как w = 2 * pi/t , измените значения соответственно, если вы используете представление f = 1/t). Ваш спектр устроен как от 0 до minfreqq - maxfreq до нуля. теперь вы можете использовать функцию fftshift для замены частот, и ваш спектр выглядит как minfreq - DC - maxfreq. теперь вы можете легко определить желаемый диапазон частот, потому что он уже отсортирован.

Шаг частоты dw = 2 * pi/(временной интервал) или максимальная частота/(N/2), где N - размер массива. Точка N/2 - постоянная или 0 частота. Nth положения максимальная частота теперь вы можете легко определить ваш диапазон

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 Гб, вы не сможете для расчета БПФ в приемлемое время. Вы можете избежать использования БПФ, если вам удобно приближать ваш частотный диапазон с помощью фильтра полосового диапазона. Обоснованием является Poisson summation formula, в котором говорится, что сумма квадратов не изменяется БПФ (или: мощность сохраняется). Пребывание во временной области позволит увеличить время обработки пропорционально длине сигнала.

Следующий код конструкции фильтра Баттерворта пути полосы, участки отклика фильтра и фильтрует сигнал выборки:

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. Алгоритм Уэлша представляет собой метод для эффективного вычисления спектральной плотности мощности сигнала (с некоторыми компромиссами).