如何使用 np.fft.fft() 正确识别 "peak energy" 并获得相关的绕组频率？

Question

我从概念上理解傅里叶变换。我写了一个简单的算法来计算变换、分解波并绘制它的各个组件。我知道它不是 'fast'，而且它也没有重建正确的振幅。它只是为了编写机器背后的数学代码，它给了我这个不错的输出：

问题

如何使用 np.fft
如何恢复 numpy 在引擎盖下选择的任何绕组频率？
如何恢复使用变换找到的分量波的振幅？

我已经尝试了一些东西。然而，当我在与上面相同的波上使用 p = np.fft.fft(signal) 时，我得到了非常古怪的图，就像这个：

f1 = 3
f2 = 5
start = 0
stop = 1
sample_rate = 0.005
x = np.arange(start, stop, sample_rate)
y = np.cos(f1 * 2 * np.pi * x) + np.sin(f2 * 2 * np.pi *x)
p = np.fft.fft(y)
plt.plot(np.real(p))

或者如果我尝试使用 np.fft.freq() 来获得横轴的正确频率：

p = np.fft.fft(y)
f = np.fft.fftfreq(y.shape[-1], d=sampling_rate)
plt.plot(f, np.real(p))

作为最近的补充，我尝试实施@wwii 的建议得到了改进，但转换中的频率功率仍然关闭：

f1 = 3
f2 = 5
start = 0
stop = 4.5
sample_rate = 0.01
x = np.arange(start, stop, sample_rate)
y = np.cos(f1 * 2 * np.pi * x) + np.sin(f2 * 2 * np.pi *x)
p = np.fft.fft(y)
freqs= np.fft.fftfreq(y.shape[-1], d=sampling_rate)
q = np.abs(p)

q = q[freqs > 0]
f = freqs[freqs > 0]
peaks, _ = find_peaks(q)
peaks

plt.plot(f, q)
plt.plot(freqs[peaks], q[peaks], 'ro')
plt.show()

那么，我的问题是，如何使用 np.fft.fft 和 np.fft.fftfreqs 来获取与我的原始方法相同的信息？其次，如何从 fft（合成波的振幅）中恢复振幅信息。

我已经阅读了文档，但它一点帮助也没有。

对于上下文，这是我的天真方法：

def wind(timescale, data, w_freq):
    """
    wrap time-series data around complex plain at given winding frequency
    """
    return data * np.exp(2 * np.pi * w_freq * timescale * 1.j)


def transform(x, y, freqs):
    """ 
    Returns center of mass of each winding frequency
    """
    ft = []
    for f in freqs:
        mapped = wind(x, y, f)
        re, im = np.real(mapped).mean(), np.imag(mapped).mean()
        mag = np.sqrt(re ** 2 + im ** 2)
        ft.append(mag)
    
    return np.array(ft)

def get_waves(parts, time):
    """
    Generate sine waves based on frequency parts.
    """
    num_waves = len(parts)
    steps = len(time)
    waves = np.zeros((num_waves, steps))
    for i in range(num_waves):
        waves[i] = np.sin(parts[i] * 2 * np.pi * time)
    
    return waves
        
def decompose(time, data, freqs, threshold=None):
    """
    Decompose and return the individual components of a composite wave form.
    Plot each component wave. 
    """
    powers   = transform(time, data, freqs)
    peaks, _ = find_peaks(powers, threshold=threshold)
    
    plt.plot(freqs, powers, 'b.--', label='Center of Mass')
    plt.plot(freqs[peaks], powers[peaks], 'ro', label='Peaks')
    plt.xlabel('Frequency')
    plt.legend(), plt.grid()
    plt.show()
    
    return get_waves(freqs[peaks], time)

以及我用来生成绘图的信号设置：

# sample data plot: sin with frequencey of 3 hz. 
f1 = 3
f2 = 5
start = 0
stop = 1
sample_rate = 0.005
x = np.arange(start, stop, sample_rate)
y = np.cos(f1 * 2 * np.pi * x) + np.sin(f2 * 2 * np.pi *x)

plt.plot(x, y, '.')
plt.xlabel('time')
plt.ylabel('amplitude')
plt.show()

freqs = np.arange(0, 20, .5)
waves = decompose(x, y, freqs, threshold=0.12)

for w in waves:
    plt.plot(x, w)
plt.show()

Answer 1

f1 = 3
f2 = 5
start = 0
stop = 1
sample_rate = 0.005
x = np.arange(start, stop, sample_rate)
y = np.cos(f1 * 2 * np.pi * x) + np.sin(f2 * 2 * np.pi *x)
p = np.fft.fft(y)
freqs = np.fft.fftfreq(y.shape[0],sample_rate)

fft returns 复数值，因此您需要平方和的平方和，就像您在 transform 中对 mag 所做的那样。

>>> p[:2]
array([-1.42663659e-14+0.00000000e+00j, -1.77635684e-15+1.38777878e-17j])

```
q = np.absolute(p)
```

>>> q[:2]
array([1.77641105e-15, 2.70861628e-14])

fft 和 fftfreqs 为您提供围绕零赫兹反射的变换的两侧。你可以在最后看到 negative 频率。

>>> freqs[-10:]
array([-10.,  -9.,  -8.,  -7.,  -6.,  -5.,  -4.,  -3.,  -2.,  -1.])

您只关心正频率，因此您可以过滤它们并绘图。

q = q[freqs > 0]
freqs = freqs[freqs > 0]
plt.bar(freqs,q)
plt.show()
plt.close()

如果有直流分量并且您想查看它，您的过滤器将是 freqs >= 0。
你的例子有 200 个数据点，所以你得到 100 (n/2) 个正频率，图形范围从 0 到 100 Hz，峰值在 3 和 5。
numpy.fft.rfft 只计算正频率。使用 numpy.fft.rfftfreq 获取频率。

对我来说，这个 was/is 很棒的资源 - 科学家和工程师指南数字信号处理 - 它在我工作的办公桌上。

如何使用 np.fft.fft() 正确识别 "peak energy" 并获得相关的绕组频率？

How do I use np.fft.fft() to correctly identify "peak energy" and also get the associated winding frequencies?

python

signals

numpy

fft