NullPointerExceptions

來為久違的痞客邦拔草，

順便記錄上學期(四下)的一堂課：ASAS數位訊號分析與合成，的最後一份作業。

這是一堂好課，

而這是關於 Linear Prediction Modeling (下稱LPM)*，

與它的演算法：Levinson-Durbin recursion method (下稱LDRM)* 的一份作業。

(*詳見最下方資料補充)

而在此作業中，除了解讀老師提供的code，

我稍微討論了程式內各種參數，對於結果的影響。

這就是本篇文章的目的，記錄我的心得與討論結果。

=========正題前的acoustic科普(?)=========

先假設發出聲音的過程，都可以用impedance(Z，猶如電學的電阻、超音波中的attenuation)來描述。

它是一個與截面積成反比的係數。

以人來說，氣流從肺部迅速地通過氣管、聲帶、喉嚨、口腔，向外吐出發聲。

這一路上的凹凸起伏，都可以看作是截面積不同的軟橡皮管，組裝而成。

(在形式上，是假設為全pole無zero的IIR filter*)

而軟管之間，氣流因為管子粗細不同會產生通過/反彈，

兩者之比名為 reflectance，

以impedance來表示即為： r12 = (Z1-Z2)/(Z1+Z2)；

相等的另一種表示法就是：Z1/Z2 = (1+r12)/(1-r12)

(此關係式之後會用到)

LPM的目的，

就是從已知的聲音資訊，拼湊出發聲軟管的原始面貌。

得知之後，通入氣流，自然就能合成出聲音了。

而方法，簡單來說，LPM就是收集聲音資訊，用 Least-Square Formulation(最小平方法)*，

以逼近得到

(1) 發聲源(source)的反彈係數

(2) 發聲管道每一截的impedance

(3) 管道末端(load)的反彈係數

而 LSF 詳細作法如何，並不是很重要，

因為matlab已經有了lpc()這個函式，短短一行就能得到error值最小的解。

…

但這不太好XD

LSF 是一種O(n^3)的苦勞，計算反矩陣，即使是電腦也是很辛苦的。

於是LDRM因應而生，用空間(記憶體)換取時間的作法，是以雙重迴圈實踐O(n^2)的方法。

親愛的matlab也支援了這個方法：LEVINSON(R,p)

R是音訊平方；

p是LPM的order，一個隨便你設的參數，後面會提到它的影響。

有了以上的基礎，就能來討論一下matlab code中，各參數對於模擬結果的影響了:D

=========科普(?)結束=========

==========實踐步驟==========

以matlab實踐LPM的步驟如下：

(1) 讀取音檔 (本份作業是要自己錄五個母音)

(2) 設定 fs 重新sampling音檔(或用matlab預設的44100)

(3) 設定 frame_length。時間上連續的音檔，必須切成分開的frame才能進行分析嘛

(4) 設定p，linear prediction order，意義上是「把發聲軟管假設為幾截」，也就是IIR filter中的pole數。

(5) 是否要作 pre-emphasis

仔細一數，能設的也就這四個參數了，接下來就是一些固定的數值。

比如：聲速c、繪圖用的dx、每個frame能得到的equation數量：L等等，因為是受上面參數影響的定值，就不討論了。

不過它們會給參數帶來一些限制，留下一階段討論。

(6) 接下來把音檔做一些處理，看是乘個window、調整長度、取平方之類的，連同p一同丟給 Levinson()函式即可。

      打得好像把食物丟到烤箱一樣

(7) 觀察結果，要看頻率、error、波形等等都可。

      Levinson丟出來的A、E、K，分別為

      A=LPM coefficient (就是IIR filter的乘項)

      E=各階段計算出來的差值

      K=可以理解為A型態變化後的結果，是LDRM特別產生、便於計算A 的數列。

      不過非常剛好地(?)，K值正是最前面假設中的reflectance。

      而在第一段提到的 Z1/Z2 = (1+r12)/(1-r12)，就可以派上用場了，

      由於Z和截面積成反比，將Z開根號，即是和截面積半徑成反比。

      於是，當我們有一堆reflectance(K)時，就能得知每截水管的截面積半徑之比!

      如此累積相乘的結果，便能觀察發聲構造的幾何分佈了XD

      不需要解剖手術，就能像用手捏著一般量測身體裡喉管的寬度，

      是不是很有趣呢(巧虎調)

==========步驟結束==========

==========參數影響==========

(1) fs

首先是 sampling frequency：fs

使用matlab預設的44100也沒什麼不好，就是計算量大了點。

(我們還是盡量避免這點吧)

fs的下限，根據可愛的Nyquist Theorem*，fs至少要fc兩倍，

就看你音檔的fc 如何了，我對自己的音檔取了-30dB的頻率範圍，fc都在5000以下，

所以fs取了10000Hz。

不同fs沒有影響，只要它在Nyquist Theorem的包蔽下。

雖然步驟很eazy，但這個fs，將是一切的起點。

***

(2) P-the Linear Prediction Order

基本上，p這個數字，決定了我們模擬出來的真實度。

就像微分切得愈細愈精準，p值代表變數的數量、計算量與不同寬度管子銜接的數量。

但它並不能無限上漲。

還記得步驟的最後，我們要繪製一個人體喉道表嗎?

人體喉道的長度是有限制的，自己摸著量量能知道，大概35~40cm吧(包含轉折長度)

如果最後繪出來，長度突破天際，就……

於是這個隨著fs而來的限制，誕生了：

      dx = 聲速/fs;   % spatial, unit:cm per sample

      代入

      35(cm)<=p*dx(喉道總長)=<40(cm)

得到：0.0001020<=(p/fs)<=0.0001166

如果用matlab預設fs=44100，則這就是p的上限：

      p<=51.4

(也可以說是matlab執行Linear Prediction的極限……)

如果用我自訂的fs=10000，則下限為：

     10.2<=p 

進整數取11；

這就是尾隨而來的第二個限制。

而關於不同P造成的影響，繪製頻譜圖可發現，

P愈大時，模擬愈接近真實情況，所以頻率曲線變得毛躁不平滑……

是好是壞就因應用而異了。

作業內容是要找出共振峰值f1、f2，所以我使用了P=11(最小值)的頻率曲線，效果十分顯著。

反觀用P=50的話，就不是那麼直觀了。

而在error方面，隨著p增大，error變少的可以預期的，成為趨近於0的一直線。

繪製喉道表的部份，是此參數的重點，當P=11時，長度與半徑都落在合理範圍內；

而P=50時，長度長達170cm、半徑寬達20、30，明顯……不是人類(也許外星來的巨嘴長脖怪物吧:D

***

(3) frame_length

相較於上述，frame_length的限制比較顯而易見。

首先，frame_length不能比音檔總長度長……(不到一個frame是要做什麼呢)

第二，LPM的計算方式中，L是等式數量、P是變數數量，

矩陣中想要解出所有的變數……等式數量必須比較多：

      L>=P

而   L=frame_length*fs

(一個frame中能產生幾個等式)

所以得到：

      frame_length>=p/fs  

不過可惜的是，frame_length對於頻率曲線的影響、波形、error都不大，

只是影響顯示的間距，非常無用。

唯一值得一提的是，在繪製喉道表時，

因為frame_length影響了迴圈的次數，產生的reflectance數量不同。

frame_length較短，疊乘出來，讓半徑的寬度變寬了 。

真是誤導人的功能啊

***

(4) Pre-Emphasis

為什麼要做pre-emphasis?

因為聲音的高頻頻率的振幅會下降，為了方便觀察高頻部份，

處理音訊時，常會特別加強其高頻部份。

但這麼做是有風險的。

基本上改變音訊的本質，會直接造成結果誤差。

從喉道表能看出，在/a/、/e/、/o/的部份，由於此三音的高頻部份較低，

經過pre-emphasis後，需要口腔張得特別寬，才能發出具有強烈高頻的聲音:D

==========參數結束==========

===========結論============

上面打的大概只是1/4部份吧，畢竟這份作業我寫了12頁。

不過我懶了，何況這篇是記錄性質，大概沒什麼人願意讀完吧，HAHA。

要感謝劉老師開了這堂非常充實又愉快的ASAS課程，

深入淺出、導讀paper、練習matlab code等等，都是我進化史中不可或缺的基石!

很可惜(?)的是，直到最後一份作業(此份)，

才從助教口中得知，應該針對各種係數進行討論，而不是單純回答作業問題。

而從討論中，了解到code與LPM、LDRM的實質意義，是很寶貴的。

如果有興趣了解作業題目、內容、講義，請連絡我。

上述有什麼錯誤，也請留言訂正，感謝:D

==========================

不是很重要的資料補充：

Linear Prediction Modeling：

http://en.wikipedia.org/wiki/Linear_prediction

Levinson-Durbin method：

http://en.wikipedia.org/wiki/Levinson_recursion

IIR：

http://en.wikipedia.org/wiki/Infinite_impulse_response

Least-Square formulation：

http://en.wikipedia.org/wiki/Least_squares

Nyquist Theorem：

http://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem