其他分享
首页 > 其他分享> > 基于MCRA-OMLSA的语音降噪(三):实现(续)

基于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的根。

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

处理后得到

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

 

2,求自然指数

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

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

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

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

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

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

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

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

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

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

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

 

令p = 2D,从而 。令x00 = Dln2,有p = ex00。因为 0 < D < 1,又ln2 = 0.69314718056,所以 0 < x00 < 0.69314718056。此时若直接用ex00的幂级数展开式求p,计算时间还很长,若适当选取x0和Δx,使得Δx << 1,且 x00 = x0 + Δx,则有 p = ex0 + Δx = ex0eΔx。可分别求ex0和eΔx,然后再相乘就得到p。论文中用查表法求ex0,用幂级数展开法求eΔx。先看怎么求ex0。将x00转换为16进制数表示,改写成x00 = 0.q1q2q3q4q5n = 0.q1q2q3 + 0.000q4q5n =  x0 + Δx,其中x0 = 0.q1q2q3 = q1 * 16-1 + q2 * 16-2 + q3 * 16-3,Δx = 0.000q4q5n = q4 * 16-4 + q5 * 16-5 + ...。所以ex0 = eq1 * 16-1 + q2 * 16-2 + q3 * 16-3 = eq1 * 16-1eq2 * 16-2eq3 * 16-3。因为x0 < x00 < 0.69314718056 < 0.75 = 12/16,所以q1的取值范围是[0, 11],q2的取值范围是[0, 15],q3的取值范围是[0, 15]。根据qx的有限个不同取值将eq1 * 16-1 、eq2 * 16-2 和eq3 * 16-3 分别预先算出做成表,计算时通过查表得到三个相应的值,再将这三个值相乘就得到ex0的值了。再来看怎么求eΔx。0 < Δx  = 0.000q4q5n < 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)     求尾码。通过查表法求ex0,通过幂级数展开式求eΔx,p = ex0eΔx即可求得,α = p – 1也可求得。尾码 =α*252就得到了。

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

 

3,求自然对数

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

又 y = ex,所以

上式两边取自然对数,得到:

即:

其中

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相加就得到自然对数

标签:OMLSA,16,二进制,降噪,ax,尾码,x0,阶码,MCRA
来源: https://www.cnblogs.com/talkaudiodev/p/15765294.html