% function F0 = f0_by_amdf(x)
% 根据 W-AMDF 求基音周期
% 输入 x 为语音数据 程序已好
clear
[x,fs] = wavread('F:\语音\语料\字\白m11.wav');
% x = x/max(abs(x));
N = 512;
win = rectwin(N);
shift = N/2;
x = enframe(x,win,shift);
n = size(x,1);
gudian = zeros(1,n);
am = zeros(n,N);
for k = 1:N % W-AMDF
am(:,k) = sum(abs(x(:,1:N+1-k)-x(:,k:N)),2)/(N-k+1);
end
N = length(yF);
rw_L = sum(abs(yF(1+L:N)-yF(1:N-L)))/(N-L-1);
th = zeros(1,n);
gudian = zeros(1,n);
for i = 1:n
th(i) = min(am(i,44:round(0.6*N)));% 寻找一个可靠的最小值,0.6是个经验值
end
for i = 2:n
for j = 32:320 % 去除50~500Hz(round(fs/500):round(fs/50))之外的野点
if th(i) == am(i,j)
gudian(i) = j;
if round((gudian(i)+1)/(gudian(i-1)+1))==2 | ...
round((gudian(i-1)+1)/(gudian(i)+1))==2
gudian(i)=abs(gudian(i)-gudian(i-1)); % 修正2倍周期点
end
end
end
end
for i = 1:n
if gudian(i) ~= 0
gudian(i) = fs/gudian(i); % 转换成相应的频率
end
end
subplot(211);
stem(gudian); title('各帧原始基音频率');
xlabel('帧数'); ylabel('F / Hz');
axis tight;
% 第二次去野
ave = mean(gudian);
index = find(abs(gudian-ave) > ave/2); % 阈值需根据实际情况调整
gudian(index) = 0;
% 后处理. 两种方法可单独使用,也可结合使用
% 3 点中值平滑
jj = 1;
while jj < n-1 % 或jj+2 <= n
average = (gudian(jj)+gudian(jj+1)+gudian(jj+2))/3;
for kk = 1:3
if abs((gudian(jj-1+kk))-average) > average/2 % 平滑精度可调
gudian(jj-1+kk) = 0; % 将野点置零
end
end
jj = jj+3; % 3点平滑
end
% 线性平滑.平滑越多,平滑段之间的阶跃性破坏越严重
% for i = 1:3 % 平滑次数人为设定
% gudian = smooth(gudian);
% end
subplot(212);
stem(gudian); title('平滑后的各帧基音频率');
xlabel('帧数'); ylabel('F / Hz');
axis tight;