1、目目 錄錄 1 引言.1 2 課程設計內容及要求.2 2.1 課程設計內容 .2 2.2 課程設計要求 .2 2.3 課程設計目的 .2 2.3 課程設計平臺.2 3 基于 matlab 的 fft 算法設計原理.3 3.1 總體設計流程圖 .3 3.2 語音信號的采集 .3 3.3 語音信號的時頻分析 .3 3.4 快速傅里葉變換 .6 3.4.1 fft的運算規律.8 3.4.2 基于 matlab 的 fft 所編寫程序的框圖 .12 3.5 自編算法與機帶算法仿真波形比較 .13 4 設計總結.16 參考文獻.17 附錄.18 課程設計說明書課程設計說明書 1 1 引言 隨著信息時代，

2、數字時代的到來，數字信號處理已經成為一門極其重要的學 科和技術領域。以 dsp 為核心芯片的處理系統日益變成了數字信號處理系統的主 流。它廣泛用于電子信息、通信、圖像處理、語音處理、生物醫學、自動控制、 地質探測等領域，受到工程設計和使用人員的青睞。 matlab，它是美國 math works 公司推出的一種面向工程和科學計算的交互 式計算軟件。它以矩陣運算為基礎，把計算、可視化、程序設計融合在一個簡單 易用的交互式工作環境中，是一款數據分析和處理功能都非常強大的工程適用軟 件。通過本次實習我們學會了分析和處理音頻信號，首先要對聲音信號進行采集， matlab 的數據采集工具箱提供了一整套命

3、令和函數，通過調用這些函數和命令， 可直接控制聲卡進行數據采集。window 自帶的錄音機程序也可驅動聲卡來采集語 音信號，并能保存為 wav 格式文件，供 matlab 相關函數直接讀取、寫入或播 放。 matlab 語言是一種數據分析和處理功能十分強大的計算機應用軟件，它可 以將聲音文件變換位離散的數據文件，然后利用其強大的矩陣運算能力處理數據， 如數據濾波、傅立葉變換、時域和頻域分析、聲音回放以及各種圖的呈現等，它 的信號處理與分析工具箱位語音信號分析提供了十分豐富的功能函數，利用這些 功能函數可以快捷而又方便的完成語音信號的處理和分析以及信號的可視化，是 人機交互更加便捷。信號處理是

4、matlab 重要應用的領域之一。 語音信號處理是研究用數字信號處理技術和語音學知識對語音信號進行處理 的新興的學科，是目前發展最為迅速的信息科學研究領域的核心技術之一。通過 語音傳遞信息是人類最重要、最有效、最常用和最方便的交換信息形式。 語音信號的處理與濾波的設計主要是用 matlab 作為工具平臺，設計中涉 及到聲音的錄制、播放、存儲和讀取，語音信號的抽樣、頻譜分析，濾波器的設 計及語音信號的濾波，通過數字信號處理課程的理論知識的綜合運用。從實踐上 初步實現對數字信號的處理。 課程設計說明書課程設計說明書 2 2 課程設計內容及要求 2.1 課程設計內容 錄制一段個人自己的語音信號，并對

5、錄制的信號進行采樣；畫出采樣后語音 信號的時域波形和頻譜圖；在 matlab 環境下編寫基 2 dit-fft 算法；利用自己編 寫的算法對已采集的語音信號進行頻譜分析，并畫出語音信號的時域與頻譜圖， 并與 matlab 數字信號處理工具箱中的 fft 函數進行對比研究，驗證自編算法的正確 性。 2.2 課程設計要求 1.完成語音信號的采集，利用 windows 自帶的錄音機或其他軟件，錄制一段 語音，時間在 1s 以內； 2.在 matlab 中編寫程序，實現輸入信號的倒序； 3.編寫程序，實現蝶形運算； 4.畫出語音信號的頻譜圖，與 matlab 數字信號處理工具箱中的 fft 函數進行對

6、 比研究，并對設計結果進行獨立思考和分析； 2.3 課程設計目的 1.學會 matlab 的使用，掌握 matlab 的程序設計方法。 2.掌握在 windows 環境下語音信號采集的方法。 3.掌握數字信號處理的基本概念、基本理論和基本方法。 4.掌握 matlab 設計 fir 和 iir 數字濾波器的方法。 5.學會用 matlab 對信號進行分析和處理。 2.4 課程設計平臺 matlab7.1 軟件 matlab 是由美國 mathworks 公司發布的主要面對科學計算、可視化以及交 互式程序設計的高科技計算環境。它將數值分析、矩陣計算、科學數據可視化以 及非線性動態系統的建模和仿真

7、等諸多強大功能集成在一個易于使用的視窗環境 中，為科學研究、工程設計以及必須進行有效數值計算的眾多科學領域提供了一 種全面的解決方案，并在很大程度上擺脫了傳統非交互式程序設計語言（如 c、fortran）的編輯模式，代表了當今國際科學計算軟件的先進水平。 課程設計說明書課程設計說明書 3 3 基于 matlab 的 fft 算法設計原理 3.1 總體設計流程圖 在一個相對較安靜的環境下，錄下 1s 左右的 wav 聲音信號，然后對聲音進行 采樣，畫出其時域波形和頻譜圖，其流程圖如圖 1 所示： 圖 1 設計流程圖 3.2 語音信號的采集 在實際工作中，我們可以利用 windows 自帶的錄音機

8、錄制語音文件。采集到 語音信號之后，需要對語音信號進行分析，如語音信號的時域分析、頻譜分析、 語譜圖分析。在 matlab 中，我們可以通過y,fs,bits=wavread(語音信號路徑 ,n1 n2)語句。用于讀取語音，采樣值放在向量 y 中，fs 表示采樣頻率(hz)，bits 表示采樣位數。n1 n2表示讀取從 n1 點到 n2 點的值（若只有一個 n 的點則表 示讀取前 n 點的采樣值）。向量 y 則就代表了一個信號（也即一個復雜的“函數表 達式”）也就是說可以像處理一個信號表達式一樣處理這個聲音信號。 3.3 語音信號的時頻分析 利用matlab中的“wavread”命令來讀入（采

9、集）語音信號，將它賦值給某一 開始 輸入聲音信號 對聲音信號采樣 蝶形運算 原信號 fft 運算 結束 課程設計說明書課程設計說明書 4 向量。再對其進行采樣，記住采樣頻率和采樣點數。下面介紹wavread 函數幾種 調用格式。 1.y=wavread（file） 功能說明：讀取 file 所規定的 wav 文件，返回采樣值放在向量 y 中。 2.y,fs,nbits=wavread(file) 功能說明：采樣值放在向量 y 中，fs 表示采樣頻率（hz），nbits 表示采樣位 數。 3.y=wavread（file，n） 功能說明：讀取錢 n 點的采樣值放在向量 y 中。 4.y=wavr

10、ead（file，n1,n2） 功能說明：讀取從 n1 到 n2 點的采樣值放在向量 y 中。 接下來，對語音信號 speech off.wav 進行采樣。其程序如下： y,fs,nbits=wavered (speech off.wav); 功能說明：把語音信號加載入 matlab 仿真軟件平臺中 然后，畫出語音信號的時域波形，再對語音信號進行頻譜分析。matlab 提供了 快速傅里葉變換算法 fft 計算 dft 的函數 fft,其調用格式如下： xk=fft(xn,n) 參數 xn 為被變換的時域序列向量，n 是 dft 變換區間長度，當 n 大于 xn 的長度 時，fft 函數自動在

11、xn 后面補零。，當 n 小于 xn 的長度時，fft 函數計算 xn 的前 n 個元素，忽略其后面的元素。 原始信號的時域波形圖如圖 3 所示： 課程設計說明書課程設計說明書 5 圖 3 原始信號的時域波形圖 原始信號的頻域特性圖如圖 4 所示： 課程設計說明書課程設計說明書 6 圖 4 原始信號的頻域特性圖 3.4 快速傅里葉變換 快速傅里葉變換（fft）是為提高 dft 運算速度而采用的一種算法。 對一個有限長度序列 x(n)的 n 點的 dft 為：x（k）=x（n）wknn (k=0,1,，n-1；n=0,1,，n-1;w=e-j2/n) 當 n=4 時，x（k）可展開為： x（0）

12、= x（0）w0*4+ x（1）w0*4 +x（2）w0*4+ x（3）w0*4 x（1）= x（0）w0*4+ x（1）w1*4 +x（2）w2*4+ x（3）w3*4 x（2）= x（0）w0*4+ x（1）w2*4 +x（2）w4*4+ x（3）w6*4 x（3）= x（0）w0*4+ x（1）w3*4 +x（2）w6*4+ x（3）w9*4 從上式可以看出,要求 4 點的 dft,需要 16 次的復數乘法運算,12 次復數乘法運 算算。由此類推，要求出 n 點的 dft，需要 n2 次復數乘法運算，n*(n-1)次復 數加法運算。當 n 值較大時，要完成的復數乘法運算和復數加法運算得次

13、數都非 常多，無論是用通用計算機還是用 dsp 芯片，都需要消耗大量的時間，不適合于 對實時處理要求高的場合。為了能實時處理 dft,要想減少 dft 的運算量可以有兩 課程設計說明書課程設計說明書 7 個途徑：第一是降 n，n 的值減小了，運算量就減少了；第二是利用旋轉因子的 周期性和對稱性，可約性。利用這兩個途徑實現 dft 的快速傅里葉變換（fft） ， fft 算法基本上可分為時域抽取法和頻域抽取法。 w=e-j2/n 的性質： （1）周期性 （2）共軛對稱性 （3）可約性 本程序是用基 2 的按時間抽取的 fft 算法（dit-fft） ，設序列 x(n)的長度為 n，且 n 滿足

14、n=2m,m 為正整數。若 n 不能滿足上述關系，可以將序列 x(n)補零 實現，則 x(n)的 n 點 dft 為： x（k）=x（n）wknn (k=0,1,，n-1；n=0,1,，n-1;w=e- j2/n)將 n 分為奇數與偶數兩部分。 按時間抽取基 2-fft 算法的基本思路是將 n 點序列按時間下標的奇偶分為兩 個 n/2 點序列，計算這兩個 n/2 點序列的 n/2 點 dft，計算量可減小約一半；每 一個 n/2 點序列按照同樣的劃分原則，可以劃分為兩個 n/4 點序列，最后，將原 序列劃分為多個 2 點序列，將計算量大大降低。 1. 按時間下標的奇偶將 n 點 x(n)分別抽

15、取組成兩個 n/2 點序列，分別記為 x1(n)和 x2(n)，將 x(n)的 dft 轉化為 x1(n)和 x2(n)的 dft 的計算。 )()(nnk n nnk n kn n www *)(*)( nk n nk n kn n www mkn mn kn n mkn mn kn n wwww / / , 1 2 , 1 , 0, )() 12( )()2( 2 1 n r rxrx rxrx 1 n 0 21 nn 0,2,4.1,3,5. 11 22 212 nn 0,10,1 11 22 212 1n2n 0,10,1 221 n nk n nn nknk nn nn rkrk r

16、r nn rkrk rr x kx n w x n wx n w xr wxrw x r wxr w 2 j 2 j2 2 2 2 ee rk n rk rkrk n nn ww 課程設計說明書課程設計說明書 8 用蝶形運算可表式為： 以此類推，還可以把 x1(n)和 x2(n)按 n 值得奇偶分為兩個序列，這樣就達到了 降 n 得目的，從而減少了運算量。 fft 對 dft 的數學運算量改進： 直接采用 dft 進行計算，運算量為 n2 次復數乘法和 n*(n-1)次復數乘法。 當采用 m 次 fft 時，由 n=2m 求得 m=logn，運算流圖有 m 級蝶形，每一 級都由 n/2 個蝶形

17、運算構成，這樣每一級蝶形運算都需要 n/2 次復數乘法和 n 次 復數加法。m 級運算共需要復數乘法次數為 c=n/2*m,復數加法次數為 c=n*m。 當 n 值較大時，fft 減少運算量的特點表現的越明顯。 3.4.1 fft 的運算規律 （1）原位運算 n=2m的fft共m級運算，每級有n/2蝶形原位計算，當數據輸入到存儲器 以后，每一組蝶形運算后，結果仍然存放在這同一組存儲器的同一位置，不需要 另辟存儲空間，直接最后輸出。 同一級的蝶形運算每個蝶形運算的輸入數據對其他級輸入沒有影響。 （2）倒序運算的規律 輸入序列先按自然順序存入存儲單元,然后經變址運算來實現倒位序排列,用 j 表示倒

18、序的十進制數，對 n=2m，m 位的二進制數從左到右各位數權值位 n/2，n/4，n/82，1。因此，最高位加 1 相當于 j+n/2。 .如果最高位為 0，則直接得到下一個倒序值，j+n/2； .如果最高位為 1，則最高位為 0(j-n/2)，次高位加 1(j+n/4)。 .以此類推，直到最后一位二進制數字。 例如 ，n=8 時如下圖 5 所示： 11 22 12 22 00 12 , 01 nn rkkrk nnn rr k n x kx r wwxr w x kw xkn （）（k） 課程設計說明書課程設計說明書 9 順序倒序 十進制二進制二進制十進制 00000000 10011004

19、 20100102 30111106 41000011 51011015 61100113 71111117 圖 5 碼位倒序（n=8） 倒序的流程圖如圖 3 所示： 課程設計說明書課程設計說明書 10 圖 3 倒序的流程圖 （3）蝶形運算的運算規律 設序列x（n）倒序后存放在數組a中，如果蝶形運算兩個輸入相距b=2(l-1). lh=n/2 n1=n-1 j=0 t=x(j+1) x(j+1)=x(i+1) x(i+1)=t i=0：n1 i=j k=lh jk j=j+k j=j-k k=k/2 y n y n 課程設計說明書課程設計說明書 11 采用原位運算：蝶形運算可表示為： x(k)

20、=x(k)+x(k+b)*wnp; x(k+b)=x(k)- x(k+b)*wnp.其中p=j*2（m-1）,(j=0,1,2,3,2(l-1)-1) 旋轉因子的確定： 第l級共有2(l-1)個旋轉因子，以(n=23=8為例)： l=1時，w=w j=0； p n j l2 l=2時，w=w j=0，1； p n j l2 l=3時，w=w j=0，1，2，3； p n j l2 w=w=w p=j*2(m-l)，用來確定第l級旋轉因子 p n j l2 )(2*lmj n 從輸入端開始，共進行 m 級運算，在進行第 l 級運算時，依次求出 2(l-1)個 旋轉因子，然后計算每個旋轉因子所對應

21、的 2(m-l)個蝶形元素。第 l 級的蝶形運 算中，每個蝶形運算的兩個輸入相距 b=2l-1,同一旋轉因子對應的蝶形運算相隔 2l 個，同一旋轉因子對應的蝶形運算有 2m-l 個。 注：1：控制第 l 級順序運算 2：控制不同種的旋轉因子 3：控制同種旋轉因子所對應的蝶形運算 蝶形運算的流程圖如圖 4 所示： ： 課程設計說明書課程設計說明書 12 1 2 3 圖 4 蝶形運算的流程圖 3.4.2 基于 matlab 的 fft 所編寫程序的框圖 基于 matlab 的 fft 所編寫程序的框圖如圖 5 所示： 開始開始開始 j=0：b-1 n=2m b=2l/2 l=1:m 輸入 x(n)

22、 k=j+1:2l:n p=j*2(m-l) t=x(k)+x(k+b)* wnp x(k+b)=x(k)-x(k+b)* wnp x(k)=t 倒序 結束 輸出 課程設計說明書課程設計說明書 13 圖 5 程序框圖 3.5 自編算法與機帶算法仿真波形比較 我們知道 matlab 軟件自帶 fft 算法，我們可以通過比較自編算法仿真結果與 機帶算法仿真結果來檢驗自編算法的正確性。 自編算法與 fft 算法幅值比較圖如圖 6 所示： 語音信號采集 完 成 語 音 信 號 時 域 圖 完 成 語 音 信 號 頻 率 特 性 圖 編寫 fft 程序，畫出語音信號 頻譜圖 實 現 輸 入 信 號 的

23、倒 序 實 現 一 級 中 不 同 種 蝶 形 算 運 實 現 一 級 中 相 同 種 蝶 形 運 算 與 m atl ab 自 帶 的 fft 比 較 課程設計說明書課程設計說明書 14 圖 6 幅值比較圖 自編算法 fft 算法頻譜分貝計較圖如圖 7 所示： 課程設計說明書課程設計說明書 15 圖 7 分貝比較圖 由以上兩圖可以看出，經過蝶形運算得出的頻譜圖和信號直接 fft 得出的頻 譜圖一致。該程序嚴格按程序框圖編寫，思路清晰、容易理解，程序的運行過程 在命令窗中一目了然。通過與 fft 函數運算的結果比對，程序編寫正確。 課程設計說明書課程設計說明書 16 4 設計總結 為期四天的課

24、設很快接近尾聲了，基于 matlab 的 fft 算法實現設計已按計劃 如期全部完成，通過這次 dsp 課程設計，我對課堂上所學到的理論知識的理解加 深了許多， 自己動腦、動手設計的能力也得到了較大提高。在這次課程設計的過 程中，我對 matlab 語言有了更深的認識。現在仔細想想，這次課程設計使得我 對 matlab 語言的理解與應用能力得到了較大的提升，也讓我認識到只要深入學 習，提升的空間永遠是存在的。在設計的過程中我遇到了一些問題，如：編寫源 程序中出現了語法錯誤等。通過查閱書本和以前設計的程序我發現了產生錯誤的 原因并解決了問題完成了設計。經過反思我發現較大一部分錯誤是因為操作的不

25、熟練造成的，這也讓我明白了要保持設計的高效率 必須經常練習。另一方面我也 發現了動手實踐的重要性。動手實踐是理論知識得以靈活運用的必要前提，也是 今后走上工作崗位之后能夠很好的完成設計工作的技術保證。只有遇到實際問題 并根據自己對課堂上獲得的專業知識的理解來解決才能真正的提高自己的能力。 這次程序設計讓我獲益匪淺，對 matlab 也有了進一步的認識：matlab 功能強、 使用靈活方便等。matlab 是在國內外廣泛使用的一種數字信號處理系統，相信除 了以上優點，還有許多我還未發現，希望能在以后的學習中有更深入的認識。 在此，我也要感謝一下董翠英老師對我這次課設的大力支持和指導，如果沒 有老

26、師的幫助我很難在這么短的時間內完成這次課設。 課程設計說明書課程設計說明書 17 參考文獻 1范壽康 dsp 技術與 dsp 芯片.北京：電子工業出版社 2程佩青.數字信號處理教程.北京：清華大學出版社出版，2001 3高西全 丁玉美等 數字信號處理. 北京：電子工業出版社，2009 4 余成波，陶紅艷。數字信號處理及 matlab 實現，北京：清華大學出版社， 2008 5 曹弋，趙陽。matlab 實用教程，北京：電子工業出版社，2007 課程設計說明書課程設計說明書 18 附錄 源程序： clear all; x,fs,bits=wavread(gl.wav,2048); %讀取聲音 x

27、1=reshape(x,1,4096); sound(x1,fs,bits); %播放語音信號 y1=fft(x1) n=length(x1); %計數讀取信號的點數 t=(1:n)/fs; %信號的時域采樣點 figure(1) plot(t, x1); %畫出聲音采樣后的時域波形 title(原聲音信號的時域波形); %給圖形加注標簽說明 xlabel(時間/t); ylabel(振幅/a); grid ; %添加網格 m=nextpow2(x1); % 求 x 的長度對應的 2 的最低冪次 m n=2m; if length(x1)n x1=x1,zeros(1,n-length(x1); % 若 x 的長度不是 2 的冪，補零到 2 的整數冪 end nv2=n/2; nm1=n-1; i=0; j=0; while inm1 if ij t=x1(j+1); x1(j+1)=x1(i+1); x1(i+1)=t; end k=nv2; 課程設計說明書課程設計說明書 19 while k=j j=j-k; k=k/2; end j=j+k; i=i+1;