1倒谱实验程序

1 倒谱实验程序
hanning=zeros(1,300);
a1=zeros(1,620);N=300;K=250;daopu=zeros(1,300);
%浊音倒谱分析
sp=wavread('D:\');%11000采样频率
qq=0;
for number=1:300
hanning(number)=(1/2)*(1-cos((2*pi)*(number-1)/299));
end
for i=1:10000
qq=qq+sp(i)/10000;
end
for i=1:10000
sp(i)=sp(i)-qq;
end
% for k=1:N
s1=abs(fft(sp(1+100:300+100).*hanning'));
for i=1:300
s2(i)=logm(s1(i));
end
%end
mm=real(ifft(s2));
for i=1:150
daopu(i)=mm(151-i);
end
for i=1:150
daopu(150+i)=mm(301-i);
end
plot(daopu)
hold on
R=zeros(1,250);
%轻音倒谱计算
sp1=wavread('D:\');%11000采样频率
qq=0;
for i=1:10000
qq=qq+sp1(i)/10000;
end
for i=1:10000
sp1(i)=sp1(i)-qq;
%for k=1:N
s1=abs(fft(sp1(1:300).*hanning'));
for i=1:300
s2(i)=logm(s1(i));
end
%end
mm=abs(ifft(s2));
for i=1:150
daopu1(i)=mm(151-i);
end
for i=1:150
daopu1(150+i)=mm(301-i);
end
%end
plot(daopu1,'r')
2．短时幅度差程序
% 短时幅度差求基因周期
a1=zeros(1,620);N=401;K=250;F=zeros(1,250);
sp=wavread('D:\语音信号处理实验3\');%11000采样频率qq=0;
for i=1:10000
qq=qq+sp(i)/10000;
end
for i=1:10000
sp(i)=sp(i)-qq;
end
sp(1:300)=sp(1:300);%.*hanning(1:300)';
for k=1:K
for m=1:N-k+1
F(k)=F(k)+abs(sp(m)-sp(m+k-1))/(N-k+1);
end
end
plot(F)
hold on
F1=zeros(1,250);
sp1=wavread('D:\语音信号处理实验3\');%11000采样频率
qq=0;
for i=1:10000
qq=qq+sp1(i)/10000;
for i=1:10000
sp1(i)=sp1(i)-qq;
end
for k=1:K
for m=1:N-k
F1(k)=F1(k)+abs(sp1(m)-sp1(m+k-1))/(N-k+1);
end
end
plot(F1,'r')
3 短时相关
a1=zeros(1,620);K=450;R=zeros(1,250);
sp=wavread('D:\语音信号处理实验3\');%11000采样频率
for i=1:250
for j=1:601-1-i
R(i)=R(i)+sp(j)*sp(j+i);
end
end
plot(R/14)
hold on
R=zeros(1,250);
sp=wavread('D:\语音信号处理实验3\');%11000采样频率f1=sp(3300:4000);%11000采样频率
for i=1:250
for j=1:601-1-i
R(i)=R(i)+f1(j)*f1(j+i);
end
end
plot(R/14,'r')
4 LPC倒谱分析
N=256;R=zeros(1,N);
p=12;a=zeros(p,p);En=zeros(1,p);
K=zeros(1,p);
sp=wavread('D:\语音信号处理实验3\');%11000采样频率
qq=0;
for i=1:10000
qq=qq+sp(i)/10000;
end
for i=1:10000
sp(i)=sp(i)-qq;
end
for j=1:p
for n=j:N-1
R(j)=R(j)+sp(n)*sp(n-j+1);
end
end
En(1)=R(1);
K(1)=R(2)/R(1);
a(1,1)=K(1);
En(2)=(1-K(1)^2)*En(1);
%K(2)=(R(3)-a(1,1)*R(2))/En(2);
%a(2,2)=K(2);
%a(2,1)=a(1,1)-K(2)*a(1,1);
for i=2:p
po=0;
for j=1:i-1
po=po+a(i-1,j)*R(i-j+1);
end
K(i)=(R(i+1)-po)/En(i);
a(i,i)=K(i);
for j=1:i-1
a(i,j)=a(i-1,j)-K(i)*a(i-1,i-j);
end
En(i+1)=(1-K(i)^2)*En(i);
end
source=sp(1:256);
for i=1:12
mm(i)=a(12,13-i);
end
mm1=zeros(1,300);
for k=1:256
if k<p
for m=1:k
mm1(k)=mm1(k)+mm(p-k+m)*source(m);
end
else
if k<256
for m=1:p
mm1(k)=mm1(k)+mm(m)*source(m+k-p);
end
else
for m=k-p+1:256
mm1(k)=mm1(k)+mm(m-k+p)*source(m);
end
end
end
end
DFTmm=abs(fft(mm1));
Ln=zeros(1,300);
for i=1:60
Ln(i)=logm(DFTmm(i));
end
ifftLn=real(ifft(Ln));
mm=ifftLn;
for i=1:150
daopu1(i)=mm(151-i);
end
for i=1:150
daopu1(150+i)=mm(301-i);
end
plot(daopu1)。