基于MCRAOMLSA的语音降噪(三)：实现(续)

上篇文章（基于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。

下面给出软件实现时的步骤：

1)     定义结构体如下，其中s放符号位，e放阶码，m放尾码，dat是自然指数运算的返回值。

typedef union {

           double dat;

           struct{

                     unsigned long m:52;

                     unsigned e:11;

                     unsigned s:1;

           }jw;

}FREXP;

2)     求符号位和阶码。因为自然指数均大于0，所以符号位均为0。对ax取整加1023就可得阶码。

3)     求尾码。通过查表法求e^x₀，通过幂级数展开式求e^Δx，p = e^x₀e^Δx即可求得，α = p – 1也可求得。尾码 =α*2⁵²就得到了。

4)     返回dat值就是自然指数的结果了。

 

3，求自然对数

自然对数是自然指数的逆运算，y = e^x，根据y求x。给定一个y，通过上面定义的结构体FREXP可以得到它的阶码E和尾码M(M_i表示每一位上的值)，表示如下式：

又 y = e^x，所以

上式两边取自然对数，得到：

即：

其中

b好求，主要看a怎么求。a是ln(1 + x)的形式，泰勒展开式如下：

所以可以用泰勒展开式求b。a和b都求出来了，一个数y的自然对数x = a + b就求出来了。

下面给出软件实现时的步骤（与自然指数共用结构体）：

1)     得到阶码E和尾码M_i

2)     根据阶码E求得b

3)     根据尾码M_i利用泰勒展开式得到a

4)     将a和b相加就得到自然对数