分别提取演唱者音频与模版音频的特征参数，通过特征参数的差值得到演唱得分。

---

源码分三个部分：
1. main.m
主文件，运行入口
音频的读入，对基音特征进行分段DTW，根据DTW的结果计算得分
计算方式为：(cos((pl/tpl)*pi)+1)*50
其中，pl为待评分音频的结果 tpl为对照音频的结果
2. solve.m
基音特征提取函数
流程：
预加重->分帧->加窗->求短时平均能量->求自相关(求基音)
主要参数：
ek ：高频提升参数
L ：帧长（ms）
RL ：帧移（ms）
kee：清浊音阈值
调整参数可以调整系统的整体特性以适应不同需求
3. DTW.m
动态时间规整函数
N^2复杂度，分段加快处理速度

---

%main.m
%西安电子科技大学 windroid
clear all;
%默认音源已提前进行统一化处理
%读取数据
file1 = 'C:\Users\97628\Desktop\音乐\music\暗香0.wma';
file2 = 'C:\Users\97628\Desktop\音乐\music\暗香1.wma';
[e1,Fs1] = audioread(file1);
[e2,Fs2] = audioread(file2);f1 = solve(e1, Fs1);
f2 = solve(e2, Fs2);figure(1);
subplot(211);plot(f1);title('file1基音频率');xlabel('帧数');ylabel('幅度');
subplot(212);plot(f2);title('file2基音频率');xlabel('帧数');ylabel('幅度');%DTW是为了小范围拟合 所以输入的音频需要提前对齐
%DTW -> n^2算法
%按帧分段
D = 100;
pl = 0;
for i=1:D:min(length(f1), length(f2))-Dpl = pl + DTW(f1(:,i:i+D), f2(:,i:i+D));
end
%空对照组
tpl = 0;
zf = zeros(1, length(f2));
for i=1:D:min(length(f1), length(zf))-Dtpl = tpl + DTW(f1(:,i:i+D), zf(:,i:i+D));
endif(tpl<pl)info = '音源差距太大，无法评分'
else%满分100%正弦曲线score = (cos((pl/tpl)*pi)+1)*50
end

---

%solve.m
function f = solve(e, Fs)
%e:音频数据; Fs:采样率
%返回基音频率 列向量
%
%函数过程：
%   预加重->分帧->加窗->求短时平均能量->求自相关(求基音)
%关键参数：
%   ek L RL kee
%西安电子科技大学 windroid%预加重
ek = -0.98; %-1 ~ 0
e = filter([1,ek],1,e);  %un为经过高频提升后的时域信号
log = '预加重完成'
un = e;
%figure(1);plot(e);title('原始语音信号');xlabel('样点数');ylabel('幅度');
%figure(2);plot(un);title('预加重后的语音信号');xlabel('样点数');ylabel('幅度');L = 30; % 帧长(ms)
RL = 20; % 帧移(ms)
t = length(e)*1000/Fs; % 时长(ms)
N = floor((t-L)/RL); % 帧数
ns = floor(Fs*L/1000); %帧采样点数
rnl = floor(Fs*RL/1000); %帧移点数un(length(un)+1:ceil(t/L)*ns, :) = 0; % 信号尾0扩充
su = size(un);
if su(2) == 2un = (un(:, 1) + un(:, 2))/2;
end%分帧
e1 = zeros(N, ns);
for i = 1:1:Ne1(i, :) = un((i-1)*rnl+1:(i-1)*rnl+ns, 1)';
end
%figure(3);plot(e1(100, :));title('分帧后的第100帧信号');xlabel('样点数');ylabel('幅度');
log = '分帧完成'%加窗
e2 = zeros(N, ns);
for i = 1:1:Ne2(i, :) = e1(i, :).*(hamming(ns))';
end
%figure(5);plot(e2(100, :));title('加窗后的第100帧信号');xlabel('样点数');ylabel('幅度');
%figure(6);plot(ee);title('加窗后的短时平均能量');xlabel('帧数');ylabel('能量');
log = '加窗完成'%{
%削波
emax = max(abs(e2'))';%中心削波
t = 0.7;
e3 = zeros(N, ns);
for i = 1:1:N for j = 1:1:nstc = t*emax(i);if e2(i, j) > tce3(i, j) = e2(i, j) - tc;elseif e2(i, j) < -tce3(i, j) = e2(i, j) + tc;elsee3(i, j) = 0;endend
end
%figure(7);plot(e3(100, :));title('中心削波后的第100帧信号');xlabel('样点数');ylabel('幅度');%TODO %三电平削波
%}
%{
%声频分离
LE = 1024;
t = zeros(LE, 1);
t(51:950, 1) = hanning(900);
ef = zeros(N, LE);
for i= 1:1:Nfe = fft(e1(i, :), LE);fe = fe.*(t');ef(i,:) = ifft(fe, LE);
end
%}%短时平均能量
ee = sum(e1.^2, 2)';
ne = sum(ee)/N;%figure(8);
%subplot(211);plot(ee);title('短时平均能量');xlabel('帧数');ylabel('能量');log = '正在求基音'
%求自相关
LE = ns;
re = e2;
tR = zeros(1, N);
kee = 1;%阈值
for i = 1:1:Nif ee(:, i) < kee*necontinueendR = zeros(1, LE);for k = 1:1:LEfor t = 1:1:LE-kR(1, t) = R(1, t) + re(i, t)*re(i, t + k);endend[~ ,tR(1, i)] = max(R);
end
%subplot(212);plot(tR);title('基音频率');xlabel('帧数');ylabel('幅度');f = tR;

---

%DTW.m
function dist = DTW(t,r)
%西安电子科技大学 windroid
%参考：http://www.cnblogs.com/luxiaoxun/archive/2013/05/09/3069036.html
%↑链接代码有问题
n = length(t);
m = length(r);
% 帧匹配距离矩阵
d = zeros(n,m);
for i = 1:nfor j = 1:md(i,j) = (t(:,i)-r(:,j))^2;end
end
% 累积距离矩阵
D = ones(n,m) * realmax;
D(1,1) = d(1,1);log = '正在DP'
% 动态规划
for i = 1:nfor j = 1:mif i>1D1 = D(i-1,j);if j>1D2 = D(i-1,j-1);D3 = D(i,j-1);elseD2 = realmax;D3 = realmax;endelseD1 = realmax;D2 = realmax;if j>1D3 = D(i,j-1);else%D3 = realmax;%此时 i == 1 && j == 1D(i, j) = d(i, j);continueendendD(i,j) = d(i,j) + min(min(D1,D2),D3);end
end
dist = D(n,m);

---

测试音频结果：

---

---

---

---

---

---