Fundamental Frequency by Cepstral Method

Question

I'm trying to find frequencies by the Cepstral method. For my tests I got the following file http://www.mediacollege.com/audio/tone/files/440Hz_44100Hz_16bit_05sec.wav, an audio signal with a frequency of 440Hz.

I've applied the following formula:

cepstrum = IFFT (log FFT (s))

I'm getting 256 chunk, But my results are always wrong ...

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

For this case how get the result = 440?

**

UPDATE:

**

Well i rewrote my source in matlab, and now everything seems to work, i did tests with frequencies of 440 Hz and 250 Hz ...

For 440Hz i get 441Hz not bad

For 250Hz i get 249.1525Hz near result

I did make one simple way to get the peaks into cepstral values.

I think I can find better results using quadract interpolation to find the maximum !

I'm plotting my results for the estimation of 440Hz

Sharing the source for Cepstral Frequency estimation:

%% 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 is probably too small to do anything useful if your sample rate is 44.1 kHz. The resolution of your FFT in this case will be 44100 / 256 = 172 Hz. If you want a resolution of the order of say 10 Hz then you might use an FFT size of 4096.

Answer 2

Cepstral methods work best with signals that have a high harmonic content, not as well on signals that are close to pure sinusoids.

The best test signal might be something more like repetitive, very near equally spaced, impulses in the time domain (the more per FFT window the better), which should produce something close to repetitive equally spaced peaks in the frequency domain, which should show up as the exciter portion of a cepstrum. The impulse response will be represented in the lower formant portion of the cepstrum.

Answer 3

I had a similar problem and so I have reused a part of your code and improvement the quality of the result by performing contiguous evaluation of the same frame and then selecting the medium value from the

I am getting a consistent result.

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]