Levinson-Durbin算法：

S为信号向量，P为预测阶数

% Levinson-DurbinËã·¨

function a = LD(s, P)

N = length(s);          %Ò»Ö¡ÓïÒô³¤¶È

% Çó×ÔÏà¹Ø¾ØÕór(j),j = 0:P;

for j = 0 : P

r(j+1) = sum(s(1:N-j).*s(j+1:N));

end

a = zeros(1, P);

E = r(1);

for i = 1 : P

temp = 0;

for j = 1 : i-1

temp = temp + a(j) * r(i-j+1);

end

ki = (r(i+1) - temp) / E;

a(i) = ki;

a_temp = a;          %Êý¾Ý½»²æÊ¹ÓÃ£¬·ÀÖ¹Êý¾Ý¸üÐÂ´íÎó

for j = i-1 : -1: 1

a(j) = a_temp(j) - ki * a_temp(i-j);

end

E = (1 - ki * ki) * E;

End

Matlab主程序：

clear all, close all,

% ¶ÁÐÅºÅ

P = 10;

% »­ÓïÒô²¨ÐÎ

figure(1);

subplot(2,2,1);plot(speech);title('Ô­Ê¼ÓïÒô');

subplot(2,2,2);plot(abs(fft(speech)));title('Ô­Ê¼ÓïÒôÆµÆ×');

% 1. ²ÎÊýÉèÖÃ

winSize = floor(0.025 * fs);  % Ö¡³¤25ms,Ò»Ö¡8000*0.025¸öµã¹²200¸öµã

nFrame = floor(length(speech) / winSize); % ¼ÆËãÖ¡Êý

speech = filter([1,-0.94], 1, speech);

subplot(2,2,3);plot(speech);title('Ô¤¼ÓÖØÓïÒô');

subplot(2,2,4);plot(abs(fft(speech)));title('Ô¤¼ÓÖØÓïÒôÆµÆ×');

% Ô¤ÏÈ·ÖÅäÄÚ´æ

a = zeros(1, nFrame*P);

err = zeros(1,length(speech));

frameData = zeros(1,P+winSize);

% ººÃ÷´°

hamwin = hamming(2*P+winSize)';

for m = 1:nFrame

% È¡Ò»Ö¡Êý¾Ý

startPos = (m-1) * winSize + 1;

endPos = startPos + winSize - 1;

ai_start = (m-1) * P +1;

ai_end = ai_start + P - 1;

if m == 1

frameData(P+1:P+winSize) = speech(startPos:endPos)';

else

frameData = speech(startPos - P:endPos)';

end

% ¼ÓººÃ÷´°

frameData = frameData .* hamwin(1:P+winSize);

a(ai_start:ai_end) = LD(frameData(P+1:P+winSize), P);

%Îó²îÊä³ö

for i = 1 : winSize

err(startPos+i) = frameData(P+i) -
sum(a(ai_start:ai_end).*frameData(P+i-1:-1:P+i-P));

end

buildspeech(startPos:endPos) =
filter(1,[1,-1*[a((m-1)*P+1:m*P)]],err(startPos:endPos));

buildspeech(startPos:endPos) = buildspeech(startPos:endPos) ./
hamwin(P+1:P+winSize);     %È¥ººÃ÷´°

end

% È¥¼ÓÖØ

buildspeech = filter(1, [1,-0.94], buildspeech);

%Ð´Êý¾ÝÎÄ¼þ

wavwrite(err,fs,nBits,'error.wav');

wavwrite(buildspeech,fs,nBits,'buildspeech.wav');

figure(2);

subplot(2,1,1);plot(err);

subplot(2,1,2);plot(abs(fft(err)));

figure(3);

subplot(2,1,1);plot(buildspeech);

subplot(2,1,2);plot(abs(fft(buildspeech)));

wavplay(buildspeech,fs,'async');

GitHub

Gitee