在Python中计算与感兴趣频率F周围每个频率带的能量。

Question

我是一个信号处理的新手，在这个问题中，我想问如何获取与感兴趣的频率F周围每个频带的能量。我找到了一个公式，但我不知道如何在Python中实现它。这是公式和我的傅里叶变换图：

x = np.linspace(0,5,100)
y = np.sin(2*np.pi*x)

## fourier transform
f = np.fft.fft(y)
## sample frequencies
freq = np.fft.fftfreq(len(y), d=x[1]-x[0])
plt.plot(freq, abs(f)**2) ## will show a peak at a frequency of 1 as it should.

- markov zain

---

你离答案很近了，sum(abs(f[F-d:F+d]) ** 2)有什么问题吗？ - Mike

Answer 1

你离成功已经很接近了，就像Mike指出的那样，但这里有一种不同的方法更容易理解。你可以设置一个变量来保存过滤后的信号，并返回Af的一维数组，然后应用上述公式，这个公式非常简单（这些振幅的平方和）。

像这样过滤信号。

from scipy.signal import butter, lfilter
def butter_bandpass(lowcut, highcut, fs, order=5):
    nyq = 0.5 * fs
    low = lowcut / nyq
    high = highcut / nyq
    b, a = butter(order, [low, high], btype='band')
    return b, a    
def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter_bandpass(lowcut, highcut, fs, order=order)
    y = lfilter(b, a, data)
    return y

现在假设y是您的原始信号，您需要信号中5Hz的能量成分，

 #let fs = 250
 #let order = 5
 oneD_array_of_amps_of_fiveHz_component = butter_bandpass_filter(y, 4, 6, 250, 5)
 #calculate energy like this
 energy_of_fiveHz_comp = sum([x*2 for x in oneD_array_of_amps_of_fiveHz_component])

- Wajdan Ali

Answer 2

使用Numpy在一维Numpy数组中查找局部最大值/最小值的信号处理模块，您可以执行以下操作：

from scipy.signal import argrelextrema
import numpy as np

delta = 0.1
F = 1
X_delta = abs(freq - F) < delta
A_f_delta = abs(f[X_delta])**2
maximum_ix = argrelextrema(A_f, np.greater)
E_F_delta = np.sum(A_f[maximum_ix])/2
E_F_delta
2453.8235776144866

- maxymoo

---

我看不出这如何回答问题：公式中的能量E_F对频率F周围的平方傅里叶系数进行求和，涵盖一个频率宽度Delta。这个答案所做的完全不同：它将总和限制为精确的（单点）局部最大值（在这个解决方案中甚至没有Delta，因此它甚至无法工作！）。 - Eric O. Lebigot

噢，好的，我在我的答案中添加了F和delta，虽然我认为这在我的第一个答案中是隐含的......回答这个问题并不需要它，因为它只有一个频率。 - maxymoo

这个很好，但为什么所有以maximum_ix开头的行：你在第A_f_delta = …行回答了问题（实际上原始问题中称为E_F）。顺便提一下：根据您关于单个频率证明只查看峰值的论点，您甚至可以将所有平方幅度加起来，因为您的argrelextrema（）只会捕捉到重要振幅。在一般情况下，由于Parseval定理的存在，这显然是相当无用的（因为结果与在原始信号上执行相同操作的结果相同）。 - Eric O. Lebigot