Расчет интенсивности пика БПФ

Question

для физической лаборатории наш профессор дал нам задание провести анализ частотного спектра щипковой струны. После получения звука нам дали скрипт для выполнения БПФ.

После БПФ у нас теперь есть несколько частотных пиков.

Затем он сказал нам, что мы должны рассчитать интенсивность каждого пика в отдельности. Я новичок в этом топи c, поэтому я прошу вашей помощи, как изменить данный код, чтобы получить в качестве пиковой интенсивности выхода, скажем, от 760 до 765 Гц.

Код здесь:

from scipy.fftpack import fft,ifft
import matplotlib.pyplot as plt
from scipy.signal import blackman

data =  np.loadtxt("mic.txt")
x = data[:,0]
y = data[:,1]

fy = fft(y)

print np.sum(y),"==",fy[0]

n = len(x)
t = x[-1]
fx = np.linspace(0,n/t,n)

plt.plot(fx[0:n/2],np.abs(fy[0:n/2]))

plt.xlabel("frequency (Hz)")
plt.show()

Буду признателен за помощь, Мэтью

Answer 1

Под интенсивностью пика я предполагаю, что вы имеете в виду величину каждого пика (конечно, это может быть нормализовано). Я использовал scipy.signal.find_peaks, чтобы получить пики в спектре. Поскольку FFT возвращает комплексные значения, я использовал его для абсолютного значения y. Это общий пример c, поэтому поиск максимумов в большинстве случаев не так уж тривиален. Существуют и другие инструменты обнаружения пиков, такие как scipy.signal.find_peaks_cwt или peakutils. Из scipy tutorials приведен минимальный рабочий пример:

import numpy as np
from scipy.fft import fft
import matplotlib.pyplot as plt

from scipy.signal import find_peaks

N = 600
T = 1.0 / 800.0
x = np.linspace(0.0, N*T, N)
y = np.sin(50.0 * 2.0 * np.pi * x) + 0.5 * np.sin(80.0 * 2.0 * np.pi * x)

yf = fft(y)
xf = np.linspace(0.0, 1.0 / (2.0 * T), N//2)

# finding the peaks in the absolute value of y
y_abs = 2.0 / N * np.abs(yf[0:N//2])
peakind, _ = find_peaks(y_abs)

plt.plot(xf, y_abs)

# plotting the peaks
plt.plot(xf[peakind], y_abs[peakind], 'k.')

plt.grid()
plt.show()

Величина пиков равна

>>> y_abs[peakind]
[0.70947072 0.4914645 ]