原始 PCM 的 FFT 图因 python 中的更高频率而出错答案

【问题标题】：FFT plot of raw PCM comes wrong for higher frequency in python原始 PCM 的 FFT 图因 python 中的更高频率而出错
【发布时间】：2020-07-15 05:24:43
【问题描述】：

这里我使用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，在该工具中，我生成了一个持续时间为 1 秒、采样率为 44100 的频率为 10000Hz 的正弦波。然后我将此音频导出为无标头的签名 16 位（意味着原始 PCM）。我可以使用这个脚本重新生成这个正弦波。我也想计算这个的FFT。所以我预计在 10000Hz 处有一个峰值，振幅为 32767。你可以看到我在 n_sa = 8 * int(freq_in_hertz) 行中更改了乘法因子 8 而不是 10。因此它奏效了。但幅度显示不正确。我将在此处附上我的新图

【问题讨论】：

标签： python-3.x numpy matplotlib plot fft

【解决方案1】：

我不确定你到底想做什么，但我怀疑 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)

【讨论】：

>> 它是否可能包含多个通道（左和右）？不是只有一个频道。 >>真的是有符号的 16 位整数吗？是的
嗨@meta4：我认为我的问题在于选择 'n_sa' 、 'T" 、 'f' 、 't_fft'
@Vishnu CS，您是否尝试绘制 30 个原始（非 fft-ed）数据样本？你能看到“周期性”签名吗？

【解决方案2】：

我正在向我构建的脚本添加一个链接，该脚本输出具有实际幅度的 FFT（对于真实信号 - 例如您的信号）。试试看它是否有效：

dt=1/frate 在你的星座中......

https://stackoverflow.com/a/53925342/4879610

【讨论】：

好的..我会检查它是否解决了我的问题
实际上有两个问题。我回答了我的问题。非常感谢

【解决方案3】：

经过长时间的家庭作业，我能够找到我的问题。正如我在更新作品中提到的：原因是我采集的样本数量错误。

我修改了代码中的两行

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 给我一些建议。

【讨论】：