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

enter image description here

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

score 4 · Accepted Answer · answered Feb 06 '12 at 15:14

4

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.

answered Feb 06 '12 at 15:14

Paul R

208,748
37
389
560

Ok i changed my chunk to 4096. but even so my results are wrong :-( Maybe I'll need to use quadratic interpolation to find the MAX – ederwander Feb 06 '12 at 15:24
1

It would probably help if you could add plots for both the initial log magnitude FFT and the final cepstrum to your post. Also I think I would probably omit the zero padding for now and keep it as simple as possible. – Paul R Feb 06 '12 at 16:12

score 2 · Answer 2 · answered Jul 30 '13 at 14:38

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]

you can use quadratic interpolation around peak(Max) position and can further improve the results ... — ederwander, Jul 30 '13 at 14:53

score 2 · Answer 3 · answered Feb 06 '12 at 19:00

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.

Cool, this is just one test i will use this in complex wave forms :-) — ederwander, Feb 06 '12 at 19:13

Fundamental Frequency by Cepstral Method

3 Answers3