线性预测分析matlab程序
注意:程序只支持单声道wav文件,可以自行录制一句话。
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,
% ¶ÁÐźÅ
[speech, fs, nBits] = wavread('speech.wav');
%[speech, fs, nBits] = wavread('a4.wav');
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');