matlab求取语音的基音频率、共振峰信息并将其标注在语谱图上_王延凯

网络投稿 02-07 3587

matlab求取语音的基音频率、共振峰信息并将其标注在语谱图上 1.matlab 求取语音的基音频率F02.matlab求取语音共振峰（F1、F2、F3）3.matlab绘制语音语谱图4.将共振峰绘制到语谱图上显示5.说明

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);

用到的以下函数

VAD.m 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 Formant_ext2.m 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 frame2time.m function frameTime=frame2time(frameNum,framelen,inc,fs) % 分帧后计算每帧对应的时间 frameTime=(((1:frameNum)-1)*inc+framelen/2)/fs; end enframe.m 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函数文件。因为时间有限，仅将代码实现部分展现了出来，后续有时间的话，也会继续将内容补上。