原始 PCM 的 FFT 图在 python 中出现更高频率的错误
FFT plot of raw PCM comes wrong for higher frequency in python
我在这里使用 numpy 的 fft 函数来绘制从 10000Hz 正弦波生成的 PCM 波的 fft。但是我得到的情节幅度是错误的。
使用我在控制台本身打印的 fftfreq 函数，频率正在变得正确。我的 python 代码在这里。
import numpy as np
import matplotlib.pyplot as plt
frate = 44100
filename = 'Sine_10000Hz.bin' #signed16 bit PCM of a 10000Hz sine wave
f = open('Sine_10000Hz.bin','rb')
y = np.fromfile(f,dtype='int16') #Extract the signed 16 bit PCM value of 10000Hz Sine wave
f.close()
####### Spectral Analysis #########
fft_value = np.fft.fft(y)
freqs = np.fft.fftfreq(len(fft_value)) # frequencies associated with the coefficients:
print("freqs.min(), freqs.max()")
idx = np.argmax(np.abs(fft_value)) # Find the peak in the coefficients
freq = freqs[idx]
freq_in_hertz = abs(freq * frate)
print("\n\n\n\n\n\nfreq_in_hertz")
print(freq_in_hertz)
for i in range(2):
print("Value at index {}:\t{}".format(i, fft_value[i + 1]), "\nValue at index {}:\t{}".format(fft_value.size -1 - i, fft_value[-1 - i]))
#####
n_sa = 8 * int(freq_in_hertz)
t_fft = np.linspace(0, 1, n_sa)
T = t_fft[1] - t_fft[0] # sampling interval
N = n_sa #Here it is n_sample
print("\nN value=")
print(N)
# 1/T = frequency
f = np.linspace(0, 1 / T, N)
plt.ylabel("Amplitude")
plt.xlabel("Frequency [Hz]")
plt.xlim(0,15000)
# 2 / N is a normalization factor Here second half of the sequence gives us no new information that the half of the FFT sequence is the output we need.
plt.bar(f[:N // 2], np.abs(fft_value)[:N // 2] * 2 / N, width=15,color="red")
输出来自控制台（我在这里粘贴的只有最少的印刷品）
freqs.min(), freqs.max()
-0.5 0.49997732426303854
freq_in_hertz
10000.0
Value at index 0: (19.949569768991054-17.456031216294324j)
Value at index 44099: (19.949569768991157+17.45603121629439j)
Value at index 1: (9.216783424692835-13.477631008179145j)
Value at index 44098: (9.216783424692792+13.477631008179262j)
N value=
80000
频率提取正确，但在情节中我做的事情不正确，我不知道。
正在更新作品：
- 当我将
n_sa = 10 * int(freq_in_hertz) 行中的乘数 10 更改为 5 时，我得到了正确的情节。对不对我看不懂
- 在
plt.xlim(0,15000) 行中，如果我再次将最大值增加到 20000 则不会绘图。直到 15000 它都在正确绘制。
- 我使用 Audacity 工具生成了此
Sine_10000Hz.bin，在该工具中我生成了频率为 10000Hz、持续时间为 1 秒、采样率为 44100 的正弦波。然后我将此音频导出为带有无标头的签名 16 位（意味着原始 PCM） .我可以使用这个脚本重新生成这个正弦波。我也想计算这个的 FFT。所以我期望在 10000Hz 处有一个峰值，振幅为 32767。你可以看到我在行 n_sa = 8 * int(freq_in_hertz) 中更改了倍增因子 8 而不是 10。因此它起作用了。但是振幅显示不正确。我会在这里附上我的新图
我不确定你到底想做什么，但我怀疑 Sine_10000Hz.bin 文件不是你想的那样。
是否可能包含多个频道（左&右）？
它真的是带符号的 16 位整数吗？
在 numpy 中创建 16 位整数的 10kHz 正弦波并不难。
import numpy as np
import matplotlib.pyplot as plt
n_samples = 2000
f_signal = 10000 # (Hz) Signal Frequency
f_sample = 44100 # (Hz) Sample Rate
amplitude = 2**3 # Arbitrary. Must be > 1. Should be > 2. Larger makes FFT results better
time = np.arange(n_samples) / f_sample # sample times
# The signal
y = (np.sin(time * f_signal * 2 * np.pi) * amplitude).astype('int16')
如果绘制 30 个信号点，您可以看到每个周期大约有 5 个点。
plt.plot(time[:30], y[:30], marker='o')
plt.xlabel('Time (s)')
plt.yticks([]); # Amplitude value is artificial. hide it
如果绘制来自 Sine_10000Hz.bin 的 30 个数据样本，它每个周期大约有 5 个点吗？
这是我尝试重新创建我理解的 FFT 工作。
fft_value = np.fft.fft(y) # compute the FFT
freqs = np.fft.fftfreq(len(fft_value)) * f_sample # frequencies for each FFT bin
N = len(y)
plt.plot(freqs[:N//2], np.abs(fft_value[:N//2]))
plt.yscale('log')
plt.ylabel("Amplitude")
plt.xlabel("Frequency [Hz]")
我得到以下情节
该图的 y 轴采用对数刻度。请注意，峰值的振幅以千为单位。其余大部分数据点的振幅都在100左右。
idx_max = np.argmax(np.abs(fft_value)) # Find the peak in the coefficients
idx_min = np.argmin(np.abs(fft_value)) # Find the peak in the coefficients
print(f'idx_max = {idx_max}, idx_min = {idx_min}')
print(f'f_max = {freqs[idx_max]}, f_min = {freqs[idx_min]}')
print(f'fft_value[idx_max] {fft_value[idx_max]}')
print(f'fft_value[idx_min] {fft_value[idx_min]}')
产生：
idx_max = 1546, idx_min = 1738
f_max = -10010.7, f_min = -5777.1
fft_value[idx_max] (-4733.232076236707+219.11718299533203j)
fft_value[idx_min] (-0.17017443966211232+0.9557200531465061j)
我正在将 link 添加到我构建的脚本中，该脚本输出具有实际幅度的 FFT（对于真实信号 - 例如您的信号）。试试看是否有效：
dt=1/frate在你的星座....
经过长时间的家庭作业，我终于找到了我的问题。正如我在 更新工作： 中提到的，原因是我采集的样本数量错误。
我改了两行代码
n_sa = 8 * int(freq_in_hertz)
t_fft = np.linspace(0, 1, n_sa)
到
n_sa = y.size //number of samples directly taken from the raw 16bits
t_fft = np.arange(n_sa)/frate //Here we need to divide each samples by the sampling rate
这解决了我的问题。
我的光谱输出是
特别感谢@meta4 和@YoniChechik 给了我一些建议。
我在这里使用 numpy 的 fft 函数来绘制从 10000Hz 正弦波生成的 PCM 波的 fft。但是我得到的情节幅度是错误的。
使用我在控制台本身打印的 fftfreq 函数，频率正在变得正确。我的 python 代码在这里。
import numpy as np
import matplotlib.pyplot as plt
frate = 44100
filename = 'Sine_10000Hz.bin' #signed16 bit PCM of a 10000Hz sine wave
f = open('Sine_10000Hz.bin','rb')
y = np.fromfile(f,dtype='int16') #Extract the signed 16 bit PCM value of 10000Hz Sine wave
f.close()
####### Spectral Analysis #########
fft_value = np.fft.fft(y)
freqs = np.fft.fftfreq(len(fft_value)) # frequencies associated with the coefficients:
print("freqs.min(), freqs.max()")
idx = np.argmax(np.abs(fft_value)) # Find the peak in the coefficients
freq = freqs[idx]
freq_in_hertz = abs(freq * frate)
print("\n\n\n\n\n\nfreq_in_hertz")
print(freq_in_hertz)
for i in range(2):
print("Value at index {}:\t{}".format(i, fft_value[i + 1]), "\nValue at index {}:\t{}".format(fft_value.size -1 - i, fft_value[-1 - i]))
#####
n_sa = 8 * int(freq_in_hertz)
t_fft = np.linspace(0, 1, n_sa)
T = t_fft[1] - t_fft[0] # sampling interval
N = n_sa #Here it is n_sample
print("\nN value=")
print(N)
# 1/T = frequency
f = np.linspace(0, 1 / T, N)
plt.ylabel("Amplitude")
plt.xlabel("Frequency [Hz]")
plt.xlim(0,15000)
# 2 / N is a normalization factor Here second half of the sequence gives us no new information that the half of the FFT sequence is the output we need.
plt.bar(f[:N // 2], np.abs(fft_value)[:N // 2] * 2 / N, width=15,color="red")
输出来自控制台（我在这里粘贴的只有最少的印刷品）
freqs.min(), freqs.max()
-0.5 0.49997732426303854
freq_in_hertz
10000.0
Value at index 0: (19.949569768991054-17.456031216294324j)
Value at index 44099: (19.949569768991157+17.45603121629439j)
Value at index 1: (9.216783424692835-13.477631008179145j)
Value at index 44098: (9.216783424692792+13.477631008179262j)
N value=
80000
频率提取正确，但在情节中我做的事情不正确，我不知道。
正在更新作品：
- 当我将
n_sa = 10 * int(freq_in_hertz)行中的乘数 10 更改为 5 时，我得到了正确的情节。对不对我看不懂 - 在
plt.xlim(0,15000)行中，如果我再次将最大值增加到 20000 则不会绘图。直到 15000 它都在正确绘制。 - 我使用 Audacity 工具生成了此
Sine_10000Hz.bin，在该工具中我生成了频率为 10000Hz、持续时间为 1 秒、采样率为 44100 的正弦波。然后我将此音频导出为带有无标头的签名 16 位（意味着原始 PCM） .我可以使用这个脚本重新生成这个正弦波。我也想计算这个的 FFT。所以我期望在 10000Hz 处有一个峰值，振幅为 32767。你可以看到我在行n_sa = 8 * int(freq_in_hertz)中更改了倍增因子8而不是10。因此它起作用了。但是振幅显示不正确。我会在这里附上我的新图
我不确定你到底想做什么，但我怀疑 Sine_10000Hz.bin 文件不是你想的那样。
是否可能包含多个频道（左&右）？
它真的是带符号的 16 位整数吗？
在 numpy 中创建 16 位整数的 10kHz 正弦波并不难。
import numpy as np
import matplotlib.pyplot as plt
n_samples = 2000
f_signal = 10000 # (Hz) Signal Frequency
f_sample = 44100 # (Hz) Sample Rate
amplitude = 2**3 # Arbitrary. Must be > 1. Should be > 2. Larger makes FFT results better
time = np.arange(n_samples) / f_sample # sample times
# The signal
y = (np.sin(time * f_signal * 2 * np.pi) * amplitude).astype('int16')
如果绘制 30 个信号点，您可以看到每个周期大约有 5 个点。
plt.plot(time[:30], y[:30], marker='o')
plt.xlabel('Time (s)')
plt.yticks([]); # Amplitude value is artificial. hide it
如果绘制来自 Sine_10000Hz.bin 的 30 个数据样本，它每个周期大约有 5 个点吗？
这是我尝试重新创建我理解的 FFT 工作。
fft_value = np.fft.fft(y) # compute the FFT
freqs = np.fft.fftfreq(len(fft_value)) * f_sample # frequencies for each FFT bin
N = len(y)
plt.plot(freqs[:N//2], np.abs(fft_value[:N//2]))
plt.yscale('log')
plt.ylabel("Amplitude")
plt.xlabel("Frequency [Hz]")
我得到以下情节
该图的 y 轴采用对数刻度。请注意，峰值的振幅以千为单位。其余大部分数据点的振幅都在100左右。
idx_max = np.argmax(np.abs(fft_value)) # Find the peak in the coefficients
idx_min = np.argmin(np.abs(fft_value)) # Find the peak in the coefficients
print(f'idx_max = {idx_max}, idx_min = {idx_min}')
print(f'f_max = {freqs[idx_max]}, f_min = {freqs[idx_min]}')
print(f'fft_value[idx_max] {fft_value[idx_max]}')
print(f'fft_value[idx_min] {fft_value[idx_min]}')
产生：
idx_max = 1546, idx_min = 1738
f_max = -10010.7, f_min = -5777.1
fft_value[idx_max] (-4733.232076236707+219.11718299533203j)
fft_value[idx_min] (-0.17017443966211232+0.9557200531465061j)
我正在将 link 添加到我构建的脚本中，该脚本输出具有实际幅度的 FFT（对于真实信号 - 例如您的信号）。试试看是否有效：
dt=1/frate在你的星座....
经过长时间的家庭作业，我终于找到了我的问题。正如我在 更新工作： 中提到的，原因是我采集的样本数量错误。
我改了两行代码
n_sa = 8 * int(freq_in_hertz)
t_fft = np.linspace(0, 1, n_sa)
到
n_sa = y.size //number of samples directly taken from the raw 16bits
t_fft = np.arange(n_sa)/frate //Here we need to divide each samples by the sampling rate
这解决了我的问题。
我的光谱输出是
特别感谢@meta4 和@YoniChechik 给了我一些建议。