不确定如何将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中计算每个频带中频谱的功率。当然，这可以优化，但我试图让代码不言自明。