Фундаментальная частота по методу Кепстрала

Я пытаюсь найти частоты по методу Кепстраля. Для моих тестов я получил следующий файл: http://www.mediacollege.com/audio/tone/files/440Hz_44100Hz_16bit_05sec.wav , аудиосигнал с частотой 440 Гц.

Я применил следующую формулу:

cepstrum = 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 Гц и 250 Гц …

Для 440 Гц я получаю 441 Гц не плохо

Для 250 Гц я получаю 249,1525 Гц вблизи результата

Я сделал один простой способ получить пики в кепстральные значения.

Я думаю, что я могу найти лучшие результаты, используя квадратичную интерполяцию, чтобы найти максимум!

Я рисую свои результаты для оценки 440 Гц

введите описание изображения здесь

Обмен источниками для Cepstral Оценка частоты:

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

256, вероятно, слишком мал, чтобы делать что-нибудь полезное, если частота дискретизации составляет 44,1 кГц. Разрешение вашего БПФ в этом случае будет 44100/256 = 172 Гц. Если вам требуется разрешение порядка 10 Гц, вы можете использовать размер FFT 4096.

Цепстральные методы лучше всего работают с сигналами с высоким содержанием гармоник, а не с сигналами, близкими к чистым синусоидам.

Лучший тестовый сигнал может быть чем-то более похожим на повторяющиеся, очень близкие друг к другу импульсы во временной области (чем больше на каждом окне FFT, тем лучше), что должно приводить к чему-то близкому к повторяющимся равноотстоящим пикам в частотной области, которые должны отображаться как возбуждающая часть кепстра. Импульсная характеристика будет представлена ​​в нижней формантной части кепстра.

У меня была аналогичная проблема, поэтому я повторно использовал часть кода и улучшал качество результата, выполняя непрерывную оценку одного и того же кадра, а затем выбирая значение среды из

Я получаю последовательный результат.

 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]