clc;
close?all;
clear?all;
[xfs]=wavread（‘song_65KHz.wav‘）;%讀入數據
y=resample（x411）;%重采樣將頻率變為8k
N=160;
y1=y（1:N）;
w1=hanning（N）;
y2=y1.*w1;%加窗?取一幀數據
p=10;%預測階數

%首先求自相關函數
r=zeros（1p+1）;
for?k=1:p+1
????sum=0;
??for?m=1:N+1-k
??????sum=sum+y2（m）.*y2（m-1+k）‘;
??end
??r（k）=sum;
end
%根據durbin算法求線性預測系數
k=zeros（1p）;
k（1）=r（2）/r（1）;
a=zeros（pp）;
a（11）=k（1）;
e=zeros（1p）;
e（1）=（1-k（1）^2）*r（1）;
%遞推過程
for?i=2:p
????c=zeros（1i）;
????sum=0;
????for?j=1:i-1
????????sum=sum+（a（i-1j）.*r（i+1-j））;
????end
????c（i）=sum;
????k（i）=（r（i+1）-c（i））/e（i-1）;%求反射系