如何通过FFT确定频率依赖幅度

Question

我正在尝试测量和计算一个音频设备的频率响应，方法是生成输入信号，然后用频率扫描来测量输出。

我想实现的目标可以用伪代码表示：

for each frequency f in (10-20k):
    generate reference signal s with frequency f
    async play s and record result r
    determine amplitude a of r using FFT
    add tuple (f,a) to result

然后在Python中：

import numpy as np
from matplotlib import pyplot as plt
import sounddevice as sd
import wave
from math import log10, ceil
from scipy.fft import fft, rfft, rfftfreq, fftfreq

SAMPLE_RATE = 44100     # Hertz
DURATION = 3            # Seconds

def generate_sine_wave(freq, sample_rate, duration):
    x = np.linspace(0, duration, int(sample_rate * duration), endpoint=False)
    frequencies = x * freq
    y = np.sin((2 * np.pi) * frequencies)
    return x, y

def main():
    # info about our device
    print(sd.query_devices(device="Studio 24c"))

    # default device settings
    sd.default.device = 'Studio 24c'
    sd.default.samplerate = SAMPLE_RATE

    f_start = 10
    f_end = 20000
    samples_per_decade = 10
    ndecades = ceil(log10(f_end) - log10(f_start))
    npoints = ndecades * samples_per_decade
    freqs = np.logspace(log10(f_start), log10(f_end), num=npoints, endpoint=True, base=10)
    
    measure_duration = 0.25 # seconds
    peaks = []
    for f in freqs:
        _, y = generate_sine_wave(f, SAMPLE_RATE, measure_duration)
        rec = sd.playrec(y, SAMPLE_RATE, input_mapping=[2], output_mapping=[2])
        sd.wait()
        yf = np.fft.rfft(rec, axis=0)
        yf_abs = 1 / rec.size * np.abs(yf)
        xf = np.fft.rfftfreq(rec.size, d=1./SAMPLE_RATE)
        peaks.append(np.max(yf_abs))

    plt.xscale("log")
    plt.scatter(freqs,peaks)
    plt.grid()
    plt.show()

if __name__ == "__main__":
    main()

在实际测量设备之前，我把输出信号循环回我的音频接口的输入，基本上是为了“校准”我正在做的事情。我本来期待在所有频率下，幅度应该是相等的。然而，我得到的结果却是这样的：

为什么即使生成的正弦波是一样的，所有的幅度还是差别这么大？在测量过程中，我没有对我的音频接口做任何更改。

Answer 1

我觉得这些差异看起来相对较小（大约1.0E-5分钟，最大1.6E-5）。你在频率范围的上限接近22.05 kHz的奈奎斯特频率，所以我想每个正弦波周期只有大约2个采样点，肯定会导致傅里叶变换不准确。

我知道你在看一个44.1 kHz的音频系统（比如CD），但也许可以试着把采样频率提高到500 kHz，看看这样能否得到更一致的结果。虽然这样会麻烦一些，但你也可以尝试在每个频率上采样n个完整周期。比如，对于10 Hz的频率，采样1秒钟以获得10个周期；对于1000 Hz的频率，采样10毫秒以获得10个周期。我觉得这样可以帮助你的基频对齐，这样你就只会得到傅里叶变换的虚部，这样可以进行更直接的比较。否则，在很多情况下，你会在正弦波的中间截断，导致引入真实的傅里叶变换成分。而且，由于你没有使用任何窗函数，你会得到超出预期频率的额外频率成分。我知道你在取绝对值，所以这应该不会让不完美的正弦波周期计数成为大问题，但在测试时可以试试，看看为什么你没有得到像预期那样平坦的频率响应。