g723源码详细分析(二) lpc转lsp

3 lpc系数转成lsf系数 AtoLsp

首先解释一下为什么需要lsf(lsp)系数

A(z)在形式上不是以因式分解给出的,如果直接设一个lpc码本表,

对lpc系数进行矢量量化,则某个lpc系数的误差对信号整个频域的影响,

为了消除这种误差影响,我们自然地想到了利用A(z)方程的因式分解形式,

也就是说对A(z)=0的根建立码本表,对其进行矢量量化,这样每个量化后根

的误差,只对某个频域的产生影响.

lsf系数即是通过对A(z)进行变形,然后求根,它是变形后的A(z)的根.

而且它的特点是,所有的根一定在单位圆上,这就对计算机处理提供了便利

(A(z)=0的根显然不会都在单位圆上,这将导致极复杂的求根过程以及量化计算量)

现在介绍A(z)的变形

我们定义

P(z) = A(z) + Z^-(P + 1)A(z^-1)

Q(z) = A(z) - Z^-(P + 1)A(z^-1) p是阶数,在lpc10阶里 p=10

而线谱对的真正含义是以下这两个形式的根即, P1(z)=0 以及Q1(z)=0的根

P1(z) = P(z)/(1 + z^(-1))

Q1(z) = Q(z)/(1 - z^(-1))

除(1 + z^(-1)) 是为了将 1 -1这样的根去除,降次

P(z) Q(z)对称多项式,特点是根一点在单位圆上,并且分别是在单位元上交错出现的

相应的证明请参考

F.K.Soong and B.H.Juang 的"Linear Spectrum Pair(LSP) and Speech Data Compression"

IEEE IC on ASSP,1984,pp.1.10.1-1.10.4

http://ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=1172448 这个为网址,文档似乎是要付费的...

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

维基百科上的解释,似乎跟P(z)是对称多项式的性质有关,导致了它的根出现在单位圆上

AtoLsp

lpc系数转成lsp系数

首先做了带宽扩展

/*

* Perform a bandwidth expansion on the LPC coefficients. This

* scales the poles of the LPC synthesis filter by a factor of

* 0.994.

*/ //lsc 理解 相当于对冲激响应序列依次乘以 0.994 0.994^n在傅里叶变换的幅值是一个碗状,在0-2派之间,

//与冲激响应进行卷积时,会让冲激响应的带宽扩展

for ( i = 0 ; i < LpcOrder ; i ++ )

LspVect[i] = mult_r( Lpc[i], BandExpTable[i] ) ;//lsc 相乘后得到的值是原值

P(z) 与 Q(z)的根是共轭对称的 很容易证出如果 z是p(z)=0的根,那个 z^-1也是p(z)的根,

在单位圆上,即说明p(z)的根是共轭对称出现

在进行长除法之前,先把 p(z) q(z)的定义详细说明一下,以及搜索时用到的一个变形也在下面给出:

P(z) = A(z) + Z^-(P + 1)A(z^-1) A(z) = 1 - a1*z^-1 - a2*z^-2 - ... - a10*z^-10 z^-n表示 z的-n次方

Q(z) = A(z) - Z^-(P + 1)A(z^-1)

P1(z) = P(z)/(1 + z^(-1))

Q1(z) = Q(z)/(1 - z^(-1)) ---- 根其实一样,做除法降次,去掉了 z=1 z=-1这样的根,因为这两个根本身不包含特征信息

P1(z)分母 = z^-11 + (-a1 - a10)z^-10 + (-a2 - a9)z^-9 + ... (-a1 - a10)z^-1 + 1

Q1(z)分母 = -z^-11 + (a1 - a10)z^-10 + (a2 - a9) z^-9 + ... (a1 - a10) z^-1 + 1

A(z) = 1/2 {P(z)(1 + z^{-1}) + Q(z)(1-z^{-1})} --- 从这里看出了itu对P1(z) Q1(z)的定义(网上流行有各种版本定义,让笔者纠结了很久...)

P1(z) = b0(z^-10 + 1) + b1(z^-9 + z^-1) + ... + b5 * z^-5

= z^-5(b(z^-5 + z^5) + b1(z^-4 + z^-4) + ... + b5) --- z= e^jw代入

= z^-5 * (b0 * 2cos(5w) + b1 * 2cos(4w) + ... + b5) b5是需要除2的,因为其它系数有个乘2因子

接下来,itu开始执行长除法,得到p1(z) q1(z)的系数

代码片段与注释

/*

* Set p_0 = q_0 = 1. The LPC coefficients are already scaled by

* 1/4. P(z) and Q(z) are scaled by an additional scaling factor of

* 1/16, for an overall factor of 1/64 = 0x02000000L.

*/

Lpq[0] = Lpq[1] = (Word32) 0x02000000L ;//lsc 0x80000000的1/64 1/4的lpc缩放源自语音信号本身缩小了1/2,而1/16的缩放见下面的代码段 这里认为是最高次 z^-11的系数

/*

* This loop computes the coefficients of P(z) and Q(z). The long

* division (to remove the real zeros) is done recursively.

*/ //lsc这个循环实际在执行多项式除法,求出pz qz的每阶的系数,求出前6个,后5是不需要求的，因为一样(可以通过多项式除法的定义，以及P(z)系数的对称性，推导出来)

for ( i = 0 ; i < LpcOrder/2 ; i ++ ) {

//lsc -11这个阶就不用求了,根据多项式除法定义,这个值一定是1

/* P(z) */

Acc0 = L_negate( Lpq[2*i+0] ) ; //lsc 根据多项式长除法,向低阶次"借"了 高阶次的系数(完成了长除法后的),应扣减

Acc1 = L_deposit_h( LspVect[i] ) ;

Acc1 = L_shr( Acc1, (Word16) 4 ) ;

Acc0 = L_sub( Acc0, Acc1 ) ;//lsc 系数的规律 如 -10这个阶,它应该是 -(a10 + a1),所以对应该阶次的系数还应再减去相应的值

Acc1 = L_deposit_h( LspVect[LpcOrder-1-i] ) ;

Acc1 = L_shr( Acc1, (Word16) 4 ) ;

Acc0 = L_sub( Acc0, Acc1 ) ;//lsc 这里完成了长除法 以-10这个阶为例 - a1 - a10 - 1(对照长除法结果, z^-11 + 1 - 1, 导到了-1的产生,而z^-10系数为-a1 - a10导致了这个结果)

Lpq[2*i+2] = Acc0 ;

/* Q(z) */ //lsc 这段同上,也是完成多项式长除法,就不详细分析了

Acc0 = Lpq[2*i+1] ;

Acc1 = L_deposit_h( LspVect[i] ) ;

Acc1 = L_shr( Acc1, (Word16) 4 ) ;

Acc0 = L_sub( Acc0, Acc1 ) ;

Acc1 = L_deposit_h( LspVect[LpcOrder-1-i] ) ;

Acc1 = L_shr( Acc1, (Word16) 4 ) ;

Acc0 = L_add( Acc0, Acc1 ) ;//lsc - a1 + a10 + 1(对照长除法结果 - z^-11 + (1 - 1)z^-10,产生-1项留给z^-10,而z^-10的系数应为a1 - a10,导致了 -a1 + a10 -1,取反,(1-z^-1)最终结果是 -a1 + a10 + 1)

Lpq[2*i+3] = Acc0 ;

最后一项没有乘2因子,所以要缩小一半

/*

* Divide p_5 and q_5 by 2 for proper weighting during polynomial

* evaluation.

*/

Lpq[LpcOrder+0] = L_shr( Lpq[LpcOrder+0], (Word16) 1 ) ;

Lpq[LpcOrder+1] = L_shr( Lpq[LpcOrder+1], (Word16) 1 ) ;//lsc 其它项都有乘2因子,唯最后一项是没有的

接下来对所有系数归一化,这里就不列举代码片段了,

然后是查cos表,在单位圆上暴力搜索过零值,求根

利用了这个等式求根

P1(z) = z^-5 * (b0 * 2cos(5w) + b1 * 2cos(4w) + ... + b5)

对cos表里的每个角度,计算当前值,与PrevVal值进行比较,如果异号,说明过零了

代码片段

CurrVal = L_mac( CurrVal, Spq[LpcOrder-2*j+k],

CosineTable[i*j%CosineTableSize] ) ;// j次方,对应复数就是角度被乘积,需要乘j

会根据PrevVal CurrVal 的比例,线性地计算小数点后7位,代码片段如下

Acc0 = L_abs( CurrVal ) ;

Acc1 = L_abs( PrevVal ) ;

Acc0 = L_add( Acc0, Acc1 ) ;

/* Normalize the sum */

Exp = norm_l( Acc0 ) ;

Acc0 = L_shl( Acc0, Exp ) ;

Acc1 = L_shl( Acc1, Exp ) ;

Acc1 = L_shr( Acc1, (Word16) 8 ) ;

LspVect[LspCnt] = div_l( Acc1, extract_h( Acc0 ) ) ;//此举的意义是保留小数点后7位,即得到的应该是 n.?? 0 < n < 256 n即当前下标i

最终得到的lsp是由9位整数部分,7位小数部分组成的,代码片段如下:

Exp = shl( (Word16) (i-1), (Word16) 7 ) ;

LspVect[LspCnt] = add( LspVect[LspCnt], Exp ) ;//得 9 + 7 即的浮点数 lsp频率

如果没有搜索到10根,则沿用之前帧的lsp系数.

致此,完成了lpc到lsp的转换

Lsp_Qnt

lsp量化,分将10个lsp系数分成 3 + 3 + 4 分别在码本里按欧式距离最小的规则进行矢量量化,待续

林绍川

2011.05.06 于杭州