聲音模型的識別與建立

1、2016長沙理工大學(xué)第二屆研究生數(shù)學(xué)建模競賽參賽承諾書我們仔細(xì)閱讀了長沙理工大學(xué)研究生數(shù)學(xué)建模競賽的競賽規(guī)則.我們完全明白，在競賽開始后參賽隊員不能以任何方式（包括電話、電子郵件、網(wǎng)上咨詢等）與隊外的任何人（包括指導(dǎo)教師）研究、討論與賽題有關(guān)的問題。我們知道，抄襲別人的成果是違反競賽規(guī)則的, 如果引用別人的成果或其他公開的資料（包括網(wǎng)上查到的資料），必須按照規(guī)定的參考文獻(xiàn)的表述方式在正文引用處和參考文獻(xiàn)中明確列出。我們鄭重承諾，嚴(yán)格遵守競賽規(guī)則，以保證競賽的公正、公平性。如有違反競賽規(guī)則的行為，我們將受到嚴(yán)肅處理。我們授權(quán)長沙理工大學(xué)研究生數(shù)學(xué)建模競賽組委會，可將我們的論文以任何形式進(jìn)行公開展

2、示（包括進(jìn)行網(wǎng)上公示，在書籍、期刊和其他媒體進(jìn)行正式或非正式發(fā)表等）。我們參賽選擇的題號是（從組委會提供的試題中選擇一項填寫）： 我們的參賽報名號為（如果組委會設(shè)置報名號的話）： 所屬學(xué)校（請?zhí)顚懲暾娜?參賽隊員 (打印并簽名) ： 1. 2. 3. 指導(dǎo)教師或指導(dǎo)教師組負(fù)責(zé)人 (打印并簽名)： 日期： 年 月 日評閱編號（由組委會評閱前進(jìn)行編號）：2016長沙理工大學(xué)研究生數(shù)學(xué)建模競賽編 號 專 用 頁評閱編號（由組委會評閱前進(jìn)行編號）：評閱記錄（可供評閱時使用）：評閱人評分備注聲音識別模型的建立與評價摘 要本文通過使用MATLAB軟件對聲音的時域和頻域特征進(jìn)行了提取，研究特征向量提

3、取方法及SVM核函數(shù)和參數(shù)選取對識別結(jié)果的影響，分析特征提取算法的優(yōu)缺點，以及不同核函數(shù)以及懲罰參數(shù)對識別性能的影響。通過使用支持向量機(jī)法建模，基本達(dá)到了區(qū)分正常聲音與非正常聲音的目的。最后提出用低通濾波濾除白噪聲。關(guān)鍵詞：特征向量，支持向量機(jī)，核函數(shù)，低通濾波，白噪聲從物理上講，聲音是由物體振動產(chǎn)生的一種波，并通過空氣作用于人的耳鼓，使人們能夠感知。聲音的具有四種性質(zhì)：1）音高：振動發(fā)出的聲波有不同的頻率，稱為“音高”；2）強弱：聲音的強弱是由振幅決定的，振幅是代表物體振動強度的特定單位，一般用分貝(dB)來表示。3）長短：一般把聲音的發(fā)展過程分為四個階段，分別是觸發(fā)、衰減、保持和消失。這四

4、個階段稱為“包絡(luò)”，包絡(luò)的發(fā)生時間，也就是一個聲音的長短。4)音質(zhì)：音質(zhì)好的聲音聽起來悅耳，相反則讓人不適。問題一利用matlab中的sound函數(shù)，播放出聲音信號，試聽并比較正常和非正常開門聲音的差別，利用plot函數(shù)繪制出具體的聲音波形圖，總結(jié)差別在哪些方面？試聽：我們使用sound函數(shù)播放聲音樣本。在仔細(xì)聽了正常開門及非正常開門的聲音后，發(fā)現(xiàn)了他們之間的差別：正常開門聲音很短促，即聲音的長短度短，且其強弱度相對較低。相反非正常開門聲音持續(xù)時間長，強弱度高。畫圖：我們選取了三組正常開門及三組非正常開門的聲音，使用plot函數(shù)畫出聲音波形圖，如下圖1所示：由圖1可以看出正常開門的聲音波形比較

5、疏松，所以音調(diào)就低。非正常開門的聲音波形比較密集，所以聲調(diào)就高。我們還可以看出前者聲波比較集中，而后者則比較散，即跨度大。這一點很好理解，正常用鑰匙開門所需時間肯定比盜賊撬鎖所需時間短，所以就造成了這種現(xiàn)象。程序源代碼見附錄一。圖1正常與非正常開門聲音波形圖問題二利用合適的時域或（和）頻域特征表達(dá)個聲音信號，建立特征向量，寫出提取特征向量的具體方法和程序代碼。首先，我們對兩種樣本求其均值。由于正常開門的第一組數(shù)據(jù)有人聲干擾，故舍去。然后用plot函數(shù)畫出聲音波形圖，如下圖2所示： 根據(jù)所分析的參數(shù)類型，語音信號分析可以分成時域分析和變換域(頻域、倒譜域)分析。其中時域分析方法是最簡單、最直觀的

6、方法，因為它直接對語音信號的時域波形進(jìn)行分析。接下來我們先進(jìn)行時域分析。圖2 正常開門聲與非正常開門聲均值短時平均能量：定義n時刻某語音信號的短時平均能量En為： （1）式中，N為窗長，可見短時能量為一幀樣點值的平方和。一般我們認(rèn)為聲音在10-30ms之內(nèi)是穩(wěn)定的，取幀長也在10-30ms之內(nèi)，而幀移通常取5-15ms之間，所以取N=55、95、125、165。如下圖3、4為正常開門及非正常開門N取不同值時短時能量函數(shù)隨幀數(shù)的變化曲線，其中橫坐標(biāo)為幀數(shù)。由圖3、圖4可以看出，N=55,N=95時的曲線不夠平滑，而N=165的曲線又過于平滑，故選取N=125時的曲線。圖3 正常開門聲音的短時能量

7、曲線圖4 非正常開門聲音的短時能量曲線通過觀察短時能量曲線，可以看出正常開門時的能量比較集中且數(shù)值小，非正常開門時的能量比較分散且數(shù)值大。容易想到，這些現(xiàn)象與前面問題一的結(jié)論是相吻合的。短時平均過零率：短時平均過零率是指每幀內(nèi)信號通過零值的次數(shù)。對有時間橫軸的連續(xù)聲音信號，可以觀察到聲音的時域波形通過橫軸的情況。在離散時間聲音信號情況下，如果相鄰的采樣具有不同的代數(shù)符號就稱為發(fā)生了過零，因此可以計算過零的次數(shù)。 （2）上式為短時平均過零率的公式，其中，sgn為符號函數(shù)，即 （3）短時平均過零率曲線如下圖5、圖6所示：圖5 正常開門聲音的短時平均過零率圖6 非正常開門聲音的短時平均過零率由圖5、

8、圖6可以看出前者比后者的過零率要高。短時自相關(guān)函數(shù)：自相關(guān)函數(shù)用于衡量信號自身時間波形的相似性。對于聲音來說，采用短時分析方法，可以定義短時自相關(guān)函數(shù)為： （4）短時自相關(guān)函數(shù)曲線如下圖7、圖8所示：圖7 正常開門聲音的自相關(guān)函數(shù)曲線圖8 非正常開門聲音的自相關(guān)函數(shù)頻域特征：頻域分析主要是對聲音波形進(jìn)行傅里葉變換，如下圖9所示：圖9 聲音波形的傅里葉變換由圖9可以看出正常開門的對數(shù)幅度要比非正常開門的對數(shù)幅度小的多。程序源代碼見附錄二。問題三建立聲音識別模型（二分類模型），利用模型區(qū)分正常和非正常聲音，評價模型的好壞。問題三利用支持向量機(jī)構(gòu)建二分類器對正常與非正常敲門進(jìn)行分類預(yù)測。支持向量機(jī)是

9、Cortes和Vapnik于1995年提出的，與傳統(tǒng)分類器比較，它在解決小樣本、非線性及高維模式識別中表現(xiàn)出許多特有的優(yōu)勢，并能夠推廣應(yīng)用到函數(shù)擬合等其他機(jī)器學(xué)習(xí)問題中。支持向量機(jī)方法是建立在統(tǒng)計學(xué)習(xí)理論的VC 維理論和結(jié)構(gòu)風(fēng)險最小原理基礎(chǔ)上的，根據(jù)有限的樣本信息在模型的復(fù)雜性和學(xué)習(xí)能力之間尋求最佳折衷，以期獲得最好的推廣能力。圖10 支持向量機(jī)分類器結(jié)構(gòu)圖采集樣本共有80組音頻信號，從中選擇60組數(shù)據(jù)做為訓(xùn)練數(shù)據(jù)訓(xùn)練網(wǎng)絡(luò)，20組數(shù)據(jù)作為測試數(shù)據(jù)測試網(wǎng)絡(luò)分類能力。下面是準(zhǔn)備數(shù)據(jù)步驟：Step1: 根據(jù)給定的數(shù)據(jù)，選定訓(xùn)練集和測試集；Step2: 為訓(xùn)練集與測試集選定標(biāo)簽集；Step3: 利用訓(xùn)

10、練集進(jìn)行訓(xùn)練分類器得到model；Step4: 根據(jù)model，對測試集進(jìn)行測試集得到accuracy rate；訓(xùn)練集：trainset(); 分別取正常開門與非正常開門數(shù)據(jù)的一半作為訓(xùn)練集；測試集：testset(); 分別取正常開門與非正常開門數(shù)據(jù)的另一半作為測試集；標(biāo)簽集：label(); 取bedroom的數(shù)據(jù)為正類標(biāo)簽為1；forse的數(shù)據(jù)為負(fù)類標(biāo)簽為-1.下面為svm二分類器的分類結(jié)果：圖11 基于支持向量機(jī)的分類器上圖將展示了訓(xùn)練與測試的結(jié)果，結(jié)果等于1表示正常聲音，結(jié)果為-1為非正常開門聲，通過學(xué)習(xí)訓(xùn)練得到上圖結(jié)果。程序源代碼見附件三。通過仿真，驗證了支持向量機(jī)在解決分類問題

11、中的優(yōu)勢，今后可以將這方法用于更深層次的語音識別中。問題四試?yán)锰卣鬟x擇或變換，對特征向量進(jìn)行優(yōu)化，并利用參數(shù)優(yōu)化技術(shù)優(yōu)化模型的參數(shù)，使識別模型的準(zhǔn)確率提高。問題四其中支持向量機(jī)的核函數(shù)主要有： 線性核函數(shù)，多項式核函數(shù) ，徑向基核函數(shù)，兩層感知器核函數(shù)。不同的核函數(shù)直接影響SVM的分類能力，使用不同的核函數(shù)SVM分類的效果不同，所需要的時間也不同，下圖是使用短時能量值和均值為特征參數(shù)，選取RBF為核函數(shù)的分類效果圖：圖 12 SVMRBF分類圖圖13 分類模型錯誤率此模型相比于問題三中建立的SVM分類模型分類正確率明顯提高，達(dá)到了85.65%，可見合適的核函數(shù)選取的重要意

12、義。程序源代碼見附錄四。問題五白噪聲，是指功率譜密度在整個頻域內(nèi)均勻分布的噪聲，即其功率譜密度為 (5)式中，為常數(shù)，。由于1和為一對傅里葉變換對，所以白噪聲的自相關(guān)函數(shù)為 (6)由式(1)、式(2)可知，白噪聲的自相關(guān)函數(shù)僅在時才不為零,而對于其他任意的，白噪聲的自相關(guān)函數(shù)都為零，即在任意兩個不同時刻上的隨機(jī)變量都是不相關(guān)的。對于一個線性時不變系統(tǒng)(如濾波器)，該系統(tǒng)是由它的沖激響應(yīng)或等效地由它的頻率響應(yīng)表征，這里的和是一對傅里葉變換對。如果令為系統(tǒng)的輸人信號，是輸出信號。系統(tǒng)的輸出可以表示成如下形式： (7)如果是平穩(wěn)隨機(jī)過程的樣本函數(shù)，那么就是隨機(jī)過程的樣本函數(shù)?？梢郧蟪鲚敵龅淖韵嚓P(guān)函數(shù)

13、是 (8)由于我們知道自相關(guān)函數(shù)和功率譜密度函數(shù)是一對傅里葉變換對，所以可以得到輸出過程的功率密度譜，即為相關(guān)函數(shù)的傅里葉變換： (9)由此可以看出，輸出信號的功率譜密度就是輸入信號的功率譜密度乘以系統(tǒng)的頻率響應(yīng)的模的平方。當(dāng)輸入隨機(jī)過程是白噪聲時，輸出隨機(jī)過程的自相關(guān)特性和功率密度譜將完全由系統(tǒng)的頻率響應(yīng)所決定。我們設(shè)計了一個低通濾波器，對加了白噪聲的聲音信號進(jìn)行濾波，即可去除白噪聲。低通濾波(Low-pass filter) 是一種過濾方式，規(guī)則為低頻信號能正常通過，而超過設(shè)定臨界值的高頻信號則被阻隔、減弱。但是阻隔、減弱的幅度則會依據(jù)不同的頻率以及不同的濾波程序（目的）而改變。在數(shù)字圖像

14、處理領(lǐng)域，從頻域看，低通濾波可以對圖像進(jìn)行平滑去噪處理。濾波器將這部分不需要的頻率濾除，就實現(xiàn)了去噪的效果。濾波前與濾波后的聲音波形對比如下圖所示：圖14 濾波前后信號波形對比實驗表明低通濾波具有良好的去噪能力，該方法對添加的噪聲處理后，與原始信號效果幾乎擬合。證明低通濾波對信號去噪是可行的。程序源代碼見附錄五附 錄附錄一：clc;clear for i=5:7 j=2*(i-4)-1; a='正' num2str(i); load(a) subplot(3,2,j); plot(y) title('正' num2str(i) xlabel('Time&

15、#39;); ylabel('Amplitude'); a='非' num2str(i); load(a) subplot(3,2,j+1) plot(y) title('非' num2str(i) xlabel('Time'); ylabel('Amplitude');endsaveas(gcf,'q1.jpg');附錄二：均值function fei=feipj()fei=zeros(55125,1);for i=1:40 a='非' num2str(i); load(a) fe

16、i=fei+y;endfei=fei/40;短時平均能量function energy=nengliang(y)s=fra(55,55,y) ; s2=s.2; %一幀內(nèi)各樣點的能量energy=sum(s2,2) ; %求一幀能量subplot(2,2,1) %定義畫圖數(shù)量和布局 plot(energy) %畫N=55時的語音能量圖xlabel('幀數(shù)') %橫坐標(biāo)ylabel('短時能量 E') %縱坐標(biāo)title('N=55') %曲線標(biāo)識 s=fra(95,95,y) ; s2=s.2;energy=sum(s2,2) ; subplot

17、(2,2,2) plot(energy) %畫N=95時的語音能量圖xlabel('幀數(shù)')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('幀數(shù)')ylabel('短時能量 E')title('N=125')s=fra(165,165,y) ; s2=s.2;energy=sum(s2,2);subplot

18、(2,2,4) plot(energy) %畫N=165時的語音能量圖xlabel('幀數(shù)')ylabel('短時能量 E')title('N=165')function f=fra(len,inc,x)fh=fix(size(x,1)-len)/inc)+1);f=zeros(fh,len);i=1;n=1;while i<=fh j=1; while j<=len f(i,j)=x(n); j=j+1;n=n+1; end n=n-len+inc; i=i+1;end短時平均過零率function f=zcro(x)f=zeros

19、(size(x,1),1); %生成全零矩陣for i=1:size(x,1) z=x(i,:); %提取一行數(shù)據(jù) for j=1:(length(z)-1); if z(j)*z(j+1)<0; f(i)=f(i)+1; end endendfunction zcr=guoling(x)s=fra(125,125,x);%分幀，幀移110zcr=zcro(s);%求過零率figure(1);subplot(2,1,1)plot(x);title('非');xlabel('Time');ylabel('Amplitude');subplot

20、(2,1,2)plot(zcr);xlabel('幀數(shù)');ylabel('過零次數(shù)');title('原始信號的過零率');短時自相關(guān)函數(shù)function A1=zixiangguan(s1)N=125;A=; %加N=125的矩形窗for k=1:125;sum=0;for m=1:N-k+1;sum=sum+s1(m)*s1(m+k-1); %計算自相關(guān)endA(k)=sum;endfor k=1:125A1(k)=A(k)/A(1); %歸一化A(k);endplot(A1);title('非正常開門的自相關(guān)函數(shù)')xla

21、bel('延時 k')ylabel('R(k)')傅里葉變換clc;clearFs=11025;zheng=zhengpj();fei=feipj();T=1/Fs;L=55125;t=(0:L-1)*T;NFFT=2nextpow2(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 =

22、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;clearzhn=zeros(40,55125);zhn(1,:)=zhengpj()'for i=2:40 name='正',num2str(i); load(name) zhn(i,:)=y'endfei=zeros(40,5512

23、5);for i=1:40 name='非',num2str(i); load(name) fei(i,:)=y'endlabel=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_d

24、ataset',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);for i=1:20 train(1,:)=trainset(i,:)+train(1,:); train(2,:)=trainset(i+20,:)+train(2,:);end for i=1:20 test(1,:)=test

25、set(i,:)+test(1,:); test(2,:)=testset(i+20,:)+test(2,:);endtrain1=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)

26、;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.TestSizesaveas(gcf,'wentisan2.jpg');附錄四：clear;clc;closeall;disp('-采用徑向基內(nèi)積函數(shù)的支持向量機(jī)的應(yīng)用-');disp('輸入樣本矩陣(每行表示一個測試數(shù)據(jù)):');%x=randn(1,20)-3randn(1,20)+3;randn(1,40)'loadpj1.matloadpj2.matdisp('樣本所屬類別:');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=%gn',T);fprintf