function c=lpc_coefficients(s,p)
%此函數(shù)的功能是用自相關(guān)法求信號(hào)s使均方預(yù)測(cè)誤差為最小的預(yù)測(cè)系數(shù)
%算法為Durbin快速遞推算法
%注：s必須是一維向量 p為聲道模型的階數(shù)

n=length(s);                                             %獲得信號(hào)長度                    
for i=1:p                                                %測(cè)試向量
    Rp(i)=sum(s(i+1:n).*s(1:n-i));                       %求向量的自相關(guān)函數(shù)
    %Rn(i)=sum(s(1:N-i).*s(1+i:N));
end
Rp_0=s'*s;                                               %即Rn(0)

Ep=zeros(p,1);                                           %Ep為p階最佳線性預(yù)測(cè)反濾波能量
k=zeros(p,1);                                            %k為偏相關(guān)系數(shù)
a=zeros(p,p);                                            %以上為初始化

%i=1的情況需要特殊處理,也是對(duì)p=1進(jìn)行處理
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

c=a(:,p);                                                %c(j)=a(j,p)
%c2 =lpc(s,p)                                            %信號(hào)處理工具箱中的函數(shù)