Как получить время/частоту из FFT в Python

Question

У меня есть небольшая проблема с управлением данными FFT. Я искал много примеров того, как делать БПФ, но я не мог получить то, что хочу от любого из них. У меня есть случайный волновой файл с частотой дискретизации 44 кГц, и я хочу получить величину N гармоник каждый X мс, скажем, 100 мс должно быть достаточно. Я пробовал этот код:

import scipy.io.wavfile as wavfile 
import numpy as np 
import pylab as pl 

rate, data = wavfile.read("sound.wav") 
t = np.arange(len(data[:,0]))*1.0/rate 
p = 20*np.log10(np.abs(np.fft.rfft(data[:2048, 0]))) 
f = np.linspace(0, rate/2.0, len(p)) 
pl.plot(f, p) 
pl.xlabel("Frequency(Hz)") 
pl.ylabel("Power(dB)") 
pl.show() 

Это был последний пример, который я использовал, я нашел его где-то в stackoverflow. Проблема в том, что это получает величину, которую я хочу, получает частоту, но времени нет. Насколько я знаю, анализ БПФ является 3D, и это «объединенный» результат всех гармоник. Я получаю это:

X-axis = Frequency, Y-axis = Magnitude, Z-axis = Time (invisible)

Из моего понимания кода, т время - и, похоже, что, но не требуется в коде - Мы, возможно, нужно хотя. p - массив мощностей (или величины), но он кажется некоторым средним из всех величин каждой частоты f, который представляет собой массив частот. Я не хочу среднего/объединенного значения, мне нужна величина для N гармоник каждые X миллисекунд.

Короче говоря, мы можем получить: 1 магнитуду всех частот.

Мы хотим: все величины N частот, включая время, когда присутствует определенная величина.

Результат должен выглядеть этот массив: [время, частота, амплитуда] Таким образом, в конце концов, если мы хотим, 3 гармоники, она будет выглядеть так:

[0,100,2.85489] #100Hz harmonic has 2.85489 amplitude on 0ms 
[0,200,1.15695] #200Hz ... 
[0,300,3.12215] 
[100,100,1.22248] #100Hz harmonic has 1.22248 amplitude on 100ms 
[100,200,1.58758] 
[100,300,2.57578] 
[200,100,5.16574] 
[200,200,3.15267] 
[200,300,0.89987] 

Визуализация не требуется, результат должен быть только массивы (или хеши/словари), как указано выше.

Answer 1

В дополнение к ответу @Paul АиР, scipy.signal.spectrogram в scipy's signal processing modulespectrogram function.

Пример по ссылке выше выглядит следующим образом:

from scipy import signal 
import matplotlib.pyplot as plt 

# Generate a test signal, a 2 Vrms sine wave whose frequency linearly 
# changes with time from 1kHz to 2kHz, corrupted by 0.001 V**2/Hz of 
# white noise sampled at 10 kHz. 

fs = 10e3 
N = 1e5 
amp = 2 * np.sqrt(2) 
noise_power = 0.001 * fs/2 
time = np.arange(N)/fs 
freq = np.linspace(1e3, 2e3, N) 
x = amp * np.sin(2*np.pi*freq*time) 
x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape) 


#Compute and plot the spectrogram. 

f, t, Sxx = signal.spectrogram(x, fs) 
plt.pcolormesh(t, f, Sxx) 
plt.ylabel('Frequency [Hz]') 
plt.xlabel('Time [sec]') 
plt.show() 

Answer 2

Похоже, вы пытаетесь реализовать spectrogram, который представляет собой последовательность оценок спектра мощности, обычно реализуемую с чередованием (обычно перекрывающихся) БПФ. Поскольку у вас только один БПФ (спектр), то у вас еще нет измерения времени. Поместите свой код FFT в цикл и обработайте один блок выборок (например, 1024) на итерацию с 50% перекрытием между последовательными блоками. Последовательность сгенерированных спектров будет представлять собой трехмерную матрицу времени v с частотой v.

Я не Python человек, но я могу дать вам некоторые псевдо-код, который должен быть достаточно, чтобы вы кодирования:

N = length of data input 
N_FFT = no of samples per block (== FFT size, e.g. 1024) 
i = 0 ;; i = index of spectrum within 3D output array 
for block_start = 0 to N - block_start 
    block_end = block_start + N_FFT 
    get samples from block_start .. block_end 
    apply window function to block (e.g. Hamming) 
    apply FFT to windowed block 
    calculate magnitude spectrum (20 * log10(re*re + im*im)) 
    store spectrum in output array at index i 
    block_start += N_FFT/2   ;; NB: 50% overlap 
    i++ 
end 

Answer 3

Edit: Ах, так, кажется, это возвращает значения, но они не подходят к звуковому файлу в все. Несмотря на то, что они могут использоваться как амплитуда на спектрограмме, они не будут работать, например, в тех классических аудио-визуализаторах, которые вы можете видеть во многих музыкальных проигрывателях. Я также попробовал pylab matplotlib для спектрограммы, но результат тот же.

import os 
import wave 
import pylab 
import math 
from numpy import amax 
from numpy import amin 

def get_wav_info(wav_file,mi,mx): 
    wav = wave.open(wav_file, 'r') 
    frames = wav.readframes(-1) 
    sound_info = pylab.fromstring(frames, 'Int16') 
    frame_rate = wav.getframerate() 
    wav.close() 
    spectrum, freqs, t, im = pylab.specgram(sound_info, NFFT=1024, Fs=frame_rate) 
    n = 0 
    while n < 20: 
     for index,power in enumerate(spectrum[n]): 
      print("%s,%s,%s" % (n,int(round(t[index]*1000)),math.ceil(power*100)/100)) 
     n += 1 

get_wav_info("wave.wav",1,20) 

Любые советы по получению дБ, которые можно использовать при визуализации? В принципе, мы, видимо, имеем все, что нам нужно, из кода выше, как заставить его возвращать нормальные значения?Игнорируйте mi и mx, так как они просто регулируют значения в массиве, чтобы вписаться в интервал mi..mx - это было бы для использования визуализации. Если я прав, spectrum в этом коде возвращает массив массивов, который содержит амплитуды для каждой частоты от freqs массива, которые присутствуют вовремя в соответствии с массивом t, но как работает значение - это действительно амплитуда, если она возвращает эти странные значения и если это так, как преобразовать его в дБ, например.

tl; dr Мне нужен вывод для визуализатора, как у музыкальных плееров, но он не должен работать в реальном времени, я хочу только данные, но значения не соответствуют файлу wav.

Редактировать 2: Я заметил, что есть еще одна проблема. За 90 секунд wav, t массив содержит время до 175.x, что кажется очень странным, учитывая, что frame_rate указан в файле wav. Итак, теперь у нас есть 2 проблемы: spectrum, похоже, не возвращает правильные значения (возможно, он будет соответствовать, если мы получим правильное время), и t, кажется, возвращает ровно двойное время wav.

Исправлено: Корпус полностью решен.

import os 
import pylab 
import math 
from numpy import amax 
from numpy import amin 
from scipy.io import wavfile 
frame_rate, snd = wavfile.read(wav_file) 
sound_info = snd[:,0] 
spectrum, freqs, t, im = pylab.specgram(sound_info,NFFT=1024,Fs=frame_rate,noverlap=5,mode='magnitude') 

Specgram нуждался в небольшой настройке, и я загрузил только один канал с библиотекой scipy.io (вместо волновой библиотеки). Также без режима, установленного на величину, он возвращает 10log10 вместо 20log10, что является причиной того, что он не вернул правильные значения.