NumPy快速傅里叶变换（FFT）无法处理在Audacity中生成的正弦波。

Question

我正在尝试使用Python的NumPy库进行频率分析。我有两个.wav文件，它们都包含一个440 Hz正弦波。一个是我用NumPy正弦函数生成的，另一个是在Audacity中生成的。FFT适用于由Python生成的那个文件，但对Audacity生成的文件无效。

以下是两个文件的链接：

无法工作的文件：440_audacity.wav 工作的文件：440_gen.wav 这是我用来进行傅里叶变换的代码：

import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as wave

infile = "440_gen.wav"
rate, data = wave.read(infile)

data = np.array(data)

data_fft = np.fft.fft(data)
frequencies = np.abs(data_fft)

plt.subplot(2,1,1)
plt.plot(data[:800])
plt.title("Original wave: " + infile)

plt.subplot(2,1,2)
plt.plot(frequencies)
plt.title("Fourier transform results")

plt.xlim(0, 1000)

plt.tight_layout()

plt.show()

我有两个16位PCM.wav文件，一个来自Audacity，另一个是使用NumPy正弦函数创建的。NumPy生成的文件给出了以下（正确的）结果，440Hz处有尖峰: 虽然Audacity生成的波形看起来相同，但傅里叶变换没有任何结果: 我承认我在这里很茫然。这两个文件应该实际上包含相同的数据。它们以相同的方式编码，并且上图中的波形看起来相同。

下面是用于生成有效文件的代码:

import numpy as np
import wave
import struct
import matplotlib.pyplot as plt
from operator import add

freq_one = 440.0
num_samples = 44100
sample_rate = 44100.0
amplitude = 12800

file = "440_gen.wav"

s1 = [np.sin(2 * np.pi * freq_one * x/sample_rate) * amplitude for x in range(num_samples)]

sine_one = np.array(s1)

nframes = num_samples
comptype = "NONE"
compname="not compressed"
nchannels = 1
sampwidth = 2

wav_file = wave.open(file, 'w')
wav_file.setparams((nchannels, sampwidth, int(sample_rate), nframes, comptype, compname))

for s in sine_one:
    wav_file.writeframes(struct.pack('h', int(s)))

- darkfire613

---

你能分享一下440_gen.wav吗？ - jwalton

在终端中尝试使用aplay 440_gen.wav录制并查看WAV文件是否正确。 - Angelo Mendes

@Ralph，我正在努力上传两个文件以进行附加。 - darkfire613

@darkfire613 谢谢。只需要由audacity生成的文件，因为你已经包含了生成python生成的.wav的代码。此外，“波形在上面的图表上看起来完全相同”并不是正确的。虽然它们的频率相同，但幅度不同（尽管这不应该影响，并且我认为这不是问题所在）。 - jwalton

@Ralph，我已经在两个文件中添加了Dropbox链接。 - darkfire613

Answer 1

让我解释一下为什么你的代码不起作用。以及为什么它在使用 [:44100] 时会起作用。

首先，你有不同的文件：

440_gen.wav      = 1 sec and 44100  samples (counts)        
440_audacity.wav = 5 sec and 220500 samples (counts)

由于在FFT中，对于440_gen.wav，您使用了参考点N = 44100和采样率44100，因此您的频率分辨率为1 Hz（bin按照1 Hz增量跟随）。 因此，在图表上，每个FFT样本对应于等于1 Hz的delta。 plt.xlim(0,1000)仅对应于0-1000 Hz范围。

但是，对于FFT中的440_audacity.wav，您使用参考点N = 220500和采样率44100。您的频率分辨率为0.2 Hz（bin按照0.2 Hz增量跟随）-在图表上，每个FFT样本对应于以0.2 Hz增量的频率（最小值 - 最大值=+(-) 22500 Hz）。 plt.xlim(0,1000)仅对应于1000x0.2 = 0-200 Hz范围。 这就是为什么结果不可见-它不在此范围内。 plt.xlim(0,5000)将纠正您的情况并将范围扩展到0-1000 Hz。 jwalton提供的解决方案[: 44100]实际上只是强制FFT使用N = 44100。这重复了计算440_gen.wav的情况。

更正确的解决方案是在代码中使用N（窗口大小）参数和np.fft.fftfreq()函数。

以下是示例代码。

我还推荐一篇优秀的文章：https://realpython.com/python-scipy-fft/

import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as wave

N = 44100 # added

infile = "440_audacity.wav"
rate, data = wave.read(infile)

data = np.array(data)

data_fft = np.fft.fft(data, N)  # added N
frequencies = np.abs(data_fft)
x_freq = np.fft.fftfreq(N, 1/44100)  # added

plt.subplot(2,1,1)
plt.plot(data[:800])
plt.title("Original wave: " + infile)

plt.subplot(2,1,2)
plt.plot(x_freq, frequencies)  # added x_freq 
plt.title("Fourier transform results")

plt.xlim(0, 1000)
plt.tight_layout()
plt.show()

- Konyukh Fyodorov

---

优秀的回答。欢迎来到 Stack Overflow！ - Cris Luengo

顺便说一句，这也可以在不仅限于仅截取第一秒信号的情况下工作，将“N”设置为信号长度即可。 - Cris Luengo

谢谢你说了那句话，Cris。 - Konyukh Fyodorov

Answer 2

自回答这个问题以来，@Konyukh Fyodorov提供了更好的、合理的解决方案（如下）。

---

以下方法适用于我，并按预期生成了图表。不幸的是，我无法完全理解为什么这样可以，但我分享这个解决方案希望它能帮助其他人跨越难关。

import numpy as np
import matplotlib.pyplot as plt
import scipy.io.wavfile as wave

infile = "440_gen.wav"
rate, data = wave.read(infile)

data = np.array(data)

# Use first 44100 datapoints in transform
data_fft = np.fft.fft(data[:44100])
frequencies = np.abs(data_fft)

plt.subplot(2,1,1)
plt.plot(data[:800])
plt.title("Original wave: " + infile)

plt.subplot(2,1,2)
plt.plot(frequencies)
plt.title("Fourier transform results")

plt.xlim(0, 1000)

plt.tight_layout()

plt.show()

- jwalton

---

谢谢！这对我也有效，使用了我之前在Audacity中创建并且遇到困难的多个不同示例文件。限制样本数量是否有什么原因呢？ - darkfire613

太好了，我很高兴这对你起作用了！但是我真的不确定为什么它有效；我已经很久没有使用傅里叶变换了，我也记不清它们的工作原理了。我需要再想一想。 - jwalton

“这对我有用”这样的回答并没有指明代码中发生了什么变化，也没有解释为什么会起作用，这并不能帮助任何人学习。这是一种修补代码但不增加 Stack Overflow 知识库的回答。 - Cris Luengo

@Cris Luengo，恕我直言，不是所有问题都可以线性解决。有时在理解“为什么”之前看到修复方法可以帮助弥合知识差距并促进学习。我无法弥合这个差距，但我分享了我的补丁，希望它能帮助其他人跨越这个鸿沟。您可以选择将解决方案的“纯度优先于实用性”，但并非每个人都会这样做。对我来说，功利主义者会因为个人哲学的限制而等待2年，但我会很高兴使用有效的解决方案。 - jwalton