分别提取演唱者音频与模版音频的特征参数,通过特征参数的差值得到演唱得分。
源码分三个部分:
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))-D
pl = 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))-D
tpl = tpl + DTW(f1(:,i:i+D), zf(:,i:i+D));
end
if(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) == 2
un = (un(:, 1) + un(:, 2))/2;
end
%分帧
e1 = zeros(N, ns);
for i = 1:1:N
e1(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:N
e2(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:ns
tc = t*emax(i);
if e2(i, j) > tc
e3(i, j) = e2(i, j) - tc;
elseif e2(i, j) < -tc
e3(i, j) = e2(i, j) + tc;
else
e3(i, j) = 0;
end
end
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:N
fe = 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:N
if ee(:, i) < kee*ne
continue
end
R = zeros(1, LE);
for k = 1:1:LE
for t = 1:1:LE-k
R(1, t) = R(1, t) + re(i, t)*re(i, t + k);
end
end
[~ ,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:n
for j = 1:m
d(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:n
for j = 1:m
if i>1
D1 = D(i-1,j);
if j>1
D2 = D(i-1,j-1);
D3 = D(i,j-1);
else
D2 = realmax;
D3 = realmax;
end
else
D1 = realmax;
D2 = realmax;
if j>1
D3 = D(i,j-1);
else
%D3 = realmax;
%此时 i == 1 && j == 1
D(i, j) = d(i, j);
continue
end
end
D(i,j) = d(i,j) + min(min(D1,D2),D3);
end
end
dist = D(n,m);
测试音频结果: