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

Следующее выполняет среднее значение по выбранному региону:

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

Для правильного масштабирования значений мощности, которые были вычислены как abs( коэффициенты fft )**2 , вам нужно умножить на (2.0 / len(data))**2 ( теорема Парсеваля )

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

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

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

Если у вас действительно есть сигнал размером 15 ГБ, вы не сможете вычислить БПФ в приемлемое время. Вы можете избежать использования БПФ, если вам удобно приближать ваш частотный диапазон с помощью фильтра полосового диапазона. Обоснованием является формула суммирования Пуассона , которая утверждает, что сумма квадратов не изменяется БПФ (или: мощность сохраняется). Пребывание во временной области позволит увеличить время обработки пропорционально длине сигнала.

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

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

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

Шаг частоты dw = 2 * pi / (временной интервал) или 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) 

Я надеюсь, это поможет

С уважением