ケプストラル法で周波数を見つけようとしています。私のテストでは、次のファイルhttp://www.mediacollege.com/audio/tone/files/440Hz_44100Hz_16bit_05sec.wav、周波数 440Hz のオーディオ信号を取得しました。
次の式を適用しました。
ケプストラム = IFFT (log FFT (s))
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 を取得するにはどうすればよいでしょうか?
**
アップデート:
**
ソースを matlab で書き直したところ、すべてが機能するようになりました。440 Hz と 250 Hz の周波数でテストを行いました ...
440Hzの場合、441Hzは悪くありません
250Hzの場合、結果に近い249.1525Hzが得られます
ピークをケプストラム値にする簡単な方法を 1 つ作成しました。
Quadract 補間を使用して最大値を見つけると、より良い結果が得られると思います。
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))