matlab求取语音的基音频率、共振峰信息并将其标注在语谱图上
matlab求取语音的基音频率、共振峰信息并将其标注在语谱图上
·
matlab求取语音的基音频率、共振峰信息并将其标注在语谱图上
1.matlab 求取语音的基音频率F0
%% 提取基音频率
[f0,idx] = pitch(x,fs, 'WindowLength',wlen,'OverlapLength',inc);%求取语音的基音频率
fn1=size(f0,1);
result_f0=f0(find(f0<400));%筛选符合条件的基音频率
figure()
plot(result_f0);%绘制基音频率
title('基音频率f0')
xlabel('帧')
ylabel('频率/Hz')
axis([0 fn1 0 400])%设置坐标轴范围
绘图展示结果：
2.matlab求取语音共振峰（F1、F2、F3）
clc;
clear all;
%% 读取语音文件
filename='bluesky1.wav';
[x,fs]=audioread(filename);%读取wav文件
%% 共振峰提取,Frmt即F1,F2,F3
wlen=512;inc=256;
[voiceseg,vsl,SF_a,NF] = VAD(x,fs,wlen,inc);
%提取共振峰特征
Doption=0;
Frmt=Formant_ext2(x,wlen,inc,fs,SF_a,Doption);
F1=(Frmt(1,:));F1_min=min(F1);F1_max=max(F1);F1_mean=mean(F1);
F2=(Frmt(2,:));F2_min=min(F2);F2_max=max(F2);F2_mean=mean(F2);
F3=(Frmt(3,:));F3_min=min(F3);F3_max=max(F3);F3_mean=mean(F3);
用到的以下函数
function [voiceseg,vsl,SF,NF] = VAD(y,fs,wlen,inc)
%UNTITLED2 此处显示有关此函数的摘要
% 此处显示详细说明
N=length(y); % 取信号长度
time=(0:N-1)/fs; % 计算时间刻度
IS=0.25; overlap=wlen-inc; % 设置前导无话段长度
NIS=fix((IS*fs-wlen)/inc +1); % 计算前导无话段帧数
frames=enframe(y,wlen,inc)'; % 分帧
fn=size(frames,2); % 帧数
amp=sum(frames.^2); % 求取短时平均能量
zcr=zc2(frames,fn); % 计算短时平均过零率
ampm = multimidfilter(amp,5); % 中值滤波平滑处理
zcrm = multimidfilter(zcr,5);
ampth=mean(ampm(1:NIS)); % 计算初始无话段区间能量和过零率的平均值
zcrth=mean(zcrm(1:NIS));
amp2=1.1*ampth; amp1=1.3*ampth; % 设置能量和过零率的阈值
zcr2=0.9*zcrth;
frameTime=frame2time(fn, wlen, inc, fs);% 计算各帧对应的时间
[voiceseg,vsl,SF,NF]=vad_param2D_revr(amp,zcr,amp2,amp1,zcr2);% 端点检测
end
%% 其他函数
function [voiceseg,vsl,SF,NF]=vad_param2D_revr(dst1,dst2,T1,T2,T3,T4)
fn=length(dst1); % 取帧数
maxsilence = 8; % 初始化
minlen = 5;
status = 0;
count = 0;
silence = 0;
%开始端点检测
xn=1;
for n=1:fn
switch status
case {0,1} % 0 = 静音, 1 = 可能开始
if dst1(n) > T2 | ... % 确信进入语音段
( nargin==6 & dst2(n) < T4 )
x1(xn) = max(n-count(xn)-1,1);
status = 2;
silence(xn) = 0;
count(xn) = count(xn) + 1;
elseif dst1(n) > T1 | ... % 可能处于语音段
dst2(n) < T3
status = 1;
count(xn) = count(xn) + 1;
else % 静音状态
status = 0;
count(xn) = 0;
x1(xn)=0;
x2(xn)=0;
end
case 2, % 2 = 语音段
if dst1(n) > T1 | ... % 保持在语音段
dst2(n) < T3
count(xn) = count(xn) + 1;
else % 语音将结束
silence(xn) = silence(xn)+1;
if silence(xn) < maxsilence % 静音还不够长，尚未结束
count(xn) = count(xn) + 1;
elseif count(xn) < minlen % 语音长度太短，认为是噪声
status = 0;
silence(xn) = 0;
count(xn) = 0;
else % 语音结束
status = 3;
x2(xn)=x1(xn)+count(xn);
end
end
case 3, % 语音结束，为下一个语音准备
status = 0;
xn=xn+1;
count(xn) = 0;
silence(xn)=0;
x1(xn)=0;
x2(xn)=0;
end
end
el=length(x1);
if x1(el)==0, el=el-1; end % 获得x1的实际长度
if x2(el)==0 % 如果x2最后一个值为0，对它设置为fn
% fprintf('Error: Not find endding point!\n');
x2(el)=fn;
end
SF=zeros(1,fn); % 按x1和x2，对SF和NF赋值
NF=ones(1,fn);
for i=1 : el
SF(x1(i):x2(i))=1;
NF(x1(i):x2(i))=0;
end
speechIndex=find(SF==1); % 计算voiceseg
voiceseg=findSegment(speechIndex);
vsl=length(voiceseg);
end
%% 其他函数
function y=multimidfilter(x,m)
a=x;
for k=1 : m
b=medfilt1(a, 5);
a=b;
end
y=b;
end
%% 其它函数
function zcr=zc2(y,fn)
if size(y,2)~=fn, y=y'; end
wlen=size(y,1); % 设置帧长
zcr=zeros(1,fn); % 初始化
delta=0.01; % 设置一个很小的阈值
for i=1:fn
yn=y(:,i); % 取来一帧
for k=1 : wlen % 中心截幅处理
if yn(k)>=delta
ym(k)=yn(k)-delta;
elseif yn(k)<-delta
ym(k)=yn(k)+delta;
else
ym(k)=0;
end
end
zcr(i)=sum(ym(1:end-1).*ym(2:end)<0); % 取得处理后的一帧数据寻找过零率
end
end
%% 其他函数
function soundSegment=findSegment(express)
if express(1)==0
voicedIndex=find(express); % 寻找express中为1的位置
else
voicedIndex=express;
end
soundSegment = [];
k = 1;
soundSegment(k).begin = voicedIndex(1); % 设置第一组有话段的起始位置
for i=1:length(voicedIndex)-1,
if voicedIndex(i+1)-voicedIndex(i)>1, % 本组有话段结束
soundSegment(k).end = voicedIndex(i); % 设置本组有话段的结束位置
soundSegment(k+1).begin = voicedIndex(i+1);% 设置下一组有话段的起始位置
k = k+1;
end
end
soundSegment(k).end = voicedIndex(end); % 最后一组有话段的结束位置
% 计算每组有话段的长度
for i=1 :k
soundSegment(i).duration=soundSegment(i).end-soundSegment(i).begin+1;
end
end
function Fmap=Formant_ext2(x,wlen,inc,fs,SF,Doption)
y=filter([1 -.99],1,x); % 预加重
xy=enframe(y,wlen,inc)'; % 分帧
fn=size(xy,2); % 求帧数
Nx=length(y); % 数据长度
time=(0:Nx-1)/fs; % 时间刻度
frameTime=frame2time(fn,wlen,inc,fs); % 每帧对应的时间刻度
Msf=repmat(SF',1,3); % 把SF扩展为3×fn的数组
Ftmp_map=zeros(fn,3);
warning off
Fsamps = 256; % 设置频域长度
Tsamps= fn; % 设置时域长度
ct = 0;
numiter = 10; % 循环10次,
iv=2.^(10-10*exp(-linspace(2,10,numiter)/1.9)); % 在0～1024之间计算出10个数
for j=1:numiter
i=iv(j);
iy=fix(length(y)/round(i)); % 计算帧数
[ft1] = seekfmts1(y,iy,fs,10); % 已知帧数提取共振峰
ct = ct+1;
ft2(:,:,ct) = interp1(linspace(1,length(y),iy)',...% 把ft1数据内插为Tsamps长
Fsamps*ft1',linspace(1,length(y),Tsamps)')';
end
ft3 = squeeze(nanmean(permute(ft2,[3 2 1]))); % 重新排列和平均处理
tmap = repmat([1:Tsamps]',1,3);
Fmap=ones(size(tmap))*nan; % 初始化
idx = find(~isnan(sum(ft3,2))); % 寻找非nan的位置
fmap = ft3(idx,:); % 存放非nan的数据
[b,a] = butter(9,0.1); % 设计低通滤波器
fmap1 = round(filtfilt(b,a,fmap)); % 低通滤波
fmap2 = (fs/2)*(fmap1/256); % 恢复到实际频率
Ftmp_map(idx,:)=fmap2; % 输出数据
Fmap=Ftmp_map';
findex=find(Fmap==nan);
Fmap(findex)=0;
% 作图
if Doption
figure(99)
nfft=512;
d=stftms(y,wlen,nfft,inc);
W2=1+nfft/2;
n2=1:W2;
freq=(n2-1)*fs/nfft;
Fmap1=Msf.*Ftmp_map; % 只取有话段的数据
findex=find(Fmap1==0); % 如果有数值为0 ,设为nan
Fmapd=Fmap1;
Fmapd(findex)=nan;
imagesc(frameTime,freq,abs(d(n2,:))); axis xy
m = 64; LightYellow = [0.6 0.6 0.6];
MidRed = [0 0 0]; Black = [0.5 0.7 1];
Colors = [LightYellow; MidRed; Black];
colormap(SpecColorMap(m,Colors)); hold on
plot(frameTime,Fmapd,'w'); % 叠加上共振峰频率曲线
title('在语谱图上标出共振峰频率');
xlabel('时间/s'); ylabel('频率/Hz')
end
end
%% 其它函数
function [fmt] = seekfmts1(sig,Nt,fs,Nlpc)
if nargin<4, Nlpc = round(fs/1000)+2; end;
ls=length(sig); % 数据长度
Nwin = floor(ls/Nt); % 帧长
for m=1:Nt,
lpcsig = sig((Nwin*(m-1)+1):min([(Nwin*m) ls]));% 取来一帧信号
if ~isempty(lpcsig),
a = lpc(lpcsig,Nlpc); % 计算LPC系数
const=fs/(2*pi); % 常数
rts=roots(a); % 求根
k=1; % 初始化
yf = [];
bandw=[];
for i=1:length(a)-1
re=real(rts(i)); % 取根之实部
im=imag(rts(i)); % 取根之虚部
formn=const*atan2(im,re); % 计算共振峰频率
bw=-2*const*log(abs(rts(i))); % 计算带宽
if formn>150 & bw <700 & formn<fs/2 % 满足条件方能成共振峰和带宽
yf(k)=formn;
bandw(k)=bw;
k=k+1;
end
end
[y, ind]=sort(yf); % 排序
bw=bandw(ind);
F = [NaN NaN NaN]; % 初始化
F(1:min(3,length(y))) = y(1:min(3,length(y))); % 输出最多三个
F = F(:); % 按列输出
fmt(:,m)=F/(fs/2); % 归一化频率
end;
end;
end
function frameTime=frame2time(frameNum,framelen,inc,fs)
% 分帧后计算每帧对应的时间
frameTime=(((1:frameNum)-1)*inc+framelen/2)/fs;
end
function f=enframe(x,win,inc)
nx=length(x(:)); % 取数据长度
nwin=length(win); % 取窗长
if (nwin == 1) % 判断窗长是否为1，若为1，即表示没有设窗函数
len = win; % 是，帧长=win
else
len = nwin; % 否，帧长=窗长
end
if (nargin < 3) % 如果只有两个参数，设帧inc=帧长
inc = len;
end
nf = fix((nx-len+inc)/inc); % 计算帧数
f=zeros(nf,len); % 初始化
indf= inc*(0:(nf-1)).'; % 设置每帧在x中的位移量位置
inds = (1:len); % 每帧数据对应1:len
f(:) = x(indf(:,ones(1,len))+inds(ones(nf,1),:)); % 对数据分帧
if (nwin > 1) % 若参数中包括窗函数，把每帧乘以窗函数
w = win(:)'; % 把win转成行数据
f = f .* w(ones(nf,1),:); % 乘窗函数
end
end
绘制共振峰F1,F2,F3
% 绘制共振峰
figure()
plot(F1,'r','LineWidth',0.3)
hold on
plot(F2,'m','LineWidth',0.3)
plot(F3,'k','LineWidth',0.3)
fn=size(F1,2);%帧数
%设置x轴
frame2time=linspace(0,fn,9);
times=0:0.5:4;
xticks(frame2time);
xticklabels(times);
xlim([0 fn]);%设置x轴范围
ylim([0 fs/2]);%设置x轴范围
%设置图例字体及大小
h=legend('F1','F2','F3');
set(h,'FontName','Times New Roman','FontSize',11,'FontWeight','normal')
xlabel('时间/s');
ylabel('频率/Hz');
title('共振峰标注');
共振峰展示展示
3.matlab绘制语音语谱图
figure()
spectrogram(x,wlen,inc,nfft,fs,'yaxis');
绘图展示：
4.将共振峰绘制到语谱图上显示
figure()
wlen=512;
inc=256;
nfft=wlen;
spec_data=spectrogram(x,wlen,inc,nfft,fs,'yaxis');
spec_data=abs(spec_data).*abs(spec_data);%取模
spec_data=10*log10(spec_data);%取对数
dim2freq=linspace(0,fs/2,size(spec_data,1));
xlabels=0:1:size(spec_data,2);
imagesc(xlabels,dim2freq,spec_data);%绘制语谱图
colorbar();%显示colorbar()
axis xy;%颠倒顺序
hold on
plot(F1,'r','LineWidth',0.3)
plot(F2,'m','LineWidth',0.3)
plot(F3,'k','LineWidth',0.3)
%设置x轴
frame2time=linspace(0,size(spec_data,2),9);
times=0:0.5:4;
xticks(frame2time);
xticklabels(times);
%设置图例字体及大小
h=legend('F1','F2','F3');
set(h,'FontName','Times New Roman','FontSize',11,'FontWeight','normal')
xlabel('时间/s');
ylabel('频率/Hz');
title('共振峰标注到语谱图');
绘制结果：
5.说明
- 实验中pitch函数和spectrogram函数为matlab自带的，如果实验者电脑无法找到这俩函数，可以在网上搜索找到该函数并复制到该文件路径下，作者的matlab版本为matlab 2019a。
- 居中表示的.m文件，需要读者将其内容全部复制，并将其命名为对应的matlab函数文件。
- 因为时间有限，仅将代码实现部分展现了出来，后续有时间的话，也会继续将内容补上。
为开发者提供学习成长、分享交流、生态实践、资源工具等服务，帮助开发者快速成长。
更多推荐
6
0- 0
已为社区贡献3条内容
所有评论(0)