声音模型的识别与建立.docx
《声音模型的识别与建立.docx》由会员分享，可在线阅读，更多相关《声音模型的识别与建立.docx（23页珍藏版）》请在冰豆网上搜索。
声音模型的识别与建立
2016理工大学第二届研究生数学建模竞赛参赛
承诺书
我们仔细阅读了理工大学研究生数学建模竞赛的竞赛规则.
我们完全明白，在竞赛开始后参赛队员不能以任何方式（包括、电子、网上咨询等）与队外的任何人（包括指导教师）研究、讨论与赛题有关的问题。
我们知道，抄袭别人的成果是违反竞赛规则的,如果引用别人的成果或其他公开的资料（包括网上查到的资料），必须按照规定的参考文献的表述方式在正文引用处和参考文献中明确列出。
我们重承诺，严格遵守竞赛规则，以保证竞赛的公正、公平性。
如有违反竞赛规则的行为，我们将受到严肃处理。
我们授权理工大学研究生数学建模竞赛组委会，可将我们的论文以任何形式进行公开展示（包括进行网上公示，在书籍、期刊和其他媒体进行正式或非正式发表等）。
我们参赛选择的题号是（从组委会提供的试题中选择一项填写）：
我们的参赛报名号为（如果组委会设置报名号的话）：
所属学校（请填写完整的全名）：
参赛队员（打印并签名）：
1.
2.
3.
指导教师或指导教师组负责人（打印并签名）：
日期：
年月日
评阅编号（由组委会评阅前进行编号）：
2016理工大学研究生数学建模竞赛
编号专用页
评阅编号（由组委会评阅前进行编号）：
评阅记录（可供评阅时使用）：
评
阅
人
评
分
备
注
声音识别模型的建立与评价
摘要
本文通过使用MATLAB软件对声音的时域和频域特征进行了提取，研究特征向量提取方法及SVM核函数和参数选取对识别结果的影响，分析特征提取算法的优缺点，以及不同核函数以及惩罚参数对识别性能的影响。
通过使用支持向量机法建模，基本达到了区分正常声音与非正常声音的目的。
最后提出用低通滤波滤除白噪声。
关键词：
特征向量，支持向量机，核函数，低通滤波，白噪声
从物理上讲，声音是由物体振动产生的一种波，并通过空气作用于人的耳鼓，使人们能够感知。
声音的具有四种性质：
1）音高：
振动发出的声波有不同的频率，称为“音高”；2）强弱：
声音的强弱是由振幅决定的，振幅是代表物体振动强度的特定单位，一般用分贝（dB）来表示。
3）长短：
一般把声音的发展过程分为四个阶段，分别是触发、衰减、保持和消失。
这四个阶段称为“包络”，包络的发生时间，也就是一个声音的长短。
4）音质：
音质好的声音听起来悦耳，相反则让人不适。
问题一
利用matlab中的sound函数，播放出声音信号，试听并比较正常和非正常开门声音的差别，利用plot函数绘制出具体的声音波形图，总结差别在哪些方面？
试听：
我们使用sound函数播放声音样本。
在仔细听了正常开门及非正常开门的声音后，发现了他们之间的差别：
正常开门声音很短促，即声音的长短度短，且其强弱度相对较低。
相反非正常开门声音持续时间长，强弱度高。
画图：
我们选取了三组正常开门及三组非正常开门的声音，使用plot函数画出声音波形图，如下图1所示：
由图1可以看出正常开门的声音波形比较疏松，所以音调就低。
非正常开门的声音波形比较密集，所以声调就高。
我们还可以看出前者声波比较集中，而后者则比较散，即跨度大。
这一点很好理解，正常用钥匙开门所需时间肯定比盗贼撬锁所需时间短，所以就造成了这种现象。
程序源代码见附录一。
图1正常与非正常开门声音波形图
问题二
利用合适的时域或（和）频域特征表达个声音信号，建立特征向量，写出提取特征向量的具体方法和程序代码。
首先，我们对两种样本求其均值。
由于正常开门的第一组数据有人声干扰，故舍去。
然后用plot函数画出声音波形图，如下图2所示：
根据所分析的参数类型，语音信号分析可以分成时域分析和变换域（频域、倒谱域）分析。
其中时域分析方法是最简单、最直观的方法，因为它直接对语音信号的时域波形进行分析。
接下来我们先进行时域分析。
图2正常开门声与非正常开门声均值
短时平均能量：
定义n时刻某语音信号的短时平均能量En为：
（1）
式中，N为窗长，可见短时能量为一帧样点值的平方和。
一般我们认为声音在10-30ms之是稳定的，取帧长也在10-30ms之，而帧移通常取5-15ms之间，所以取N=55、95、125、165。
如下图3、4为正常开门及非正常开门N取不同值时短时能量函数随帧数的变化曲线，其中横坐标为帧数。
由图3、图4可以看出，N=55,N=95时的曲线不够平滑，而N=165的曲线又过于平滑，故选取N=125时的曲线。
图3正常开门声音的短时能量曲线
图4非正常开门声音的短时能量曲线
通过观察短时能量曲线，可以看出正常开门时的能量比较集中且数值小，非正常开门时的能量比较分散且数值大。
容易想到，这些现象与前面问题一的结论是相吻合的。
短时平均过零率：
短时平均过零率是指每帧信号通过零值的次数。
对有时间横轴的连续声音信号，可以观察到声音的时域波形通过横轴的情况。
在离散时间声音信号情况下，如果相邻的采样具有不同的代数符号就称为发生了过零，因此可以计算过零的次数。
（2）
上式为短时平均过零率的公式，其中，sgn[]为符号函数，即
（3）
短时平均过零率曲线如下图5、图6所示：
图5正常开门声音的短时平均过零率
图6非正常开门声音的短时平均过零率
由图5、图6可以看出前者比后者的过零率要高。
短时自相关函数：
自相关函数用于衡量信号自身时间波形的相似性。
对于声音来说，采用短时分析方法，可以定义短时自相关函数为：
（4）
短时自相关函数曲线如下图7、图8所示：
图7正常开门声音的自相关函数曲线
图8非正常开门声音的自相关函数
频域特征：
频域分析主要是对声音波形进行傅里叶变换，如下图9所示：
图9声音波形的傅里叶变换
由图9可以看出正常开门的对数幅度要比非正常开门的对数幅度小的多。
程序源代码见附录二。
问题三
建立声音识别模型（二分类模型），利用模型区分正常和非正常声音，评价模型的好坏。
问题三利用支持向量机构建二分类器对正常与非正常敲门进行分类预测。
支持向量机是Cortes和Vapnik于1995年提出的，与传统分类器比较，它在解决小样本、非线性及高维模式识别中表现出许多特有的优势，并能够推广应用到函数拟合等其他机器学习问题中。
支持向量机方法是建立在统计学习理论的VC维理论和结构风险最小原理基础上的，根据有限的样本信息在模型的复杂性和学习能力之间寻求最佳折衷，以期获得最好的推广能力。
图10支持向量机分类器结构图
采集样本共有80组音频信号，从中选择60组数据做为训练数据训练网络，20组数据作为测试数据测试网络分类能力。
下面是准备数据步骤：
Step1:
根据给定的数据，选定训练集和测试集；
Step2:
为训练集与测试集选定标签集；
Step3:
利用训练集进行训练分类器得到model；
Step4:
根据model，对测试集进行测试集得到accuracyrate；
训练集：
trainset（）;分别取正常开门与非正常开门数据的一半作为训练集；
测试集：
testset（）;分别取正常开门与非正常开门数据的另一半作为测试集；
标签集：
label（）;取bedroom的数据为正类标签为1；forse的数据为负类标签为-1.
下面为svm二分类器的分类结果：
图11基于支持向量机的分类器
上图将展示了训练与测试的结果，结果等于1表示正常声音，结果为-1为非正常开门声，通过学习训练得到上图结果。
程序源代码见附件三。
通过仿真，验证了支持向量机在解决分类问题中的优势，今后可以将这方法用于更深层次的语音识别中。
问题四
试利用特征选择或变换，对特征向量进行优化，并利用参数优化技术优化模型的参数，使识别模型的准确率提高。
问题四其中支持向量机的核函数主要有：
线性核函数，多项式核函数 ，径向基核函数，两层感知器核函数。
不同的核函数直接影响SVM的分类能力，使用不同的核函数SVM分类的效果不同，所需要的时间也不同，下图是使用短时能量值和均值为特征参数，选取RBF为核函数的分类效果图：
图12SVMRBF分类图
图13分类模型错误率
此模型相比于问题三中建立的SVM分类模型分类正确率明显提高，达到了85.65%，可见合适的核函数选取的重要意义。
程序源代码见附录四。
问题五
白噪声，是指功率谱密度在整个频域均匀分布的噪声，即其功率谱密度为
（5）
式中，
为常数，
。
由于1和
为一对傅里叶变换对，所以白噪声的自相关函数为
（6）
由式
（1）、式
（2）可知，白噪声的自相关函数仅在
时才不为零,而对于其他任意的
，白噪声的自相关函数都为零，即在任意两个不同时刻上的随机变量都是不相关的。
对于一个线性时不变系统（如滤波器），该系统是由它的冲激响应
或等效地由它的频率响应
表征，这里的
和
是一对傅里叶变换对。
如果令
为系统的输人信号，
是输出信号。
系统的输出可以表示成如下形式：
（7）
如果
是平稳随机过程
的样本函数，那么
就是随机过程
的样本函数。
可以求出输出的自相关函数是
（8）
由于我们知道自相关函数和功率谱密度函数是一对傅里叶变换对，所以可以得到输出过程的功率密度谱，即为相关函数的傅里叶变换：
（9）
由此可以看出，输出信号的功率谱密度就是输入信号的功率谱密度乘以系统的频率响应的模的平方。
当输入随机过程是白噪声时，输出随机过程的自相关特性和功率密度谱将完全由系统的频率响应所决定。
我们设计了一个低通滤波器，对加了白噪声的声音信号进行滤波，即可去除白噪声。
低通滤波（Low-passfilter）是一种过滤方式，规则为低频信号能正常通过，而超过设定临界值的高频信号则被阻隔、减弱。
但是阻隔、减弱的幅度则会依据不同的频率以及不同的滤波程序（目的）而改变。
在数字图像处理领域，从频域看，低通滤波可以对图像进行平滑去噪处理。
滤波器将这部分不需要的频率滤除，就实现了去噪的效果。
滤波前与滤波后的声音波形对比如下图所示：
图14滤波前后信号波形对比
实验表明低通滤波具有良好的去噪能力，该方法对添加的噪声处理后，与原始信号效果几乎拟合。
证明低通滤波对信号去噪是可行的。
程序源代码见附录五
附录
附录一：
clc;clear
fori=5:
7
j=2*（i-4）-1;
a=['正'num2str（i）];
load（a）
subplot（3,2,j）;
plot（y）
title（['正'num2str（i）]）
xlabel（'Time'）;
ylabel（'Amplitude'）;
a=['非'num2str（i）];
load（a）
subplot（3,2,j+1）
plot（y）
title（['非'num2str（i）]）
xlabel（'Time'）;
ylabel（'Amplitude'）;
end
saveas（gcf,'q1.jpg'）;
附录二：
均值
functionfei=feipj（）
fei=zeros（55125,1）;
fori=1:
40
a=['非'num2str（i）];
load（a）
fei=fei+y;
end
fei=fei/40;
短时平均能量
functionenergy=nengliang（y）
s=fra（55,55,y）;
s2=s.^2;%一帧各样点的能量
energy=sum（s2,2）;%求一帧能量
subplot（2,2,1）%定义画图数量和布局
plot（energy）%画N=55时的语音能量图
xlabel（'帧数'）%横坐标
ylabel（'短时能量E'）%纵坐标
title（'N=55'）%曲线标识
s=fra（95,95,y）;
s2=s.^2;
energy=sum（s2,2）;
subplot（2,2,2）
plot（energy）%画N=95时的语音能量图
xlabel（'帧数'）
ylabel（'短时能量E'）
title（'N=95'）
s=fra（125,125,y）;
s2=s.^2;
energy=sum（s2,2）;
subplot（2,2,3）
plot（energy）%画N=125时的语音能量图
xlabel（'帧数'）
ylabel（'短时能量E'）
title（'N=125'）
s=fra（165,165,y）;
s2=s.^2;
energy=sum（s2,2）;
subplot（2,2,4）
plot（energy）%画N=165时的语音能量图
xlabel（'帧数'）
ylabel（'短时能量E'）
title（'N=165'）
functionf=fra（len,inc,x）
fh=fix（（（size（x,1）-len）/inc）+1）;
f=zeros（fh,len）;
i=1;n=1;
whilei<=fh
j=1;
whilej<=len
f（i,j）=x（n）;
j=j+1;n=n+1;
end
n=n-len+inc;
i=i+1;
end
短时平均过零率
functionf=zcro（x）
f=zeros（size（x,1）,1）;%生成全零矩阵
fori=1:
size（x,1）
z=x（i,:
）;%提取一行数据
forj=1:
（length（z）-1）;
ifz（j）*z（j+1）<0;
f（i）=f（i）+1;
end
end
end
functionzcr=guoling（x）
s=fra（125,125,x）;%分帧，帧移110
zcr=zcro（s）;%求过零率
figure
（1）;
subplot（2,1,1）
plot（x）;
title（'非'）;
xlabel（'Time'）;
ylabel（'Amplitude'）;
subplot（2,1,2）
plot（zcr）;
xlabel（'帧数'）;
ylabel（'过零次数'）;
title（'原始信号的过零率'）;
短时自相关函数
functionA1=zixiangguan（s1）
N=125;
A=[];%加N=125的矩形窗
fork=1:
125;
sum=0;
form=1:
N-k+1;
sum=sum+s1（m）*s1（m+k-1）;%计算自相关
end
A（k）=sum;
end
fork=1:
125
A1（k）=A（k）/A
（1）;%归一化A（k）;
end
plot（A1）;
title（'非正常开门的自相关函数'）
xlabel（'延时k'）
ylabel（'R（k）'）
傅里叶变换
clc;clear
Fs=11025;
zheng=zhengpj（）;
fei=feipj（）;
T=1/Fs;
L=55125;
t=（0:
L-1）*T;
NFFT=2^nextpow2（L）;
subplot（2,1,1）
Y=fft（zheng,NFFT）/L;
f=Fs/2*linspace（0,1,NFFT/2+1）;
y1=2*abs（Y（1:
NFFT/2+1））;
plot（f,2*abs（Y（1:
NFFT/2+1）））
title（'正常开门的傅里叶变换'）
subplot（2,1,2）
Y=fft（fei,NFFT）/L;
f=Fs/2*linspace（0,1,NFFT/2+1）;
y2=2*abs（Y（1:
NFFT/2+1））;
plot（f,2*abs（Y（1:
NFFT/2+1）））
title（'非正常开门的傅里叶变换'）
saveas（gcf,'fuliye.jpg'）
附录三：
clc;clear
zhn=zeros（40,55125）;
zhn（1,:
）=zhengpj（）';
fori=2:
40
name=['正',num2str（i）];
load（name）
zhn（i,:
）=y';
end
fei=zeros（40,55125）;
fori=1:
40
name=['非',num2str（i）];
load（name）
fei（i,:
）=y';
end
label=[ones（1,641）,-ones（1,641）]';
trainset=[zhn（1:
20,:
）;fei（1:
20,:
）];
testset=[zhn（21:
40,:
）;fei（21:
40,:
）];
[mtrain,ntrain]=size（trainset）;
[mtest,ntest]=size（testset）;
test_dataset=[trainset;testset];
[dataset_scale,ps]=mapminmax（test_dataset',0,1）;
dataset_scale=dataset_scale';
trainset=dataset_scale（1:
mtrain,:
）;
testset=dataset_scale（（mtrain+1）:
（mtrain+mtest）,:
）;
train=zeros（2,55125）;
test=zeros（2,55125）;
fori=1:
20
train（1,:
）=trainset（i,:
）+train（1,:
）;
train（2,:
）=trainset（i+20,:
）+train（2,:
）;
end
fori=1:
20
test（1,:
）=testset（i,:
）+test（1,:
）;
test（2,:
）=testset（i+20,:
）+test（2,:
）;
end
train1=zeros（2,1282）;
test1=zeros（2,1282）;
train1（1,:
）=train（1,1:
43:
55125）;
train1（2,:
）=train（2,1:
43:
55125）;
test1（1,:
）=test（1,1:
43:
55125）;
test1（2,:
）=test（2,1:
43:
55125）;
train1=train1';
test1=test1';
p=cvpartition（label,'Holdout',0.2）;
model=svmtrain（train1（p.training,:
）,label（p.training）,'showplot',true）;
C=svmclassify（model,train1（p.test,:
）,'showplot',true）;
err_rate=sum（label（p.test）~=C）/p.TestSize
saveas（gcf,'wentisan2.jpg'）;
附录四：
clear;clc;closeall;
disp（'-------------采用径向基积函数的支持向量机的应用-------------'）;
disp（'输入样本矩阵（每行表示一个测试数据）:
'）;
%x=[randn（1,20）-3randn（1,20）+3;randn（1,40）]'
loadpj1.mat
loadpj2.mat
disp（'样本所属类别:
'）;
a1=pj1（1:
100:
55125）;
a2=pj2（1:
100:
55125）;
a3=[a1a2];
y=[-ones（1,276）ones（1,276）]';
[nsv,alpha,bias,T]=svm168（a3,y,'rbf',10,4）;
fprintf（'运行时间T=%g\n',T）;
fprintf（'支持向量个数：
%g\n',nsv）;
svcplot（a3,y,'linear',alpha,bias）;
disp（'测试数据:
'）;
%X_test=[randn（1,20）-3randn（1,20）+3;randn（1,40）]'
X_test=a3;
disp（'测试数据类别:
'）;
Y_test=[-ones（1,276）ones（1,276）]';
disp（'测试数据类别判决:
'）;
Output=kernel（'rbf',X_test,a3（find（alpha）,:
）,4）*alpha（find（alpha））+bias
disp（'错误率:
'）;
number_error=sum（abs（Y_test-sign（Output）））/2;
error=number_error/20
saveas（gcf,'wentisan.jpg'）;
附录五：
%低通滤波器的设计
fp=1500;fc=1700;As=100;Ap=1;fs=11025;n=55125;
wc=2*pi*fc/fs;wp=2*pi*fp/fs;
wdel=wc-wp;
beta=0.112*（As-8.7）;
N=ceil（（As-8）/2.285/wdel）;
wn=kaiser（N+1,beta）;
ws=（wp+wc）/2/pi;
b=fir1（N,ws,wn）;
figure（3）;
freqz（b,1）;
%对语音信号进行滤波
loadpj1.mat
y1=pj1+（-0.15+0.3*rand（size（pj1）））*max（pj1）;
x=fftfilt（b,y1）;
X=fft（x,n）;
figure（4）;
subplot（2,1,1）;plot（y1）;
title（'滤波前信号的波形'）
subplot（2,1,2）;plot（x）;
title（'滤波后信号的波形'）
saveas（gcf,'quchubaizaosheng.jpg'）;
升级会员 