clear;clc;
fid=fopen('female.dat');
x1=fread(fid,16000,'short');
figure;plot(x1);title('读信号');
grid
采用fread函数以二进制数据读入数据。设定信号长度为16000,精度为16位整数型。语音信号图如下:
%分帧
wlen=8000*20*10.^(-3); inc=160; % 给出帧长和帧移
win=boxcar(wlen); % 给出矩形窗
N=length(x1); % 信号长度
time=(0:N-1)/8000; % 计算出信号的时间刻度
X=enframe(x1,win,inc); % 分帧
X=X';
fn=size(X,2); % 求出帧数
%去均值
for p=1:89
ave=mean(X(:,p));
X(:,p)=X(:,p)-ave;
for o=1:160
XX((p-1)*160+o)=X(o,p);
end
end
figure;plot(1:14240,XX,'k'); axis([5000 5400 ,-inf,inf]);xlabel('时间(s)');ylabel('时间波形');title('展开信号');
hold on;
b1=[0.008233 -0.004879 0.007632 -0.004879 0.008233];
a1=[1 -3.6868 5.8926 -5.0085 2.2518 -0.4271];
[H,w]=freqz(b1,a1,512);
for r=1:89
sf_1(:,r)=filter(b1,a1,X(:,r));
end
plot(1:14240,sf,'b');xlabel('时间(s)');ylabel('幅值');title('椭圆滤波');axis([5000 5400 ,-inf,inf]) %改变横轴范围
hold on;
b2=[1 1 1 1 1 1 1 1 1];
a2=[9 0 0 0 0 0 0 0 0];
[H,w]=freqz(b2,a2,512);
for r=1:89
sf_2(:,r)=filter(b2,a2,sf_1(:,r));
end
plot(1:14240,sf,'r');xlabel('时间(s)');ylabel('幅值');title('数值滤波');axis([5000 5400 ,-inf,inf]) %改变横轴范围
legend('去均值','低通滤波','数值滤波');
由于
是一个极小值,因此在此程序里任取为0.0000000001。在此给出各帧短时平均能量ELP计算代码和结果图:
for i=1 : fn
for z=2:wlen
S_LPN(1,1)=sf(1,(1+(i-1)*inc))*sf(1,(1+(i-1)*inc));
S_LPN(1,z)=S_LPN(1,(z-1))+sf(1,(z+(i-1)*inc))*sf(1,(z+(i-1)*inc));
end
En(i)=10*log10(S_LPN(1,160)/160+0.0000000001);
end
为了避免选择实际基音的倍数作为估计值,可采用如下后处理措施。
在三个基音区域内计算ρ(τ),分别得到各自区域的最大值及对应的延迟,记为ρmax1,ρmax2,ρmax3和T1,T2,T3。三个区域分别为:
区域1:80 ~ 147
区域2:40 ~ 79
区域3:20 ~ 39
按如下方法确定最优基音延迟Topt
Topt = T1
ρmax = ρmax1
if ρmax2 ≥ cρmax
ρmax = ρmax2
Topt = T2
if ρmax3 ≥ cρmax
ρmax = ρmax3
Topt = T3
end
end
c为经验因子,经过测试,发现男性取值为0.92,女性取值为0.94时比较合适。
for m=2: fn
%取一帧
if En(1,m)>32.62
for tao=20:39 %对于不同的延迟tao
rxx3(1,1)=sf(1+(m-1)*inc)*sf(1-tao+(m-1)*inc);
s_n(1,1)=sf(1+(m-1)*inc)*sf(1+(m-1)*inc);
s_n_tao3(1,1)=sf(1-tao+(m-1)*inc)*sf(1-tao+(m-1)*inc);
for j=2:wlen %分别计算分子和分母
rxx3(1,j)=rxx3(j-1)+sf(j+(m-1)*inc)*sf(j-tao+(m-1)*inc); %分子
s_n(1,j)=s_n(j-1)+sf(j+(m-1)*inc)*sf(j+(m-1)*inc);
s_n_tao3(1,j)=s_n_tao3(j-1)+sf(j-tao+(m-1)*inc)*sf(j-tao+(m-1)*inc);
end
Rxx3(1,tao-19)=rxx3(1,160)/((s_n(1,160)*s_n_tao3(1,160))^0.5);
end
[Rxx_tao_max3,index3]=max(Rxx3); %最大互相关函数峰值=Rxx_tao-max、延迟tao=index
for tao=40:79 %对于不同的延迟tao
rxx2(1,1)=sf(1+(m-1)*inc)*sf(1-tao+(m-1)*inc);
s_n_tao2(1,1)=sf(1-tao+(m-1)*inc)*sf(1-tao+(m-1)*inc);
for j=2:wlen %分别计算分子和分母
rxx2(j)=rxx2(j-1)+sf(j+(m-1)*inc)*sf(j-tao+(m-1)*inc); %分子
s_n_tao2(j)=s_n_tao2(j-1)+sf(j-tao+(m-1)*inc)*sf(j-tao+(m-1)*inc);
end
Rxx2(tao-39)=rxx2(160)/((s_n(160)*s_n_tao2(160))^0.5);
end
[Rxx_tao_max2,index2]=max(Rxx2); %最大互相关函数峰值=Rxx_tao-max、延迟tao=index
for tao=80:147 %对于不同的延迟tao
rxx1(1,1)=sf(1+(m-1)*inc)*sf(1-tao+(m-1)*inc);
s_n_tao1(1,1)=sf(1-tao+(m-1)*inc)*sf(1-tao+(m-1)*inc);
for j=2:wlen %分别计算分子和分母
rxx1(j)=rxx1(j-1)+sf(j+(m-1)*inc)*sf(j-tao+(m-1)*inc); %分子
s_n_tao1(j)=s_n_tao1(j-1)+sf(j-tao+(m-1)*inc)*sf(j-tao+(m-1)*inc);
end
Rxx1(tao-79)=rxx1(160)/((s_n(160)*s_n_tao1(160))^0.5);
end
[Rxx_tao_max1,index1]=max(Rxx1); %最大互相关函数峰值=Rxx_tao-max、延迟tao=index
%最优基因延迟Topt
Rxx_tao_max=Rxx_tao_max1;Topt=index1+79;
c=0.92;
if (Rxx_tao_max2>=c*Rxx_tao_max)
Rxx_tao_max=Rxx_tao_max2;Topt=index2+39;
end
if (Rxx_tao_max3>=c*Rxx_tao_max)
Rxx_tao_max=Rxx_tao_max3;Topt=index3+19;
end
P_max(1,m)=Rxx_tao_max;
%清浊判决
P_avr=(Rxx_tao_max1+Rxx_tao_max2+Rxx_tao_max3)/3; %前三个峰的平均值
Z_period=Rxx_tao_max+P_avr; %周期性水平量
if (Z_period>0.79) %Topt在允许的基频范围内
Topt_tao(1,m)=Topt;
else
Topt_tao(1,m)=0;
end
else
Topt_tao(1,m)=0;
end
end
figure;plot(1:97,Topt_tao) % 画出基音周期
title('基音周期');ylabel('基音周期值'); xlabel('帧数');
axis([0 100,0 120]) %改变横轴范围
hold on;