2009-05-31 75 views
2

音頻處理對我來說是非常新穎的。目前正在使用Python Numpy來處理wave文件。計算FFT矩陣後,我得到了不存在頻率的噪聲功率值。我對可視化數據感興趣,準確性不是高優先級。有沒有一種安全的方法來計算削波值以去除這些值,還是應該使用每個樣本集的所有FFT矩陣來計算平均值?剪切FFT矩陣

問候

編輯:

from numpy import * 
    import wave 
    import pymedia.audio.sound as sound 
    import time, struct 
    from pylab import ion, plot, draw, show 

    fp = wave.open("500-200f.wav", "rb") 
    sample_rate = fp.getframerate() 
    total_num_samps = fp.getnframes() 
    fft_length = 2048. 
    num_fft = (total_num_samps/fft_length) - 2 
    temp = zeros((num_fft,fft_length), float) 

    for i in range(num_fft): 
     tempb = fp.readframes(fft_length); 
     data = struct.unpack("%dH"%(fft_length), tempb) 
     temp[i,:] = array(data, short) 
    pts = fft_length/2+1 
    data = (abs(fft.rfft(temp, fft_length))/(pts))[:pts] 

    x_axis = arange(pts)*sample_rate*.5/pts 
    spec_range = pts 
    plot(x_axis, data[0]) 
    show() 

以下是在非對數標度的曲線圖,對於包含500HZ(淡出)+ 200Hz的使用GOLDWAVE正弦波創建合成的波形文件。

+3

您是否使用已知良好的FFT輸出驗證了您的輸出? (matlab或fftw將是很好的來源)。另外,嘗試輸入純色調即以已知頻率出現的正弦波,並驗證您的輸出是否有不同的fft尺寸。 – basszero 2009-06-01 00:18:56

回答

3

模擬波形不應該顯示像你的數字的FFT,所以有些東西是非常錯誤的,可能不是FFT,而是輸入波形。你的情節中的主要問題不是漣漪,而是1000赫茲左右的諧波和500赫茲的次諧波。模擬的波形不應該顯示任何這個(例如,見下面我的情節)。

首先,您可能只想嘗試繪製原始波形,這可能會指出一個明顯的問題。另外,將波解包成無符號短路看起來很奇怪,即「H」,特別是在此之後沒有大的零頻分量。

正如分諧波和高次諧波(以及特雷弗)所建議的那樣,通過對波形應用削波,我能夠得到與FFT相當接近的副本。您可以在模擬或拆包中引入裁剪。無論哪種方式,我通過在numpy中創建波形來繞過這一點。

下面是正確的FFT應該看起來像什麼(即基本完善,除了峯的加寬由於窗口)

alt text http://i43.tinypic.com/1rvsqx.png

這裏有一個從一個已經修剪(和非常相似的FFT,從諧波到的精確模式的波形1000Hz附近的高3次諧波)

alt text http://i44.tinypic.com/mt4avd.png 這是我用來產生這些

from numpy import * 
from pylab import ion, plot, draw, show, xlabel, ylabel, figure 

sample_rate = 20000. 
times = arange(0, 10., 1./sample_rate) 
wfm0 = sin(2*pi*200.*times) 
wfm1 = sin(2*pi*500.*times) *(10.-times)/10. 
wfm = wfm0+wfm1 
# int test 
#wfm *= 2**8 
#wfm = wfm.astype(int16) 
#wfm = wfm.astype(float) 
# abs test 
#wfm = abs(wfm) 
# clip test 
#wfm = clip(wfm, -1.2, 1.2) 

fft_length = 5*2048. 
total_num_samps = len(times) 
num_fft = (total_num_samps/fft_length) - 2 
temp = zeros((num_fft,fft_length), float) 

for i in range(num_fft): 
    temp[i,:] = wfm[i*fft_length:(i+1)*fft_length] 
pts = fft_length/2+1 
data = (abs(fft.rfft(temp, fft_length))/(pts))[:pts] 

x_axis = arange(pts)*sample_rate*.5/pts 
spec_range = pts 
plot(x_axis, data[2], linewidth=3) 
xlabel("freq (Hz)") 
ylabel('abs(FFT)') 
show() 
+0

您提出了一個很好的觀點,即FFT是否對信號進行削波,混疊或增加噪聲的最快速和最簡單的測試僅僅是對信號進行FFT,然後對其進行反變換並查看是否恢復原始信號。我的直覺是,FFT工作正常,而這僅僅是誤解FFT輸出的一種情況。 – las3rjock 2009-06-01 14:37:53

2

FFT的,因爲它們是窗口和sampled原因aliasing和抽樣頻域爲好。在時域進行濾波只是在頻域中進行乘法運算,因此您可能只需要應用一個濾波器,該濾波器僅將每個頻率乘以所用濾波器的函數值。例如,在通帶中乘以1,其餘的都是零。意想不到的值可能是由於頻率較高的頻率摺疊到您看到的頻率而造成的。 original signal needs to be band limited to half your sampling rate或者您將得到aliasing。更令人擔憂的是混淆是扭曲了感興趣的區域,因爲對於這個頻段你想知道頻率是來自預期頻率。

要記住的另一件事是,當你從波形文件中獲取一段數據時,你用數字乘以方波。這會導致sinx/x與頻率響應卷積以最小化這種情況,您可以將原始窗口信號乘以Hanning window之類的值。

1

值得一提的一維FFT的代碼,該第一個元素(索引[0])包含DC(零頻)項,元素[1:N/2]包含正頻率,元素[N/2+1:N-1]包含負頻率。由於您未提供代碼示例或有關FFT輸出的附加信息,因此我不能排除「不存在頻率下的噪聲功率值」不僅僅是頻譜的負頻率的可能性。


EDITHere在純Python用一個簡單的測試例程認定的矩形脈衝,[1.,1.,1.,1.,0.,0.,0.,0.]的FFT實現的基數爲2的FFT的一個例子。您可以在codepad運行的例子,看到該序列的FFT是

[0j,     Negative frequencies 
(1+0.414213562373j), ^
0j,      | 
(1+2.41421356237j),  | 
(4+0j),    <= DC term 
(1-2.41421356237j),  | 
0j,      v 
(1-0.414213562373j)] Positive frequencies 

注意,代碼打印出來的傅里葉係數在上升頻率的順序,即從最高負頻率高達DC,然後達到最高正頻率。

+0

好點。此外,在numpy中,函數fftfreq以從fft函數返回的順序返回頻率倉。 – tom10 2009-06-01 13:44:13

+0

謝謝,我的計算波似乎很好,但我仍然遇到處理合成波文件的問題,要麼是我的文件處理例程有問題,要麼是波編輯器本身。即使沒有,這是我無法解決的泄漏,我仍然需要找到一種安全的方式排除非重要的值。當然,它的對數尺度看起來更糟糕。 – 2009-06-02 01:45:29

1

我不知道你的問題是否能夠真正回答任何具體問題。

但這裏有幾件事情,從我自己的經歷寫的FFT嘗試:

  • 確保您按照奈奎斯特規則
  • 如果您正在查看的FFT的線性輸出...你將無法看到自己的信號,並認爲一切都被打破。確保您正在查看FFT幅度的dB。 (即「plot(10 * log10(abs(fft(x))))」)
  • 通過提供生成的數據(如純音)爲FFT()函數創建unitTest。然後將相同的生成數據送入Matlab的FFT()。做兩個輸出數據系列之間的絕對值差異,並確保最大絕對值差異是像10^-6(即唯一的區別是由小浮點錯誤引起的)
  • 確保您windowing your data

如果所有這三件事情都起作用,那麼你的fft很好。而你的輸入數據可能是問題。

時間doamin剪輯顯示爲在特定規則的間隔更小振幅的頻域的信號的鏡像。