БПФ данных, полученных от PyAudio, дает неправильную частоту - PullRequest
Новая
       58

БПФ данных, полученных от PyAudio, дает неправильную частоту

0 голосов
/ 08 февраля 2019

Моя главная задача - распознавать гул человека в микрофоне в режиме реального времени.В качестве первого шага к общему распознаванию сигналов я сделал 5-секундную запись сигнала 440 Гц, генерируемого приложением на моем телефоне, и попытался обнаружить ту же частоту.

Я использовал Audacity для построения и проверкиспектр из того же wav-файла 440 Гц, и я получил это, которое показывает, что 440 Гц действительно является доминирующей частотой: (* ​​1003 *https://i.imgur.com/2UImEkR.png)

Чтобы сделать это с python, я использую библиотеку PyAudio и ссылаюсь на этот блог . У меня есть код, который я запускаю с файлом wav: 

"""PyAudio Example: Play a WAVE file."""

import pyaudio
import wave
import sys
import struct
import numpy as np
import matplotlib.pyplot as plt

CHUNK = 1024

if len(sys.argv) < 2:
    print("Plays a wave file.\n\nUsage: %s filename.wav" % sys.argv[0])
    sys.exit(-1)

wf = wave.open(sys.argv[1], 'rb')

p = pyaudio.PyAudio()
stream = p.open(format=p.get_format_from_width(wf.getsampwidth()),
                channels=wf.getnchannels(),
                rate=wf.getframerate(),
                output=True)

data = wf.readframes(CHUNK)

i = 0
while data != '':
    i += 1
    data_unpacked = struct.unpack('{n}h'.format(n= len(data)/2 ), data) 
    data_np = np.array(data_unpacked) 
    data_fft = np.fft.fft(data_np)
    data_freq = np.abs(data_fft)/len(data_fft) # Dividing by length to normalize the amplitude as per https://www.mathworks.com/matlabcentral/answers/162846-amplitude-of-signal-after-fft-operation
    print("Chunk: {} max_freq: {}".format(i,np.argmax(data_freq)))

    fig = plt.figure()
    ax = fig.add_subplot(1,1,1)
    ax.plot(data_freq)
    ax.set_xscale('log')
    plt.show()

    stream.write(data)
    data = wf.readframes(CHUNK)

stream.stop_stream()
stream.close()

p.terminate()

В выводе я получаю, что максимальная частота равна 10 для всехкуски и пример одного из графиков: (https://i.imgur.com/zsAXME5.png)

Я ожидал, что это значение будет 440 вместо 10 для всех кусков. Я признаю, что очень мало знаю о теории БПФ иЯ ценю любую помощь, позволившую мне решить эту проблему.

РЕДАКТИРОВАТЬ: частота дискретизации составляет 44100. количество каналов равно 2 и ширина выборки также 2.

1 Ответ

0 голосов
/ 08 февраля 2019

Предисловия

Как указывалось xdurch0, вы читаете некий индекс, а не частоту.Если вы собираетесь делать все вычисления самостоятельно, вам нужно вычислить собственный частотный вектор перед построением графика, если вы хотите получить согласованный результат.Чтение этого ответа может помочь вам найти решение.

Вектор частоты для БПФ (полуплоскость): 

 f = np.linspace(0, rate/2, N_fft/2)

Или (полная плоскость): 

 f = np.linspace(-rate/2, rate/2, N_fft)

С другой стороны, мы можем делегировать большую часть работы превосходному набору инструментов scipy.signal, который нацелен на решение подобных проблем (и многих других).

MCVE

Используя пакет scipy, можно получить желаемый результат для простого файла WAV с одной частотой ( source ): 

import numpy as np
from scipy import signal
from scipy.io import wavfile
import matplotlib.pyplot as plt

# Read the file (rate and data):
rate, data = wavfile.read('tone.wav') # See source

# Compute PSD:
f, P = signal.periodogram(data, rate) # Frequencies and PSD

# Display PSD:
fig, axe = plt.subplots()
axe.semilogy(f, P)
axe.set_xlim([0,500])
axe.set_ylim([1e-8, 1e10])
axe.set_xlabel(r'Frequency, $\nu$ $[\mathrm{Hz}]$')
axe.set_ylabel(r'PSD, $P$ $[\mathrm{AU^2Hz}^{-1}]$')
axe.set_title('Periodogram')
axe.grid(which='both')

В основном:

. Вывод:

enter image description here

Find Peak

Тогда мы можем найти частоту первого наивысшего пика (P>1e-2, этот критерий подлежит настройке), используя find_peaks: 

idx = signal.find_peaks(P, height=1e-2)[0][0]
f[idx] # 440.0 Hz

Вводвсе вместе это сводится к следующему: 

def freq(filename, setup={'height': 1e-2}):
    rate, data = wavfile.read(filename)
    f, P = signal.periodogram(data, rate)
    return f[signal.find_peaks(P, **setup)[0][0]]

Обработка нескольких каналов

Я попробовал этот код с моим WAV-файлом, и получил ошибку для строки axe.semilogy (f, Pxx_den) следующим образом: ValueError: x и y должны иметь одинаковое первое измерение.Я проверил формы и f имеет (2,), в то время как Pxx_den имеет (220160,2).Кроме того, кажется, что массив Pxx_den имеет все нули.

WAV-файл может содержать несколько каналов, в основном это монофонические или стереофайлы (макс. 2**16 - 1 каналы).Проблема, которую вы подчеркнули, возникает из-за файла с несколькими каналами ( стерео семпл ). 

rate, data = wavfile.read('aaaah.wav') # Shape: (46447, 2), Rate: 48 kHz

enter image description here

Этоне очень хорошо документировано , но метод signal.periodogram также работает с матрицей, и его ввод не согласуется напрямую с выводом wavfile.read (по умолчанию они работают на разных осях).Поэтому нам нужно тщательно ориентировать размеры (используя переключатель axis) при выполнении PSD: 

f, P = signal.periodogram(data, rate, axis=0, detrend='linear')

Он также работает с транспозицией data.T, но затем нам нужно транспонировать результат обратно.

Задание оси решает проблему: частотный вектор является правильным, и PSD не везде равен нулю (до того, как он выполнялся на axis=1 длиной 2, в вашем случае он выполнил 220160 PSD для сигналов с 2 выборками, которые мы хотелинаоборот).

Переключатель detrend гарантирует, что сигнал имеет нулевое среднее значение, а его линейный тренд удален.

Реальное приложение

Этот подход должен работать для реальных фрагментированных выборок,при условии, что фрагменты содержат достаточно данных (см. Теорема выборки Найквиста-Шеннона ).Затем данные представляют собой подвыборки сигнала (куски), и скорость поддерживается постоянной, поскольку она не изменяется во время процесса.

Имея куски размером 2**10, кажется, работает, мы можем определить конкретные частоты из них: 

f, P = signal.periodogram(data[:2**10,:], rate, axis=0, detrend='linear') # Shapes: (513,) (513, 2)
idx0 = signal.find_peaks(P[:,0], threshold=0.01, distance=50)[0] # Peaks: [46.875, 2625., 13312.5, 16921.875] Hz

fig, axe = plt.subplots(2, 1, sharex=True, sharey=True)
axe[0].loglog(f, P[:,0])
axe[0].loglog(f[idx0], P[idx0,0], '.')
# [...]

enter image description here

На данный момент самая сложная часть - это точная настройка метода find-peaks для захвата желаемых частот.Возможно, вам придется подумать о предварительной фильтрации вашего сигнала или последующей обработке PSD, чтобы упростить идентификацию.

...