FFT意想不到的頻移

我得到的聲音信號的FFT計算這個Python代碼：FFT意想不到的頻移

from math import * 
from cmath import exp, pi 

def fft(x): 
    N = len(x) 
    if N <= 1: return x 
    even = fft(x[0::2]) 
    odd = fft(x[1::2]) 
    return ([even[k] + exp(-2j * pi * k/N) * odd[k] for k in xrange(N/2)] + 
      [even[k] - exp(-2j * pi * k/N) * odd[k] for k in xrange(N/2)]) 

N = 64 

res = [sin(k) for k in xrange(N)] 

# Window function 
a = 2*pi/(N-1) 
for k in xrange(N): 
    z = a * res[k] 
    res[k] = 0.42659 - 0.49656*cos(z) + 0.076849*cos(2*z) 


res = fft(res) 

for k in xrange(N/2): 
    # get the amplitude... 
    sqr = sqrt(res[k].real * res[k].real + res[k].imag * res[k].imag) 
    if sqr > 0: 
     print 20 * log10(sqr) # ...in decibels 
    else: 
     print "-INF"

我得到了以下的結果：

WITHOUT窗函數（註釋的）：

-20.3017238269 
-16.9192604108 
-12.5089302395 
-8.97999530657 
-5.96033201086 
-3.12975820108 
-0.242090896634 
2.97021879504 
6.95134203457 
12.8752188937 
29.5096108632 <-- PEAK 
17.1668404562 
10.6485650284 
7.24321329787 
4.98448122464 
3.3242079033 
2.03154022635 
0.987966110459 
0.124898554197 
-0.600705355004 
-1.21748302238 
-1.74534177237 
-2.1985940834 
-2.5878009699 
-2.92091067118 
-3.20399051424 
-3.44171254421 
-3.63768393032 
-3.79467588076 
-3.91478386211 
-3.99953964778 
-4.04998822971

有窗函數

：

-6.55943077129 
-65.8567720414 
-65.7987645827 
-65.7012678903 
-65.5625673034 
-65.380788761 
-65.1529344157 
-64.8750852394 
-64.5420675211 
-64.1470597764 
-63.6810080181 
-63.131731575 
-62.4825087571 
-61.7097419947 
-60.7788888801 
-59.6368610687 
-58.1964601495 
-56.3001921054 
-53.6185951634 
-49.2554491173 
-38.3322646561 <-- PEAK 
-43.3318138698 
-52.0838904305 
-56.7277347745 
-60.2038755771 
-62.9772322874 
-65.442363488 
-67.7550361967 
-70.0212827894 
-72.3056579688 
-74.5822818952 
-76.5522909937

由於某種原因，峯值出現偏移。這是一個2倍頻移！我試過這個Java applet：
http://www.random-science-tools.com/maths/FFT.htm

看來，沒有任何窗口函數的結果是正確的（譜圖的1/3的峯值）。相反，如果我在我的python腳本中應用窗口函數，峯值將顯示在頻譜的2/3處。

這是應該發生的嗎？我如何解決它？

來源

2013-10-04 user2464424

好的，同時我意識到什麼是錯的。我在問題中寫的窗口函數完全沒有意義。

這是正確的：

a = 2*pi/(N-1) 

for k in xrange(N): 
    z = a * k 
    res[k] *= 0.42659 - 0.49656*cos(z) + 0.076849*cos(2*z) # Blackman

結果：

-63.8888312044 
-62.1859660802 
-59.4560808775 
-57.5235455007 
-57.0010514385 
-59.4284419437 
-66.6535724743 
-46.1441434426 
-2.31562840406 
16.0873761957 
22.4136439765 <-- PEAK 
19.5784749467 
6.43274013629 
-28.3842042716 
-55.5273291654 
-68.8982705127 
-53.3843989911 
-49.731974213 
-48.3131204305 
-47.6953570892 
-47.4386151256 
-47.361972079 
-47.3787962267 
-47.4434419084 
-47.530228024 
-47.6240076874 
-47.7155325706 
-47.799012933 
-47.870764286 
-47.9284264139 
-47.9705003855 
-47.9960714351

峯現在也正是它應該是。

你可能想嘗試一些其他的窗口：

res[k] *= 0.355768 - 0.487396*cos(z) + 0.144232*cos(2*z) - 0.012604*cos(3*z) 
res[k] *= 1 - 1.93*cos(z) + 1.29*cos(2*z) - 0.388*cos(3*z) + 0.028*cos(4*z) 
res[k] *= 1 - 1.985844164102*cos(z) + 1.791176438506*cos(2*z) - 1.282075284005*cos(3*z) + 0.667777530266*cos(4*z) - 0.240160796576*cos(5*z) + 0.056656381764*cos(6*z) - 0.008134974479*cos(7*z) + 0.000624544650*cos(8*z) - 0.000019808998*cos(9*z) + 0.000000132974*cos(10*z)

爲了：納托爾，FTSRS，HFT248D。

https://holometer.fnal.gov/GH_FFT.pdf

來源

2013-10-04 19:36:15 user2464424

FFT意想不到的頻移

回答

相關問題