%程序6.1：lpc_coefficients.m
%此程序的功能是用自相关法求使信号s均方预测误差为最小的预测系数 
%算法为Levinson-Durbin快速递推算法 
%首先对输入语音进行分帧，并给出LPC分析阶次
  fid=fopen('sx86.txt','r');
  p1=fscanf(fid,'%f')
  fclose(fid);
  p=filter([1 -0.68], 1, p1)             %预加重滤波
  x=fra(320,160,p);                   %将预加重后语音分帧，每帧320个样点，帧重叠160
  x=x(60,:);                         %取第60帧输入信号进行处理,x为行向量
  s=x';                              %x为行向量，s为列向量
  N=16;                            %LPC阶次N=16
  p=N;                             %获得LPC阶次
 n=length(s);                         %获得信号长度                     

  for i=1:p                                                
    Rp(i,1)=sum(s(i+1:n).*s(1:n-i))    %求向量的自相关函数 .*两个同维矩阵相应元素相乘
    %Rn(i)=sum(s(1:N-i).*s(1+i:N)); 
  end
  Rp=Rp(:)                         %将自相关函数变为列向量
  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,1)=Rp(1,1)/Rp_0; 
  a(1,1)=k(1,1); 
  Ep(1,1)=(1-k(1,1)^2)*Ep_0; 
  %i>=2以后使用递归算法 
  if p>1                                                   
    for i=2:p 
        k(i,1)=(Rp(i,1)-sum(a(1:i-1,i-1).*Rp(i-1:-1:1)))/Ep(i-1,1);
  % 求（6-21）式k(i)
        a(i,i)=k(i,1);  %求（6-24）式a(i)
        Ep(i,1)=(1-k(i,1)^2)*Ep(i-1,1); %求（6-23）式Ei
        for j=1:i-1 
            a(j,i)=a(j,i-1)-k(i,1)*a(i-j,i-1);%求（6-25）式a(j,i)
        end 
    end
  end
  c=-a(:,p); %将a矩阵从第1到最后一行的第p列元素乘以(-1)赋给c，c 即为最后求得的LPC系数，但不包括第一个%系数1
  %得到最终的LPC系数a1，此处a1为行向量
  a1(1,1)=1.0;%赋上第一个LPC系数1
  for i=2:p+1
    a1(1,i)= c(i-1,1);%得到第2到第p+1个LPC系数
  end     
