График БПФ необработанного PCM неверен для более высокой частоты в python

Question

Здесь я использую fft функцию numpy, чтобы построить FFT волны PCM, генерируемой из синусоидальной волны 10000 Гц. Но амплитуда графика, который я получаю, неверна.

Частота получается правильной с использованием функции 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

Извлечение частоты идет правильно, но на графике что-то, что я делаю, неверно, чего я не знаю.

Обновление работы:

1. Когда я меняю коэффициент умножения 10 в строке n_sa = 10 * int(freq_in_hertz) до 5 дает мне правильный сюжет. Правильно ли это или нет, я не могу понять,
2. В строке plt.xlim(0,15000), если я увеличу максимальное значение до 20000 , снова не построение. До 15000 оно строится правильно.
3. Я сгенерировал это Sine_10000Hz.bin, используя инструмент Audacity, где я генерирую синусоидальную волну с частотой 10000 Гц длительностью 1 с c и частотой дискретизации 44100. Затем я экспортировал этот звук в 16-битный бит со знаком без заголовка (означает, что PCM). Я мог бы восстановить эту синусоидальную волну, используя этот скрипт. Также я хочу рассчитать БПФ этого. Поэтому я ожидаю пика на частоте 10000 Гц с амплитудой 32767. Вы можете видеть, что я изменил коэффициент умножения 8 вместо 10 в строке n_sa = 8 * int(freq_in_hertz). Следовательно, это сработало. Но амплитуда показывает неверно. Я приложу мою новую фигуру здесь

Answer 1

Я не совсем уверен, что именно вы пытаетесь сделать, но я подозреваю, что файл Sine_10000Hz.bin - это не то, что вы думаете.

Возможно, он содержит более одного канала (слева и справа)?

Реально ли подписаны 16-битные целые числа?

Нетрудно создать синусоидальную волну 10 кГц в 16-битных целых числах в numpy.

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

Если вы строите 30 образцы данных из Sine_10000Hz.bin имеет ли он около 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)

Answer 2

После долгой домашней работы я смог найти свою проблему. Как я упоминал в Обновлении работы: причина была в том, что количество образцов, которые я взял, было неверным.

Я изменил две строки в коде

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 за предоставленные мне советы.

Answer 3

Я добавляю ссылку на созданный мной скрипт, который выводит БПФ с АКТУАЛЬНОЙ амплитудой (для реальных сигналов - например, вашего сигнала). Имейте go и посмотрите, работает ли он:

dt=1/frate в вашем созвездии ....

{ ссылка }