python - 倒谱法的基频

Question

我正在尝试通过倒谱法查找频率。对于我的测试，我得到了以下文件 http://www.mediacollege.com/audio/tone/files/440Hz_44100Hz_16bit_05sec.wav , 频率为 440Hz 的音频信号。

我应用了以下公式:

倒谱 = IFFT(对数 FFT (s))

我得到了 256 个 block ，但我的结果总是错误的......

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

如果您的采样率为 44.1 kHz，256 可能太小而无法做任何有用的事情。在这种情况下，您的 FFT 分辨率将为 44100/256 = 172 Hz。如果您想要大约 10 Hz 的分辨率，那么您可以使用 4096 的 FFT 大小。