lchw2.m

这是一个用matlab实现的语音信号预测编码中的得宾算法部分。还有实验报告。

w=ones(1,240);
n=[0:239];
w(n+1)=0.54-0.46*cos(2*n*pi/239); %汉明窗

p=15;
 temp=WAVREAD('lj',[15001 15000+240]);%一帧240个点  浊音
%temp=WAVREAD('lj',[57501 57500+240]);%一帧240个点  清音
temp=reshape(temp,1,240);
temp=temp.*w; %利用汉明窗对输入语音帧进行加窗
r=zeros(1,p+1);%短时自相关函数
for nn=[0:p],
    for ii=[1:240-nn],
        r(nn+1) =r(nn+1)+ temp(ii)*temp(nn+ii);
    end
end
%%%%%%%%%%%%%%%%%%下面是Durbin算法
k=0;%k表示km
a=ones(p,p);%存储a参数的p*p矩阵
o=zeros(1,p+1);%o表示arpha

m=0;
o(0+1)=r(0+1);%m=0时，o(0)=r(0)

%m=1
m=1;
k=-r(1+1)/o(0+1);
a(1,1)=k;
o(1+1)=o(0+1)*(1-k*k);

for m=2:p
     s=r(m+1);
    for t=1:m-1
        s=s+a(m-1,t)*r(m-t+1);
    end    
    k=-s/o(m-1+1);
    
    for i=1:m-1        
        a(m,i)=a(m-1,i)+k*a(m-1,m-i);
    end
    a(m,m)=k;
    
    o(m+1)=o(m-1+1)*(1-k*k);
end

% 
figure(5)%apha随n的变化曲线
i=[1:p];
plot(i,o(i),'b');
grid



aa=ones(1,p);
for i=1:p        
    aa(i)=a(p,i);
end
h=zeros(1,p+1);
h(1)=1;
for i=1:p        
    h(i+1)=aa(i);
end
e=conv(temp,h);

figure(2)%残差信号
set(1,'Position',[10,35,350,650])
i=[1:240];
plot(i,e(i),'r');
grid
axis([1,240,-0.1,0.1])


figure(3)%加窗后的信号
set(1,'Position',[10,35,350,650])
i=[1:240];
plot(i,temp(i),'b');
grid