2012-02-06 108 views
3

我试图通过倒谱法找到频率。对于我的测试,我得到了以下文件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的

enter image description here

估计为共享倒谱系频率估计来源:

%% 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)) 

回答

4

256可能太小,如果做任何有用的采样率是44.1 kHz。在这种情况下,FFT的分辨率将是44100/256 = 172 Hz。如果你想要的分辨率为10赫兹的顺序,那么你可以使用FFT大小为4096.

+0

好吧,我改变了我的块到4096。 但即便如此,我的结果是错误的:-( 也许我会需要使用二次插值找到MAX – ederwander 2012-02-06 15:24:06

+1

它可能会帮助,如果你能为这两个最初的数幅度FFT和最终倒谱加图你的帖子,我想我现在可能会忽略零填充,并尽可能简单地保持它。 – 2012-02-06 16:12:47

+0

感谢保罗你看到更新:-) – ederwander 2012-02-06 18:53:17

2

倒谱方法对谐波含量较高的信号效果最好,而对靠近纯正弦波的信号。

最好的测试信号可能更像是重复的,非常接近均等间隔的时域脉冲(每个FFT窗口越多越好),这应该会产生一些接近于频域中重复的等间隔峰值的信号,这应该显示为倒谱的激励部分。脉冲响应将在倒谱的较低共振峰部分中表示。

+0

很酷,这只是一个测试,我会用这个复杂的波形:-) – ederwander 2012-02-06 19:13:52

2

我也有类似的问题,所以我从

我得到一个执行相同的帧的连续评估,然后选择中间值重用你的代码和改进结果的质量的一部分一致的结果。

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] 
+0

你可以使用峰值(最大)位置的二次插值,并可以进一步提高结果... – ederwander 2013-07-30 14:53:33

+0

什么是海明? – David 2014-04-29 20:07:05

+0

@David在这里看看https://en.wikipedia.org/wiki/Window_function – ederwander 2016-05-06 10:44:24