[x,fs]=wavread('11.wav');%读入数据
y=resample(x,4,11);%重采样将频率变为8k
y=double(y); % Convert to double
N=512;
y1=y(1:N);
y1=y1;
w1=hanning(N);
s=y1.*w1;%加窗 取一帧数据
subplot(3,1,1),plot(s);
title('原始语音波形')
p=20;%预测阶数
n=length(s); %获得信号长度
for i=1:p %测试向量
Rp(i)=sum(s(i+1:n).*s(1:n-i)); %求向量的自相关函数
%Rn(i)=sum(s(1:N-i).*s(1+i:N));
end
Rp_0=s'*s; %即Rn(0)
Ep=zeros(p,1); %Ep为p阶最佳线性预测反滤波能量
k=zeros(p,1); %k为偏相关系数
a=zeros(p,p); %以上为初始化
%i=1的情况需要特殊处理,也是对p=1进行处理
Ep_0=Rp_0;
k(1)=Rp(1)/Rp_0;
a(1,1)=k(1);
Ep(1)=(1-k(1)^2)*Ep_0;
%i=2起使用递归算法
if p>1
for i=2:p
k(i)=(Rp(i)-sum( a(1:i-1,i-1).*Rp(i-1:-1:1)'))/Ep(i-1);
a(i,i)=k(i);
Ep(i)=(1-k(i)^2)*Ep(i-1);
for j=1:i-1
a(j,i)=a(j,i-1)-k(i)*a(i-j,i-1);
end
end
end
ai=a(:,p);
LP = filter([1 -ai(2:end)'],1,s); % 建立语音帧的正则方程
FFTlp = fft(LP);
E = s - LP; % 预测误差
subplot(3,1,2),plot(1:N,s,1:N,LP,'-r');grid;
title('原始语音和预测语音波形')
subplot(3,1,3),plot(E);grid;
title('预测误差');