信息处理 - 频率变化且幅度恒定的正弦波的 FFT 不显示平台，为什么？ - 吾爱随笔录

频率变化且幅度恒定的正弦波的 FFT 不显示平台，为什么？

信息处理 fft 频谱 Python

2022-02-03 14:18:39

我想使用 python 来查看扬声器的音频响应是什么。这可以通过将用麦克风测量的扬声器输出与给定输入进行比较来完成。这意味着对于将在时域中的这两个信号，我需要执行 FFT。最终我想播放接近哈曼曲线的目标曲线。我试图生成信号以将其放入 wav 文件中，当检查信号是否符合我对 FFT 的要求时，我得到了一些奇怪的行为。因此，我尝试对具有恒定幅度并在 5 秒内以 44100 Hz 的采样率从 21-19700 Hz 变化的正弦波执行 FFT。我想使用 welch 方法，我将窗口作为样本的全长，因为在这种情况下，我生成的信号没有任何噪声。这给了我以下结果：

我使用了以下python代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import welch

samplerate = 44100
duration = 5
time = np.linspace(0,duration,samplerate*duration, endpoint= False)
frequency = np.linspace(21,19700,len(time),endpoint= False)
Const_wave = np.sin(np.pi*2*frequency*time)
f, check = welch(Const_wave,samplerate,nperseg=len(Wave))
plt.plot(f,check)
plt.xscale('log')
plt.xlabel('Frequency [Hz]')
plt.yscale('log')
plt.ylabel('Power')
plt.grid()
plt.show()

之后，我尝试了 scipy.fft 模块，它产生了不同但也不是预期的结果：

为此，我使用了以下代码：

f = rfftfreq(samplerate*duration,1/samplerate)
check = rfft(Const_wave)
plt.plot(f,check) 
plt.xscale('log')
plt.xlabel('Frequency [Hz]')
plt.yscale('log')
plt.ylabel('Power')
plt.grid()
plt.show()

我希望在所有频率上看到一个平坦的等功率平台，但这并没有发生。我也尝试对 REW 生成的线性扫描 wav 文件做同样的事情并得到类似的结果，那么我对平坦高原的期望是不正确的，还是这里有其他东西在起作用？

2个回答

您的 Python 代码有两个弱点，而不是拼写错误，但这些并不是您未能在 PSD 图中“显示平台”的原因。首先我们回顾一下这些“错别字”。

线

f, check = welch(Const_wave,samplerate,nperseg=len(Wave))

必读

f, check = welch(Const_wave,samplerate,nperseg=len(Const_wave))

否则，变量Const_wave必须是Wave。如果不更正您将代码复制到帖子时可能出现的“错字”，您的代码将无法编译。

然后，

time = np.linspace(0,duration,samplerate*duration, endpoint= False)
frequency = np.linspace(21,19700,len(time),endpoint= False)
Const_wave = np.sin(np.pi*2*frequency*time)

导致频率扫描高达 19700*2，而不是 19700，但也许你是故意这样做的。如果不是，请注意瞬时频率是相位的导数。达到一个频率 $f_{max}$ 在频率扫描中，公式为

p H 一种 s e (吨) = F_{米 一世 n} \cdot 吨 + (F_{米 一种 X} - F_{米 一世 n}) / d 你 r 一种 吨 一世 ○ n \cdot 吨^{2} / 2 s w e e p_w 一种 v e (吨) = 罪 (2 π \cdot p H 一种 s e (吨)) 一世 n s 吨 一种 n 吨 一种 n e ○ 你 s_F r e q (吨) = F_{米 一世 n} + ((F_{米 一种 X} - F_{米 一世 n}) / d 你 r 一种 吨 一世 ○ n) \cdot 吨

$phase(t) = f_{min}·t + (f_{max} - f_{min})/duration·t^2/2 \\ sweep\_wave(t) = \sin(2π·phase(t)) \\ instantaneous\_freq(t) = f_{min} + ((f_{max} - f_{min})/duration)·t$

你应该写代码

from scipy import signal 
...
fmin = 21
fmax = 19700
time = np.linspace(0,duration,samplerate*duration, endpoint= False)
Const_wave = np.sin(np.pi*2*(fmin*time + (fmax-fmin)/duration*time*time/2.0))

回到您的问题，为什么“具有不同频率和恒定幅度的正弦波不会显示一个平台”：答案可以在您选择 Welch 方法的参数中找到。您将 Welch 方法应用于使用 Hann 窗口处理的整个范围内的信号（请参阅 scipy.signal.welch 参考），并且您的 Hann 窗口“泄漏”到 PSD 中。您有两种选择来达到所需的“高原”。

1st：window用类似数组的值指定命名参数，像这样

f, check = signal.welch(Const_wave,samplerate,window=np.ones(len(Const_wave)))

2：保留默认（Hann）窗口（省略命名参数window），将命名参数的值减小nperseg到与扫描中的最大波周期（最小频率）相称的较低值，像这样

f, check = signal.welch(Const_wave,samplerate,nperseg=samplerate/frequency[0])

您也可以尝试调整noverlap参数值（默认为nperseg值的一半）。

在这两种情况下，图表都不会成为理想的“高原”，要理解差异——在第二种情况下相当重要——你应该更深入地研究韦尔奇的方法，特别是研究信号处理算法.

无论如何，您必须这样做，因为音频工程需要非常扎实的信号处理背景。将来，您参与更高级的 AE 项目可能（并且通常会）需要订阅Journal of the Audio Engineering Society。

同时，让我向您推荐与您在问题中提出的问题相关的参考资料（包括哈曼目标曲线）：

Swen Müller 和 Paulo Massarani 的扫描传递函数测量

Angelo Farina 使用扫描正弦技术同时测量脉冲响应和失真

Github 的项目 AutoEQ

这是产生频率斜坡的常见陷阱。对于给出的斜坡 $\cos(2\pi f(t) t)$ 斜坡函数为 $f(t)$ , $f(t)$ 不是所需的瞬时频率。请参阅这篇文章，其中详细说明了这一点：

频率斜坡的模拟

对于用于信道估计的优化频率啁啾生成，以及最小化与 FFT 相关的混叠效应，请参阅这篇文章：

如何使用快速傅里叶变换在波特图上绘制频率响应？

其它你可能感兴趣的问题

上一篇Scipy 和 Matlab 的巴特沃斯滤波器的传递函数似乎与理论不符下一篇为什么 ASK 调制会产生傅立叶边带？