【语音分析】基于matlab短时自相关基音周期检测+LPC预测增益计算【含Matlab源码 1517期】

原创

Matlab科研 2021-11-16 13:38:47 ©著作权

文章标签 matlab 开发语言语音识别参考文献深度学习 文章分类 代码人生

©著作权归作者所有：来自51CTO博客作者Matlab科研的原创作品，请联系作者获取转载授权，否则将追究法律责任

一、案例简介

1 实验要求
1.录一段声音，内容为“间隔读出数字0-4”，转为单声道、8khz、wav格式；
2.每帧帧长20ms160个样点，使用汉明窗hamming加窗；
3.设定能量阈值，判断语音中无声、有声区间；
4.计算短时自相关 R（k），k取20-100，计算基音周期和基音频率；
5.信息预测：数字“0”的预测增益E0/Ep（原始信号能量/残差信号能量），并求预测系数a1-a10，阶数p取10。

2 思路
先对声音进行窗函数划分帧，160个样点一帧，计算各帧短时能量，设定能量阈值筛选出5段浊音段，每段分别进行以下操作：段中的每一帧计算自相关函数（k为自相关滞后时间，取20-100），通过波峰最大值点计算得到该帧浊音的基音频率，进而计算该段的基音频率，进而计算5段浊音段的基音频率。
对于信息预测，对于数字“0”的语音段的某一帧（如第101帧），利用lpc 函数得到阶数取10下的预测系数，进而计算预测信号及残差信号，利用预测增益公式E0/Ep（原始信号能量/残差信号能量）计算得到预测增益。

二、部分源代码

clc;
clear all;
close all;        

% 读取语料
[y,fs,wmode,fidx]=readwav('number.wav','p',-1,-1);
figure(1);stem(y,'.');title('原始音频波形');%显示声音的波形
%l=length(y);

%分帧
%f=enframe(y,hamming(160));                   %分帧长:20ms(160样点)
%figure(2);stem(f,'.');title('分帧后的音频波形');%显示声音的波形


%加窗，计算短时能量
for k = 1:nFrames

    E(k) = sum(y_sub.^2); 
end
figure(2);stem(E,'.');title('全区间的短时能量图');

%如何设定能量阈值为0.001得到安静、有声的区间？
for i = 1:nFrames
    if E(i)<0.001
        E(i)=0
    end
end
figure(3);stem(E,'.');title('有声区间的短时能量图');%处理后的帧的短时能量图像

%计算短时自相关
n=160;
for m=1:length(y)/n            %对每一帧求短时自相关函数，每帧的Rm最大值存在N(m)里

        Rm(k)=0;
     
            Rm(k)=Rm(k)+y(i+(m-1)*n)*y(i+k+(m-1)*n);
        end
    end

end                            

%计算基音周期和基音频率
%N=N+10;                        %补回前边去掉的10个点
%T=N/8;                         %算出对应的周期

figure(4);stem(T,'.');axis([0 length(T) 0 20]);
xlabel('帧数(n)');ylabel('周期(ms)');title('初始未处理的基音周期');

T1= medfilt1(T,7);             %去除野点，中值平滑
for k = 1:nFrames
    if E(k)==0
        T1(k)=0;
    end
end
figure(5);stem(T1,'.');axis([0 length(T1) 0 20]);
xlabel('帧数(n)');ylabel('周期(ms)');title('中值滤波后的基音周期');


for k = 1:nFrames
    if E(k)==0
        F1(k)=0;
    end
end
figure(6);stem(F1,'.');
xlabel('帧数(n)');ylabel('频率(Hz)');title('基音频率');

%计算数字“0”的差分方程系数和预测增益

y1=y(99*160:100*160-1);  %取数字“0”的一帧（第100帧），进行汉明加窗
A=lpc(y1.*w,p);  %得到系数（11个系数中第一个系数值为1）
est_Frame=filter([0 -A(2:end)],1,y1);%estimate frame(lp)预测第100/101帧?
FFT_est=fft(est_Frame);
%y2=y(100*160:101*160-1);

figure(7);
%subplot(221);plot(1:3360,y1,1:3360,est_Frame,'-r');grid;title('原始语音帧 vs.预测后的语音帧');
subplot(221);plot(y1);grid;title('原始语音段');
%subplot(222);plot(est_Frame);grid;title('预测语音段');
subplot(222);plot(est_Frame);grid;title({'预测语音段及10个预测系数',[num2str(A(2)),',',num2str(A(3)),',',num2str(A(4)),',',num2str(A(5)),',',num2str(A(6)),',',num2str(A(7)),',',num2str(A(8)),',',num2str(A(9)),',',num2str(A(10)),',',num2str(A(11))]});
subplot(223);plot(err);grid;title('误差');
subplot(224);plot(prodictive_gain);grid;title('预测增益');