基於MCRA-OMLSA的語音降噪(三)：實現(續)

2022 年 1 月 5 日
筆記
傳統音頻

上篇文章（基於MCRA-OMLSA的語音降噪(二)：實現）講了基於MCRA-OMLSA的語音降噪的軟體實現。本篇繼續講，主要講C語言下怎麼對數學庫里的求平方根（sqrt()）、求自然指數（exp()）、求自然對數（log()）的函數做替換。

1，求平方根

求平方根最常用的方法是牛頓迭代法。下圖是y = f(x)的曲線，當f(x) =0時的值(α)就是該方程的根。

可以通過多次迭代逼近的方法求得這個根，原理如下：

任取一個x₀，這個值對應的y值為f(x₀)。在x₀處畫y = f(x)的切線，與x軸交點為x₁。根據斜率的定義，在x₀處的斜率如下：

又斜率是函數的一次導數f』(x₀)，所以

可求得

基於x₁再畫一條切線，運用上面的求法得到與x軸交點為x₂，一直迭代下去可得x_3，…….，x_n，x_n+1等，從而求得x_n+1與x_n的關係如下式：

這些值會向方程的根α無限逼近。當| x_n+1– x_n| < ε (ε是事先設定的一個精度)時就停止迭代，這時x_n+1就是方程f(x) = 0的根。

具體到求平方根，x² = v (v是一個大於等於0的實數值)，x² – v = 0，令f(x) = x² – v ，得到f』(x) = 2x，把f(x)和f』(x)帶入上式得到

處理後得到

上式就是求平方根的迭代數學表達式。設定好精度後就可求出平方根，與C數學庫的sqrt()結果比較，值是非常接近的。

2，求自然指數

求自然指數是基於論文《指數函數e^x的快速計算方法》。用這個方法前得搞清楚浮點數的二進位存儲表示方法，浮點數包括單精度浮點數（float）和雙精度浮點數（double）。先看float的二進位存儲表示，float的搞明白了，double的類似，也好懂。

float佔4個位元組，32比特，存儲格式如下圖：

其中第0-22位共23位表示尾碼M，第23-31位共8位表示階碼E，第31位共1位表示符號位S。符號位好理解，0表示正數，1表示負數。以0.625為例，是正數，所以符號位是0。至於階碼和尾碼，方便理解，依舊以0.625為例。0.625 = 1.25 * 2^-1= (1 + 0.25) * 2^-1= (1 + x) * 2^y，其中x表示小數部分，y表示指數。

階碼E = y + 127 的二進位表示。這裡y = -1，所以E = -1 + 127 = 126，表示成二進位就是1111110，用8位二進位表示就是01111110。

尾碼M = x * 2²³的二進位表示。這裡x = 0.25，所以0.25 * 2²³= 2097152，用23位的二進位表示，M = 01000000000000000000000。

最終0.625的二進位存儲表示如下圖：

double佔8個位元組，64比特，存儲格式如下圖：

它的二進位表示跟float類似，不同的是階碼E = y + 1023。依舊以0.625為例,

階碼E = -1 + 1023 = 1022，表示成二進位就是1111111110，用11位二進位表示就是01111111110。

尾碼M = x * 2⁵²的二進位表示。這裡x = 0.25，所以0.25 * 2⁵²= 1125899906842624，用52位的二進位表示，M = 0100000000000000000000000000000000000000000000000000。符號位還是0。最終0.625的二進位存儲表示如下圖：

浮點數的存儲機制搞明白了，現在看怎麼求自然指數。求自然指數的傳統方法是用指數函數的冪級數展開式，如下式：

該論文用了一種計算速度更快的方法。下面具體看怎麼做的。為簡單起見，令x > 0，當x < 0時，只要用1除就可以了。

令 y = e^x，所以。log₂e是個定值1.4426950408889634，這裡令為a，即a = log₂e = 1.4426950408889634。從而log₂y = ax，即 y = 2^ax。令n是ax的整數部分，即 n = [ax]，從而ax的小數部分為ax – n，令其為D，即D = ax – n。所以 ax = n + D，y = 2^ax = 2^n+D = 2^D2ⁿ 。因為 0 < D < 1，所以1 < 2^D < 2，從而可以寫成1 + α（0 < α < 1）的形式，所以 y = (1 + α)2ⁿ。對標C數學庫里exp()用的是double型，這裡也用double型。根據上文double型的二進位存儲形式，可知n+1023就是階碼，α*2⁵²就是尾碼。n很好求，ax取整就可以了。下面看α怎麼求。α = 2^D – 1，2^D求出，α就有了。

令p = 2^D，從而。令x₀₀ = Dln2，有p = e^x₀₀。因為 0 < D < 1，又ln2 = 0.69314718056，所以 0 < x₀₀ < 0.69314718056。此時若直接用e^x₀₀的冪級數展開式求p，計算時間還很長，若適當選取x₀和Δx，使得Δx << 1，且 x₀₀ = x₀ + Δx，則有 p = e^{x₀ + Δx} = e^x₀e^Δx。可分別求e^x₀和e^Δx，然後再相乘就得到p。論文中用查表法求e^x₀，用冪級數展開法求e^Δx。先看怎麼求e^x₀。將x₀₀轉換為16進位數表示，改寫成x₀₀ = 0.q₁q₂q₃q₄q₅n = 0.q₁q₂q₃ + 0.000q₄q₅n = x₀ + Δx，其中x₀ = 0.q₁q₂q₃ = q₁ * 16^-1 + q₂ * 16^-2 + q₃ * 16^-3，Δx = 0.000q₄q₅n = q₄ * 16^-4 + q₅ * 16^-5 + …。所以e^x₀ = e^{q₁ * 16^-1 + q₂ * 16^-2 + q₃ * 16^-3}= e^{q₁ * 16^-1}e^{q₂ * 16^-2}e^{q₃ * 16^-3}。因為x₀ < x₀₀ < 0.69314718056 < 0.75 = 12/16，所以q₁的取值範圍是[0, 11]，q₂的取值範圍是[0, 15]，q₃的取值範圍是[0, 15]。根據q_x的有限個不同取值將e^{q₁ * 16^-1} 、e^{q₂ * 16^-2} 和e^{q₃ * 16^-3} 分別預先算出做成表，計算時通過查表得到三個相應的值，再將這三個值相乘就得到e^x₀的值了。再來看怎麼求e^Δx。0 < Δx = 0.000q₄q₅n < 16^-3 = 1/4096 << 1，用冪級數展開式求e^Δx只要取前面4項即可保證精度了，所以用冪級數展開式求e^Δx。