基于MATLAB的音频特征提取

 音频文件无法上传。就只上传代码。

 分帧：

 

分帧是指按照一定长度将输入的声音数据分割成若干数据段，通常取N=512个采样点作为一帧，帧与帧之间如果相互连接会使得在后续的计算中丢失数据，所以取256个采样点的一个帧移。

汉宁窗：

 

 

 梅尔滤波器，取个数为64个：

 

 

 

 

 

 

 

 MFSC:

 

 

 

clc;
clear all;
close all;
%读取音频文件的第2-10s
fs = 16000;
[y,fs]=audioread('C:\Users\dengw\Desktop\D32_999.wav',[2*fs,10*fs]);
figure('name','音频文件幅度-采样点波形');
plot(y);
title('音频文件幅度-采样点波形');
ylabel('幅度');
xlabel('采样点');
%横坐标换算为时间；
t0=[0:length(y)-1]/fs;
figure('name','音频文件幅度-时间波形');
plot(t0,y);
title('音频文件幅度-时间波形');
ylabel('幅度');
xlabel('时间/s');
%对音频信号进行分帧
y0=buffer(y,256); %将矩阵等分为256行，最后一列不够的补0；
y1=zeros(512,160);
for j=1:160;
        y1(:,j)=[y0(:,j);y0(:,j+1)];
end
%分帧音频数据幅度-采样点的波形
n=0:511;
figure('name','分帧音频数据幅度-采样点的波形');
for j=1:160;
    plot(n,y1(:,j));
    hold on;
end   
title('分帧音频数据幅度-采样点的波形');
ylabel('幅度');
xlabel('采样点');
%选第10帧、50帧、100帧、150帧数据查看波形；
%第10帧数据波形
figure('name','第10帧数据幅度-采样点波形');
plot(y1(:,10));
title('第10帧数据幅度-采样点波形');
ylabel('幅度');
xlabel('采样点');
%第50帧数据波形
figure('name','第50帧数据幅度-采样点波形');
plot(y1(:,50));
title('第50帧数据幅度-采样点波形');
ylabel('幅度');
xlabel('采样点');
%第100帧数据波形
figure('name','第100帧数据幅度-采样点波形');
plot(y1(:,100));
title('第100帧数据幅度-采样点波形');
ylabel('幅度');
xlabel('采样点');
%第150帧数据波形
figure('name','第150帧数据幅度-采样点波形');
plot(y1(:,150));
title('第150帧数据幅度-采样点波形');
ylabel('幅度');
xlabel('采样点');
%加汉窗
n=0:511;
wn=0.5*(1-cos(2*pi*n/511));
%汉宁窗波形
figure('name','汉宁窗波形');
plot(wn);
title('汉宁窗波形');
ylabel('幅度');
xlabel('采样点');
%加窗
y2=zeros(512,160);
for j=1:160
    for i=1:512
    y2(i,j)=wn(i)*y1(i,j);
    end
end
%加窗后的波形
%音频数据波形
figure('name','分帧音频数据+汉宁窗幅度-采样点的波形');
for j=1:160;
    plot(n,y2(:,j));hold on;
     
end   
title('分帧音频数据+汉宁窗幅度-采样点的波形');
ylabel('幅度');
xlabel('采样点');
%选第10帧、50帧、100帧、150帧数据查看波形；
%第10帧数据波形
figure('name','第10帧数据+汉宁窗幅度-采样点波形');
plot(n,y2(:,10));
title('第10帧数据+汉宁窗幅度-采样点波形');
ylabel('幅度');
xlabel('采样点');
%第50帧数据波形
figure('name','第50帧数据+汉宁窗幅度-采样点波形');
plot(n,y2(:,50));
title('第50帧数据+汉宁窗幅度-采样点波形');
ylabel('幅度');
xlabel('采样点');
%第100帧数据波形
figure('name','第100帧数据+汉宁窗幅度-采样点波形');
plot(n,y2(:,100));
title('第100帧数据+汉宁窗幅度-采样点波形');
ylabel('幅度');
xlabel('采样点');
%第150帧数据波形
figure('name','第150帧数据+汉宁窗幅度-采样点波形');
plot(n,y2(:,150));
title('第150帧数据+汉宁窗幅度-采样点波形');
ylabel('幅度');
xlabel('采样点');
%离散傅里叶变换+求模
f_abs=zeros(512,160);
for j=1:160
    f_abs(:,j)=abs(fft(y2(:,j)));
end
%只取前257个点
f=zeros(257,160);
f=f_abs(1:257,:);
nf=0:256;
%第10帧数据波形
figure('name','第10帧数据+汉宁窗+FFT+求模幅度-采样点波形');
plot(nf,f(:,10));
title('第10帧数据+汉宁窗+FFT+求模幅度-采样点波形');
ylabel('幅度');
xlabel('采样点');    
%第50帧数据波形
figure('name','第50帧数据+汉宁窗+FFT+求模幅度-采样点波形');
plot(nf,f(:,50));
title('第50帧数据+汉宁窗+FFT+求模幅度-采样点波形');
ylabel('幅度');
xlabel('采样点');   
%第100帧数据波形
figure('name','第100帧数据+汉宁窗+FFT+求模幅度-采样点波形');
plot(nf,f(:,100));
title('第100帧数据+汉宁窗+FFT+求模幅度-采样点波形');
ylabel('幅度');
xlabel('采样点');   
%第150帧数据波形
figure('name','第150帧数据+汉宁窗+FFT+求模幅度-采样点波形');
plot(nf,f(:,150));
title('第150帧数据+汉宁窗+FFT+求模幅度-采样点波形');
ylabel('幅度');
xlabel('采样点');   

nff=(0:256)*fs/512;
%第10帧数据波形
figure('name','第10帧数据+汉宁窗+FFT+求模幅度-频率波形');
plot(nff,f(:,10));
title('第10帧数据+汉宁窗+FFT+求模幅度-频率波形');
ylabel('幅度');
xlabel('频率/Hz');   
%第50帧数据波形
figure('name','第50帧数据+汉宁窗+FFT+求模幅度-频率波形');
plot(nff,f(:,50));
title('第50帧数据+汉宁窗+FFT+求模幅度-频率波形');
ylabel('幅度');
xlabel('频率/Hz'); 
%第100帧数据波形
figure('name','第100帧数据+汉宁窗+FFT+求模幅度-频率波形');
plot(nff,f(:,100));
title('第100帧数据+汉宁窗+FFT+求模幅度-频率波形');
ylabel('幅度');
xlabel('频率/Hz'); 
%第150帧数据波形
figure('name','第150帧数据+汉宁窗+FFT+求模幅度-频率波形');
plot(nff,f(:,150));
title('第150帧数据+汉宁窗+FFT+求模幅度-频率波形');
ylabel('幅度');
xlabel('频率/Hz'); 
%梅尔滤波器组
fs=16000;
fl=0; fh=fs/2;
bl=1125*log(1+fl/700);%将频率转换为Mel频率
bh=1125*log(1+fh/700);
p=64;%滤波器个数
nfft=512;%FFT点数
B=bh-bl;
y_mel=linspace(0,B,p+2);%产生0到B之间p+2个数
Fb=700*(exp(y_mel/1125)-1);%将Mel频率转换为频率
W2=nfft/2+1;%fs/2内对应的FFT点数
df=fs/nfft;
freq=(0:W2-1)*df;%采样频率值
bank=zeros(64,W2);%生成一个24行W2列的全零数组
figure('name','梅尔滤波器组');
for k=2:p+1%why从2开始？因为k-1
    f1=Fb(k-1); 
    f2=Fb(k+1); 
    f0=Fb(k);
    n1=floor(f1/df)+1;%f(m-1)在频域中的谱线索引号
    n2=floor(f2/df)+1;%f(m+1)在频域中的谱线索引号
    n0=floor(f0/df)+1;%f(m)在频域中的谱线索引号。f(m)是从0开始，而在MATLAB中数组的索引是从1开始，所以要加1，否则会出现index=0的错误
    for i=1 : W2
        if i>=n1 & i<=n0
            bank(k-1,i)=(i-n1)/(n0-n1);
        elseif i>n0 & i<=n2
            bank(k-1,i)=(n2-i)/(n2-n0);
        end
    end

    plot(freq,bank(k-1,:),'r','linewidth',2); 
    hold on;
    title('梅尔滤波器组');
    ylabel('幅度');
    xlabel('频率'); 
end
grid; 
%第10帧数据的音频特征
fm10=log(bank*f(:,10));
figure('name','第10帧数据的音频特征');
plot(fm10);
title('第10帧数据的音频特征');
ylabel('MFSC');
xlabel('梅尔滤波器组次序');
%第50帧数据的音频特征
fm50=log(bank*f(:,50));
figure('name','第50帧数据的音频特征');
plot(fm50);
title('第50帧数据的音频特征');
ylabel('MFSC');
xlabel('梅尔滤波器组次序');
%第100帧数据的音频特征
fm100=log(bank*f(:,100));
figure('name','第100帧数据的音频特征');
plot(fm100);
title('第100帧数据的音频特征');
ylabel('MFSC');
xlabel('梅尔滤波器组次序');
%第150帧数据的音频特征
fm150=log(bank*f(:,150));
figure('name','第150帧数据的音频特征');
plot(fm150);
title('第150帧数据的音频特征');
ylabel('MFSC');
xlabel('梅尔滤波器组次序');
fm=zeros(64,160);
fm=log(bank*f);

figure('name','音频特征');
x=(1:160)*0.016;%每帧数据512个点，512/16000=0.032s，但帧移=256个点，
y11=1:64;
pcolor(x,y11,fm);
% surf(x,y,z); view(0,90); %等效的写法
shading interp; 
colorbar; colormap(jet);
title('音频特征')
xlabel('时间/s');ylabel('滤波器次序');