通過倒譜法的基本頻率

Question

我試圖通過倒譜法找到頻率。對於我的測試，我得到了以下文件http://www.mediacollege.com/audio/tone/files/440Hz_44100Hz_16bit_05sec.wav，一個頻率爲440Hz的音頻信號。

我已經應用了以下公式：

倒= IFFT（FFT記錄（一個或多個））

我得到256塊，但我的成績永遠是錯的...

from numpy.fft import fft, ifft 
import math 
import wave 
import numpy as np 
from scipy.signal import hamming 

index1=15000; 
frameSize=256; 
spf = wave.open('440.wav','r'); 
fs = spf.getframerate(); 
signal = spf.readframes(-1); 
signal = np.fromstring(signal, 'Int16'); 
index2=index1+frameSize-1; 
frames=signal[index1:int(index2)+1] 

zeroPaddedFrameSize=16*frameSize; 

frames2=frames*hamming(len(frames)); 
frameSize=len(frames); 

if (zeroPaddedFrameSize>frameSize): 
    zrs= np.zeros(zeroPaddedFrameSize-frameSize); 
    frames2=np.concatenate((frames2, zrs), axis=0) 

fftResult=np.log(abs(fft(frames2))); 
ceps=ifft(fftResult); 

posmax = ceps.argmax(); 

result = fs/zeroPaddedFrameSize*(posmax-1) 

print result 

對於這種情況如何得到結果= 440？

**

UPDATE：

**

嗯，我改寫了我的源代碼在MATLAB，現在一切似乎工作，我的440的頻率做了測試Hz和250 Hz ...

對於440Hz我得到441Hz不壞

對於250Hz我得到249.1525Hz接近結果

我做了一個簡單的方法來獲得峯值倒譜值。

我想我可以找到更好的結果使用四角插值找到最大值！

我繪製我的結果440Hz的

估計爲共享倒譜系頻率估計來源：

%% ederwander Cepstral Frequency (Matlab) 
waveFile='440.wav'; 
[y, fs, nbits]=wavread(waveFile); 

subplot(4,2,1); plot(y); legend('Original signal'); 

startIndex=15000; 
frameSize=4096; 
endIndex=startIndex+frameSize-1; 
frame = y(startIndex:endIndex); 

subplot(4,2,2); plot(frame); legend('4096 CHUNK signal'); 

%make hamming window 
win = hamming(length(frame)); 


%samples multplied by hamming window 
windowedSignal = frame.*win; 


fftResult=log(abs(fft(windowedSignal))); 
subplot(4,2,3); plot(fftResult); legend('FFT signal'); 

ceps=ifft(fftResult); 

subplot(4,2,4); plot(ceps); legend('ceps signal'); 

nceps=length(ceps) 

%find the peaks in ceps 

peaks = zeros(nceps,1); 

k=3; 

while(k <= nceps - 1) 
    y1 = ceps(k - 1); 
    y2 = ceps(k); 
    y3 = ceps(k + 1); 
    if (y2 > y1 && y2 >= y3) 
     peaks(k)=ceps(k); 
    end 
k=k+1; 
end 

subplot(4,2,5); plot(peaks); legend('PEAKS'); 

%get the maximum ... 
[maxivalue, maxi]=max(peaks) 



result = fs/(maxi+1) 


subplot(4,2,6); plot(result); %legend('Frequency is' result); 

legend(sprintf('Final Result Frequency =====>>> (%8.3f)',result)) 

Answer 1

256可能太小，如果做任何有用的採樣率是44.1 kHz。在這種情況下，FFT的分辨率將是44100/256 = 172 Hz。如果你想要的分辨率爲10赫茲的順序，那麼你可以使用FFT大小爲4096.

Answer 2

倒譜方法對諧波含量較高的信號效果最好，而對靠近純正弦波的信號。

最好的測試信號可能更像是重複的，非常接近均等間隔的時域脈衝（每個FFT窗口越多越好），這應該會產生一些接近於頻域中重複的等間隔峯值的信號，這應該顯示爲倒譜的激勵部分。脈衝響應將在倒譜的較低共振峯部分中表示。

Answer 3

我也有類似的問題，所以我從

我得到一個執行相同的幀的連續評估，然後選擇中間值重用你的代碼和改進結果的質量的一部分一致的結果。

def fondamentals(frames0, samplerate): 
    mid = 16 
    sample = mid*2+1 
    res = [] 
    for first in xrange(sample): 
     last = first-sample 
     frames = frames0[first:last] 
     res.append(_fondamentals(frames, samplerate)) 
    res = sorted(res) 
    return res[mid] # We use the medium value 

def _fondamentals(frames, samplerate):  
    frames2=frames*hamming(len(frames)); 
    frameSize=len(frames); 
    ceps=ifft(np.log(np.abs(fft(frames2)))) 
    nceps=ceps.shape[-1]*2/3 
    peaks = [] 
    k=3 
    while(k < nceps - 1): 
     y1 = (ceps[k - 1]) 
     y2 = (ceps[k]) 
     y3 = (ceps[k + 1]) 
     if (y2 > y1 and y2 >= y3): peaks.append([float(samplerate)/(k+2),abs(y2), k, nceps]) 
     k=k+1 
    maxi=max(peaks, key=lambda x: x[1]) 
    return maxi[0]