不確定如何將FFT數據用於頻譜分析儀

Question

我正在嘗試創建一個帶8個LED條的自制頻譜分析儀。

我正在努力的部分是執行FFT並瞭解如何使用結果。

到目前爲止，這是我：

import opc 
import time 
import pyaudio 
import wave 
import sys 
import numpy 
import math 

CHUNK = 1024 

# Gets the pitch from the audio 
def pitch(signal): 
    # NOT SURE IF ANY OF THIS IS CORRECT 
    signal = numpy.fromstring(signal, 'Int16'); 
    print "signal = ", signal 

    testing = numpy.fft.fft(signal) 
    print "testing = ", testing 

wf = wave.open(sys.argv[1], 'rb') 
RATE = wf.getframerate() 
p = pyaudio.PyAudio() # Instantiate PyAudio 

# Open Stream 
stream = p.open(format=p.get_format_from_width(wf.getsampwidth()), 
       channels=wf.getnchannels(), 
       rate=wf.getframerate(), 
       output=True) 

# Read data 
data = wf.readframes(CHUNK) 

# Play Stream 
while data != '': 
    stream.write(data) 
    data = wf.readframes(CHUNK) 
    frequency = pitch(data) 
    print "%f frequency" %frequency 

我與在pitch方法做什麼掙扎。我知道我需要對傳入的數據執行FFT，但我真的不確定如何執行此操作。

也應該用this函數？

Answer 1

由於np.fft.fft的工作原理，如果使用1024個數據點，您將獲得512個頻率的值（加零值Hz，DC偏移量）。如果您只需要8個頻率，則必須爲其提供16個數據點。

您可能可以按照64倍因子進行降採樣 - 然後16個降採樣點將是時間等效到1024個原始點。我從來沒有探索過這個，所以我不知道這會帶來什麼或者什麼是缺陷。

你將不得不做一些學習 - The Scientist and Engineer's Guide to Digital Signal Processing真的是一個優秀的資源，至少它是我的。

請記住，對於音頻cd .wav文件，採樣頻率爲44100 Hz - 1024樣本塊只有23 mS的聲音。

scipy.io.wavfile.read使獲取數據變得簡單。

samp_rate, data = scipy.io.wavfile.read(filename) 

data是2- d numpy的陣列與一個信道在零列，數據[：，0]，和其他在第1列，數據[：，1]

Matplotlib的specgram和psd函數可以給你你想要的數據。一個圖表模擬您正在嘗試做什麼。

from matplotlib import pyplot as plt 
import scipy.io.wavfile 
samp_rate, data = scipy.io.wavfile.read(filename) 
Pxx, freqs, bins, im = plt.specgram(data[:1024,0], NFFT = 16, noverlap = 0, Fs = samp_rate) 
plt.show() 
plt.close() 

由於您沒有進行任何繪圖，只是使用matplolib.mlab.specgram。

Pxx, freqs, t = matplolib.mlab.specgram(data[:1024,0], NFFT = 16, noverlap = 0, Fs = samp_rate) 

它的返回值（地址Pxx，freqs，噸）是

 - *Pxx*: 2-D array, columns are the periodograms of successive segments 

    - *freqs*: 1-D array of frequencies corresponding to the rows in Pxx 

    - *t*: 1-D array of times corresponding to midpoints of segments. 

Pxx[1:, 0]將是的值的頻率爲T0，Pxx[1:, 1]爲T1，Pxx[1:, 2]爲T2， ...這就是你要顯示的東西。您不使用Pxx[0, :]，因爲它用於0 Hz。

功率譜密度 - matplotlib.mlab.psd()

---

也許另一個戰略，以獲得降至8個帶是使用大塊和規範值。然後你可以將這些值分成8個部分並得到每個部分的總和。我認爲這是有效的 - 可能只是針對功率譜密度。 sklearn.preprocessing.normalize

w = sklearn.preprocessing.normalize(Pxx[1:,:], norm = 'l1', axis = 0) 

---

不過話又說回來，我只是做一切了。

Answer 2

我不知道@wwii在他的回答中提到的scipy.io.wavfile.read函數，但似乎他的建議是處理信號加載的方法。不過，我只想評論一下傅立葉變換。

我想你打算用你的LED設置做什麼，是根據你打算使用的8個頻段中每一個頻段的光譜功率來改變每個LED的亮度。因此，我理解你需要的是隨着時間的推移以某種方式計算權力。第一個複雜因素是「如何計算光譜功率？」

要做到這一點的最佳方法是使用numpy.fft.rfft計算只有實數（不是複數）的信號的傅里葉變換。另一方面，函數numpy.fft.fft是一個通用函數，可以計算具有複數的信號的快速傅立葉變換。概念上的區別是numpy.fft.fft可用於研究行波及其傳播方向。這是因爲返回的振幅對應於positive or negative frequencies，表明波如何傳播。 numpy.fft.rfft會產生實數頻率的振幅，如numpy.fft.rfftfreq所示，這正是您所需要的。

最後一個問題是選擇合適的頻段來計算頻譜功率。人耳具有巨大的頻率響應範圍，每個頻帶的寬度將變化很大，低頻帶非常窄，高頻帶非常寬。谷歌搜索周圍，我發現this好的資源，它定義7個相關頻帶

1. 子低音：20至60赫茲
2. 低音：60至250Hz
3. 低中檔：250至500Hz
4. 中型：500Hz至2kHz的
5. 上中音：2至4kHz
6. 存在：4至6千赫
7. 光輝：6至20kHz

我會建議使用這些頻段，但將高頻中頻分成2-3 kHz和3-4 kHz。這樣你就可以使用8個LED設置。我上傳的更新的間距的功能讓你使用

wf = wave.open(sys.argv[1], 'rb') 
CHUNK = 1024 
RATE = wf.getframerate() 
DT = 1./float(RATE) # time between two successive audio frames 
FFT_FREQS = numpy.fft.nfftfreq(CHUNCK,DT) 
FFT_FREQS_INDS = -numpy.ones_like(FFT_FREQS) 
bands_bounds = [[20,60],  # Sub-bass 
       [60,250],  # Bass 
       [250,500], # Low midrange 
       [500,2000], # Midrange 
       [2000,3000], # Upper midrange 0 
       [3000,4000], # Upper midrange 1 
       [4000,6000], # Presence 
       [6000,20000]] # Brilliance 

for f_ind,freq in enumerate(FFT_FREQS): 
    for led_ind,bounds in enumerate(bands_bounds): 
     if freq<bounds[1] and freq>=bounds[0]: 
      FFT_FREQS_INDS[ind] = led_ind 

# Returns the spectral power in each of the 8 bands assigned to the LEDs 
def pitch(signal): 
    # CONSIDER SWITCHING TO scipy.io.wavfile.read TO GET SIGNAL 
    signal = numpy.fromstring(signal, 'Int16'); 
    amplitude = numpy.fft.rfft(signal.astype(numpy.float)) 
    power = [np.sum(np.abs(amplitude[FFT_FREQS_INDS==led_ind])**2) for led_ind in range(len(bands_bounds))] 
    return power 

代碼的第一部分計算FFT頻率和構建體，用於指示所述8個頻帶的FFT頻率對應於所述陣列FFT_FREQS_INDS。然後，在pitch中計算每個頻帶中頻譜的功率。當然，這可以優化，但我試圖讓代碼不言自明。