基于MCRA-OMLSA的语音降噪(三)：实现(续)

上篇文章（基于MCRA-OMLSA的语音降噪(二)：实现 ）讲了基于MCRA-OMLSA的语音降噪的软件实现。本篇继续讲，主要讲C语言下怎么对数学库里的求平方根（sqrt()）、求自然指数（exp()）、求自然对数（log()）的函数做替换。

1，求平方根

求平方根最常用的方法是牛顿迭代法。下图是y = f(x)的曲线，当f(x) =0时的值(α)就是该方程的根。

可以通过多次迭代逼近的方法求得这个根，原理如下：

任取一个x0，这个值对应的y值为f(x0)。在x0处画y = f(x)的切线，与x轴交点为x1。根据斜率的定义，在x0处的斜率如下：

又斜率是函数的一次导数f’(x0)，所以

可求得

基于x1再画一条切线，运用上面的求法得到与x轴交点为x2，一直迭代下去可得x3，…….，xn，xn+1等，从而求得xn+1与xn的关系如下式：

这些值会向方程的根α无限逼近。当| xn+1 - xn| < ε (ε是事先设定的一个精度)时就停止迭代，这时xn+1就是方程f(x) = 0的根。

具体到求平方根， = v (v是一个大于等于0的实数值)， – v = 0，令f(x) =  – v ，得到f’(x) = 2x，把f(x)和f’(x)带入上式得到

处理后得到

上式就是求平方根的迭代数学表达式。设定好精度后就可求出平方根，与C数学库的sqrt()结果比较，值是非常接近的。

2，求自然指数

求自然指数是基于论文《指数函数的快速计算方法》。用这个方法前得搞清楚浮点数的二进制存储表示方法，浮点数包括单精度浮点数（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 *   = (1 + 0.25) *  = (1 + x) * ，其中x表示小数部分，y表示指数。

阶码E = y + 127 的二进制表示。这里y = -1，所以E = -1 + 127 = 126，表示成二进制就是1111110，用8位二进制表示就是01111110。

尾码M = x * 的二进制表示。这里x = 0.25，所以0.25 * = 2097152，用23位的二进制表示，M = 01000000000000000000000。

最终0.625的二进制存储表示如下图：

double占8个字节，64比特，存储格式如下图：

它的二进制表示跟float类似，不同的是阶码E = y + 1023。依旧以0.625为例,

阶码E = -1 + 1023 = 1022，表示成二进制就是1111111110，用11位二进制表示就是01111111110。

尾码M = x * 的二进制表示。这里x = 0.25，所以0.25 * = 1125899906842624，用52位的二进制表示，M = 0100000000000000000000000000000000000000000000000000。符号位还是0。最终0.625的二进制存储表示如下图：

浮点数的存储机制搞明白了，现在看怎么求自然指数。求自然指数的传统方法是用指数函数的幂级数展开式，如下式：

该论文用了一种计算速度更快的方法。下面具体看怎么做的。为简单起见，令x > 0，当x < 0时

 ，只要用1除就可以了。

令 y = ，所以

log2e是个定值1.4426950408889634，这里令为a，即a = log2e = 1.4426950408889634。从而log2y = ax，即 y = 。令n是ax的整数部分，即 n = [ax]，从而ax的小数部分为ax – n，令其为D，即D = ax – n。所以 ax = n + D，y = =  =  。因为 0 < D < 1，所以1 < < 2，从而可以写成1 + α（0 < α < 1）的形式，所以 y = (1 + α)。对标C数学库里exp()用的是double型，这里也用double型。根据上文double型的二进制存储形式，可知n+1023就是阶码，α*就是尾码。n很好求，ax取整就可以了。下面看α怎么求。α =  – 1，求出，α就有了。

令p = ，从而 

令 = Dln2，有p = 。因为 0 < D < 1，又ln2 = 0.69314718056，所以 0 <  < 0.69314718056。此时若直接用的幂级数展开式求p，计算时间还很长，若适当选取和Δx，使得Δx << 1，且  =  + Δx，则有 p =  = 。可分别求和，然后再相乘就得到p。论文中用查表法求，用幂级数展开法求。先看怎么求。将转换为16进制数表示，改写成 = 0.n = 0. + 0.000n =   + Δx，其中 = 0. =  *  +  *  + * ，Δx = 0.000n =  *  +  *  + ...。所以 = =  。因为x0 <  < 0.69314718056 < 0.75 = 12/16，所以的取值范围是[0, 11]，的取值范围是[0, 15]，的取值范围是[0, 15]。根据的有限个不同取值将 、 和 分别预先算出做成表，计算时通过查表得到三个相应的值，再将这三个值相乘就得到的值了。再来看怎么求。0 < Δx  = 0.000n <  = 1/4096 << 1，用幂级数展开式求只要取前面4项即可保证精度了，所以用幂级数展开式求。

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

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)     求尾码。通过查表法求，通过幂级数展开式求，p = 即可求得，α = p – 1也可求得。尾码 =α*就得到了。

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

3，求自然对数

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

又 y = ，所以

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

即：

其中

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

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

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

1)     得到阶码E和尾码Mi

2)     根据阶码E求得b

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

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