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

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

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

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

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

6、方法，因為它直接對語音信號的時域波形進行分析。接下來我們先進行時域分析。圖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 正常開門聲音的短時能量

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

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

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

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

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

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

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

14、處理領域，從頻域看，低通濾波可以對圖像進行平滑去噪處理。濾波器將這部分不需要的頻率濾除，就實現了去噪的效果。濾波前與濾波后的聲音波形對比如下圖所示：圖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; %一幀內各樣點的能量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

17、(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

18、(2,2,4) plot(energy) %畫N=165時的語音能量圖xlabel('幀數')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,:); %提取一行數據 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('幀數');ylabel('過零次數');title('原始信號的過零率');短時自相關函數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); %計算自相關endA(k)=sum;endfor k=1:125A1(k)=A(k)/A(1); %歸一化A(k);endplot(A1);title('非正常開門的自相關函數')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('-采用徑向基內積函數的支持向量機的應用-');disp('輸入樣本矩陣(每行表示一個測試數據):');%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