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])%设置坐标轴范围

绘图展示结果:
matlab求取语音的基音频率、共振峰信息并将其标注在语谱图上

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));       %01024之间计算出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('共振峰标注');

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

3.matlab绘制语音语谱图

figure()
spectrogram(x,wlen,inc,nfft,fs,'yaxis');

绘图展示:
matlab求取语音的基音频率、共振峰信息并将其标注在语谱图上

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('共振峰标注到语谱图');

绘制结果:
matlab求取语音的基音频率、共振峰信息并将其标注在语谱图上文章来源地址https://www.toymoban.com/news/detail-472630.html

5.说明

  1. 实验中pitch函数和spectrogram函数为matlab自带的,如果实验者电脑无法找到这俩函数,可以在网上搜索找到该函数并复制到该文件路径下,作者的matlab版本为matlab 2019a。
  2. 居中表示的.m文件,需要读者将其内容全部复制,并将其命名为对应的matlab函数文件。
  3. 因为时间有限,仅将代码实现部分展现了出来,后续有时间的话,也会继续将内容补上。

到了这里,关于matlab求取语音的基音频率、共振峰信息并将其标注在语谱图上的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处: 如若内容造成侵权/违法违规/事实不符,请点击违法举报进行投诉反馈,一经查实,立即删除!

领支付宝红包 赞助服务器费用

相关文章

  • 代码片记录-使用DDC对两路信号相位差进行求取(matlab实现)

    参考链接: 1、FPGA综合系统设计(七)基于DDC的两路信号相位差检测_fpga相位测量-CSDN博客 (真的太牛了) 2、双通道中频信号数字下变频及相位差估计(FPGA)_根据 i、q 两路信号以如何求出中频信号的频率和相位差-CSDN博客 (通过这个篇博客注意到了反正切的值域问题,里面

    2024年04月14日
    浏览(37)
  • 中文语音标注工具FunASR(语音识别)

    全称  A Fundamental End-to-End Speech Recognition Toolkit (一个语音识别工具) 可能大家用过 whisper (openAi),它【标注英语的确很完美】,【但中文会出现标注错误】或搞了个没说的词替换上去,所以要人工核对,麻烦。 FunASR作用 :能【准确】识别语音,并转成【文字、标出声调】

    2024年02月04日
    浏览(40)
  • Matlab数据处理:用离散数据根据自定义多变量非线性方程拟合解析方程求取参数

    问题:已知xlsx表格[X,Y,Z]的离散取值,希望用  来拟合,用matlab求得[C1,C2,C3,C5,C6]的值 解答: 运行结果:  备注: 1.rsquare=0.8668认为接近1,拟合效果不错 2.fill函数的startpoint如何设置[C1,...C6]得到一个收敛点?(我找了没找到什么设置startpoint好方法,摸索用如下方法找到了一个

    2024年02月11日
    浏览(51)
  • 【matlab进阶学习-7】matlab 图表标注操作

    本文参考:MATLAB04:基础绘图-CSDN博客 plot(x,y,LineSpec) 各参数意义如下: x : 图线上点的x坐标 y : 图线上点的y坐标 LineSpec : 图线的线条设定,三个指定 线型 , 标记符号 和 颜色 的 设定符 组成一个字符串,设定符不区分先后.具体细节请参考 官方文档 . 线型设定符 线型 标记设定符 标

    2024年02月22日
    浏览(33)
  • Matlab FFT变换细节(信号采样频率,FFT变换点数,频率分辨率)

    在做深度学习的故障诊断中,发现代码直接将原始信号fft之后直接将实频域信号输入网络中进行诊断,虽说效果比较不错95% 但因为输入的是双边谱且频率范围远超故障特征频率同时由于单个样本的点数只有1024点,信号的采样频率又特别高12k,导致频率分辨率极低,输入网络的序列

    2023年04月08日
    浏览(97)
  • MATLAB中滤波函数、频率响应函数以及频率响应函数不同表达形式的转换

            频率响应函数的表达式:         对应的z变换的多项表达式:         Z变换的零极点表达形式:         Z变换的二阶因子级联形式: filter函数,仅可以用于零状态响应系统。         y=filter(b,a,x) ;                %b为z变换多项表达式公式中[b0,b1...bM]的矩阵

    2024年02月07日
    浏览(47)
  • 自己编译RustDesk,并将自建ID服务器和key信息写入客户端

            今天总算是把编译环境给折腾清楚了,编译出来了至少能用,但说不上好用,问题还不少,官方的客户端就是要手工填写ID服务器地址和key才可以用,而且还容易被别人白嫖你搭建的服务器,当然如果拿到你编译后的客户端,也是存在被白嫖的可能。这方面还没有找

    2024年04月14日
    浏览(153)
  • 微信小程序实现登录授权,并将获取到的用户授权信息存储到数据库中

    官方开发文档 注意:在实现授权登录时,不要使用测试号进行 wx.getUserProfile使用文档 config文件代码如下 路由模块 主要代码(如果不想多个文件,可以将wxuserHandle.wxuser这个位置内容替换为exports.wxuser的内容) 注意: 在返回token的字符串拼接时,Bearer 后面必须有一个空格 数据

    2024年02月10日
    浏览(59)
  • matlab柱状图的绘制及数值的标注

    目标是绘制以下的柱状图:  代码: 其中:state_x的数据为6*3矩阵  解释下循环中的程序:text(xx(i)-0.5,state_x(i,1),num2str(state_x(i,1)),\\\'HorizontalAlignment\\\',\\\'center\\\',\\\'VerticalAlignment\\\',\\\'bottom\\\'); %  xx(i)-0.5,state_x(i,1)表示为所要标注的位置,一个为横轴上的,一个为纵轴上的; %  num2str(state_x

    2024年02月16日
    浏览(36)
  • 【NLP pytorch】基于标注信息从句子中提取命名实体内容

    给定一个句子和已经通过模型训练标注好的信息,从而提取出句子中的实体内容,如下 输入: (1)句子信息 (2)标注信息

    2024年02月14日
    浏览(63)

觉得文章有用就打赏一下文章作者

支付宝扫一扫打赏

博客赞助

微信扫一扫打赏

请作者喝杯咖啡吧~博客赞助

支付宝扫一扫领取红包,优惠每天领

二维码1

领取红包

二维码2

领红包