python中的短时傅里叶变换答案

【问题标题】：Short Time Fourier Transform in pythonpython中的短时傅里叶变换
【发布时间】：2012-03-13 01:31:26
【问题描述】：

我想在 wav 文件中的每一刻都获得最大功率的频率。所以我使用来自 scipy 的 fft 在 Python 中编写了 STFT。我使用了 scipy 的 kaiser 窗口函数。一切看起来都很棒，但我的输出看起来很奇怪。它有一些非常小的数字和一些非常高的数字。

这是一个 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]

【问题讨论】：

您是否尝试针对已知波形（如正弦波）对其进行测试，看看您是否获得了预期的输出？
我刚刚发现了这个：stackoverflow.com/questions/2459295/stft-and-istft-in-python 它看起来很相似，我看到在窦图中是 2 条线，而不是 1 条。我的窦输出中有相同的。我不知道为什么...

标签： python fft

【解决方案1】：

正弦波的实际 FT 是一对与 0 频率等距的 delta 函数。对于离散函数（样本），在频域中每隔fs（采样率）重复一次。 FFT 计算中的小错误意味着这两个增量（正弦波的 FT）的高度不会完全相同，因此您的算法只是选择较高的一个。

scipy FFT 函数将为您提供域[0, fs] 的频率分量。由于（正如我上面提到的）这是周期性的，因此这些值也可以通过在中心点交换结果来重新映射为 [-fs/2, fs/2] - 考虑使用 fftshift 来执行此操作。但是，听起来您可能只对正频率感兴趣，因此您可以简单地丢弃 FFT 结果的后半部分。

来自scipy.fftpack.fft的笔记：

结果的打包是“标准的”：如果 A = fft(a, n)，则 A[0] 包含零频率项，A[1:n/2+1] 包含正频率项项，并且 A[n/2+1:] 包含负频率项，按负频率递减的顺序排列。所以对于 8 点变换，结果的频率是 [0, 1, 2, 3, 4, -3, -2, -1]。

【讨论】：