python - 如何使用 np.fft.fft() 正确识别“峰值能量”并获得相关的绕组频率？

Question

我从概念上理解傅里叶变换。我写了一个简单的算法来计算变换，分解一个波并绘制它的各个组件。我知道它不是“快”，它也不能重建正确的幅度。它只是为了对机器背后的数学进行编码，它给了我很好的输出：

问题

1. 我该如何做类似的事情np.fft
2. 如何恢复 numpy 在引擎盖下选择的任何绕组频率？
3. 如何恢复使用变换找到的分量波的幅度？

我已经尝试了几件事。然而，当我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()

标签： pythonnumpysignalsfft

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 返回复数值，因此您需要平方和的 sqrt，就像您对 in 所做的mag那样transform。

• >>> 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为您提供反映在零赫兹附近的变换两侧。你可以在最后看到负频率。

>>> 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()

---

• 如果有一个 dc 组件并且您想查看它，您的过滤器将是freqs >= 0.
• 您的示例有 200 个数据点，因此您得到 100 (n/2) 个正频率，图表范围从零到一百赫兹，峰值在三和五。
• numpy.fft.rfft仅计算正频率。使用 numpy.fft.rfftfreq 获取频率。

---

对我来说，这是/现在是一个很棒的资源—— 《科学家和工程师数字信号处理指南》 ——它在我的办公桌上。