Короткое преобразование Фурье в питоне

Question

Я хочу получить частоту с максимальной мощностью для каждого момента в wav-файле. Итак, я написал STFT в Python, используя fft из scipy. Я использовал функцию окна кайзера из scipy. Все выглядит великолепно, но мой результат выглядит странно. У этого есть некоторые очень маленькие числа и некоторые очень высокие.

здесь есть выход на один WAV файл: http://pastebin.com/5Ryd2uXj и вот код в Python:

import scipy, pylab 
import wave 
import struct 
import sys 

def stft(data, cp, do, hop): 
    dos = int(do*cp) 
    w = scipy.kaiser(dos,12) //12 is very high for kaiser window 
    temp=[] 
    wyn=[] 
    for i in range(0, len(data)-dos, hop): 
     temp=scipy.fft(w*data[i:i+dos]) 
     max=-1 
     for j in range(0, len(temp),1): 
      licz=temp[j].real**2+temp[j].imag**2 
      if(licz>max): 
       max = licz 
       maxj = j 
     wyn.append(maxj) 
    #wyn = scipy.array([scipy.fft(w*data[i:i+dos]) 
     #for i in range(0, len(data)-dos, 1)]) 
    return wyn 

file = wave.open(sys.argv[1]) 
bity = file.readframes(file.getnframes()) 
data=struct.unpack('{n}h'.format(n=file.getnframes()), bity) 
file.close() 

cp=44100 #sampling frequency 
do=0.05 #window size 
hop = 5 

wyn=stft(data,cp,do,hop) 
print len(wyn) 
for i in range(0, len(wyn), 1): 
    print wyn[i] 

Answer 1

Действительное FT синусоидальной волны пара дельта-функций, равноудаленных от 0-частоты. С дискретной функцией (выборки) это повторяется каждые fs (частота дискретизации) в частотной области. Небольшие ошибки при вычислении FFT будут означать, что эти две дельта (FT вашей синусоидальной волны) не будут точно такой же высоты, поэтому ваш алгоритм просто выбирает более высокий.

Функция Scipy FFT даст вам частотные составляющие с доменом [0, fs]. Поскольку (как я уже упоминал выше) это периодическое значение, эти значения также можно переустановить как [-fs/2, fs/2] путем замены результата в центральной точке. Для этого воспользуйтесь fftshift. Похоже, вас интересует только положительных частот, поэтому вы можете просто отбросить вторую половину результата вашего БПФ.

Из записок scipy.fftpack.fft:

Упаковка результата является «стандартным»: Если A = FFT (а, п), то А [0] содержит термин нулевой частоты, A [ 1: n/2 + 1] содержит положительно-частотные члены, а A [n/2 + 1:] содержит отрицательно-частотные члены в порядке убывающей частоты. Таким образом, для 8-точечного преобразования частоты результата равны [0, 1, 2, 3, 4, -3, -2, -1].