完成了g723源代码的分析,现作一些整理
1 信号高通滤波 Rem_Dc:
这个函数做高通滤波用的.将低频噪声滤除
滤波器的系统函数为
H(z)=(1-z^(-1)) / (1 - (127/128)*z^(-1))
将 单位圆上的值代入(即:cos(a) + sin(a)*i)
可以看出这个函数在π达到峰值,在 0 和 2π处是谷值,明显是高通滤波
迭代公式如下
y[n] = x[n] - x[n-1] + (127/128) * y[n-1]
为了避免出现除法运算,所有的值都扩大了
x[n] - x[n - 1] 即为itu代码中的fir部分
(127/128) * y[n-1] 即iir部分
我们来看看滤波代码:
/* Do the Fir and scale by 2 */
Acc0 = L_mult( Dpnt[i], (Word16) 0x4000 ) ; //lsc 乘 128 * 128 * 2
Acc0 = L_mac ( Acc0, CodStat.HpfZdl, (Word16) 0xc000 ) ; //lsc 乘-128 * 128 * 2 + acc0 0xc000是 -128 * 128的补码
CodStat.HpfZdl = Dpnt[i] ; //lsc 完成滤波器零点部分运算 acc0 = 256 * 128 * x[n] - 256 * 128 * x[n-1]
/* Do the Iir part */
Acc1 = L_mls( CodStat.HpfPdl, (Word16) 0x7f00 ) ; //lsc y[n-1] * (127 * 256 / 128 * 256) y[n - 1]在第一轮计算中已经和x[n]处于一个数量级了
Acc0 = L_add( Acc0, Acc1 ) ; //lsc 完成零极点和,完成了整个滤波运算, 显然,这里的运算结果扩大了 256 * 128倍
CodStat.HpfPdl = Acc0 ; //lsc 保存y[n] 至寄存器,形成下一步运算的y[n-1]
Dpnt[i] = round(Acc0) ; //lsc 取出运算结果的高16位
最终的结果是 输入信号缩小了一半
从另一个分支也可以看出
Dpnt[i] = shr( Dpnt[i], (Word16) 1 ) ;
如果不做高通滤波,就直接把信号缩小一半
补充说明
L_mult 计算两个相乘，并将结扩大2倍
L_mac(var1, var2, var3) 令后两个参数var2, var3相乘，并扩大两倍,再与第一个参数var1相加,
为了使加法运算在同一个数量级 var1必须预先扩大2倍
L_mls(var1,var2) 完成var1 var2的乘法,运算是这么进行的,先计算var1低16与 var2的乘法,再计算var1高16与var2乘法结果
对低16的计算结果做一些舍弃,最终的结果是缩小至真实乘法值的2^(-15)
round 取出一个数的高16位
2 计算lpc系数 Comp_Lpc
程序会保留之前帧的120个样点,与当前帧的240个样点,组合成一个新数组,lpc系统是针对这个数组进行计算的
lpc系数分四组运算,每组取180个样点,组与组之间存在着120个样点的交叠
ShAcf_sf[k] = Vec_Norm( Vect, (Word16) LpcFrame ) ;
这段代码对信号进行归一化处理,即找出绝对值最大的输入信号,将其归一化,即计算该输入信号移到"最大"值所需要
左移次数n,然后将180个样点进行左移n-3,完成了归一化,该函数同时返回 n-3,(实际上是保留了"3空位",可能是为了防上后继计算出现溢出)
可以看到itu使用了海明窗对信号进行截取,
窗函数的作用:截取信号,并尽量不影响信号的频域特性,所以窗函数的频域特性,应该是主瓣窄,并且旁瓣低,由于
同时满足两者有冲突,比如频域上是个冲激(这是最理想的,不会改变信号的任何频域特性),但在时域则为无穷个点,这种窗显然没有任何实际意义
于是数学家们设计了一些窗函数,分别适用于某种特定应用场合,海明窗是其中的一个
itu中的海明窗是放大了2^15的,对应的乘法则缩小2^(-15),所以加窗过程没有数值做任何放大
加窗代码如下
/* Apply the Hamming window */ //对信号加一个海明窗
for ( i = 0 ; i < LpcFrame ; i ++ )
Vect[i] = mult_r(Vect[i], HammingWindowTable[i]) ;
然后计算零时自相关,即能量,根据柯西定理,也可以知道,零时自相关是所有自相关里绝对值最大的
代码如下
/* Compute the zeroth-order coefficient (energy) */ //lsc 零时自相关,即能量
Acc1 = (Word32) 0 ;
for ( i = 0 ; i < LpcFrame ; i ++ ) {
Acc0 = L_mult( Vect[i], Vect[i] ) ;
Acc0 = L_shr( Acc0, (Word16) 1 ) ;
Acc1 = L_add( Acc1, Acc0 ) ;
}
可以看出,计算结果没有做任何放大,能量还会加1/1024噪声能量
然后将能量归一化
/* Normalize the energy */ //acc1,能量归一化,即左移至顶,需要多少次乘2
Exp = norm_l( Acc1 ) ;
Acc1 = L_shl( Acc1, Exp ) ;//能量归一化
取出归一化后能量的高16位
curAcf[0] = round( Acc1 ) ; //取出能量的高16位
如果能量的高16位为零(即能量为零),所有的自相关都认为是零
计算10个自相关系数,用于莱文森-德宾递推公式使用,代码片段如下:
for ( i = 1 ; i <= LpcOrder ; i ++ ) {
Acc1 = (Word32) 0 ;
for ( j = i ; j < LpcFrame ; j ++ ) {
Acc0 = L_mult( Vect[j], Vect[j-i] ) ;
Acc0 = L_shr( Acc0, (Word16) 1 ) ;
Acc1 = L_add( Acc1, Acc0 ) ;//lsc 计算自相关,没有任何放大
}
Acc0 = L_shl( Acc1, Exp ) ;/* lsc 已经归一化 */
Acc0 = L_mls( Acc0, BinomialWindowTable[i-1] ) ;/* 一个二项式窗函数 为何要加??? --- 加窗后的值应该没有做改变的 */
curAcf[i] = round(Acc0) ;//lsc 取出高16位
}
/* Save Acf scaling factor */
ShAcf_sf[k] = add(Exp, shl(ShAcf_sf[k], 1));
//lsc归一化的时候,本身移动了ShAcf_sf[k],而自相关为归一化后的信号相乘,意味着
//得到的结果本身已经扩了2^(ShAcf_sf[k]*2),归一化又移动了Exp,自然要更新相应的ShAcf_sf[k] = ShAcf_sf[k] * 2 + Exp
Durbin --- 莱文森-德宾递推算法,求出10个lpc系数
呼叫时传入的参数
Durbin( &UnqLpc[k*LpcOrder], &curAcf[1], curAcf[0], &Pk2 );
第一个参数不解释了,将记录Durbin函数计算得到的10个lpc系数
&curAcf[1] 为自相关系数组,用于莱文森-德宾递推,curAcf[0]为零时自相关,即能量,做为
递推时分母的第一个因子,Durbin函数的返回值,为残差信号的能量 &Pk2记录部分相关系数k(该参数似乎是在生成舒适噪音时才使用的)
莱文森-德宾递推公式证明:
首先从lpc系数的求解开始说吧
我们假定 s[n]是输入的语音信号
s`[n]是10阶预测信号
s`[n] = a10 * s[n - 10] + a9 * s[n - 9] + ... + a1 * s[n - 1]
我们取s`[n] 与 s[n] 方差最小值
(s[n] - s`[n])^2 ^2表示平方
然后对每个 a[i]求偏导，自然就得到了一个10元一次方程组,表示
10
Σa[i] * R(|k-i|) = R(k) k = 1,2,...,10 R(k)表示输入信号的自相关 --- 方程组1
i=1
即莱文森-德宾递推公式就是一种适合于计算机实现的解这个10元一次方程组的算法
下面提到内积和正交的概念
A(z) B(z)分别表示前向预测与后向预测的逆滤波器的系统函数
则它们关于输入s[n]的内积这么定义
10 11
Σ Σ a[i]*b[k] R(|i-k|) R(n)同样表示输入信号s[n]的自相关函数
i=1 k=1
记作 <A(z),B(z)>
如果内积为零，则被称之为正交
这里可以立即得出这么一个结论
<A(z), z^(-l)>一定为零,这个参见方程组1
也就是说A(z)与z^(-l) l=1,2,...10正交时,A(z)是该阶次下最优估计的逆滤波器
开始递推,首先从零开始,零阶时,啥都没有
最优的逆滤波器自然就是
A0(z)=1
B0(z)=z^-1
构造出递推公式
A[i](z) = A[i-1](z) + ki * B[i - 1](z)
B[i](z) = z^-1{ B[i - 1](z) + ki * A[i-1](z) } ----- 等式2
其中 ki = -{ <A[i-1](z),B[i - 1](z)> } / { B[i - 1](z), B[i - 1](z) } ---- 等式1
现在只需要证 A[i](z) 与 z^-i正交即可 (当然B[i](z)也要与z^-i,两个证明的过程差不多)
将A[i](z) = A[i-1](z) + ki * B[i - 1](z)代入 <A[i](z), z^-i>
我们立该得到
<A[i-i](z), z^-i> + ki<B[i - 1](z), z^-i> = 0;
求出满足这个条件的ki就是了
<A[i-i](z), z^-i> 实际上与 <A[i-i](z), B[i-1](z)>是相等的,为什么呢？
因为A[i-1](z)是最优估计，那它一定与 z^-l (l=1,2,...,i-i)正交,而 B[i-i](z)的最高阶系数是1...
同理<B[i - 1](z), z^-i> 与<B[i-i](z), B[i-1](z)>也是相等的...
自然ki就是等式1的那种形式,
莱文森-德宾递推公式至此证明完毕
ki递推公式化简
i-1
分子可以化简为 R[i] + Σ a(m-1)[n] * R(|l-i|) 这个可以用<A[i-i](z), z^-i> 直接推导出
l=1
分母可以化简为 (1-k(l)^2)*(1-k(l-1)^2) ... (1-k(1)^2)*R(0)
这是由于 <z^-1 * F(z), z^-1 * G(z)> = <F(z), G(z)> = <F(1/z), G(1/z)> 由内积的定理可以直接证出
用 <B[i](z), B[i](z)> 等式2进行代换,再进行相应的合并同类项处理,就可以得到化简后的分子
递推过程中,相应的更新每一轮的自相关系数:
a(i+1)(j) = ai(j) + k ai(m-j)
这个是从 A(i+1)(z) = Ai(z) + k Bi(z),推出的
我们可以证明 Bi(z) 的系数 与 Ai(z)的系数时"相反"的, 即把 Ai(z)的系数逆序排序其实就是Bi(z)的系数
Durbin函数代码如下:
Word16 Durbin( Word16 *Lpc, Word16 *Corr, Word16 Err, Word16 *Pk2 )
{
int i,j ;
Word16 Temp[LpcOrder] ;
Word16 Pk ;
Word32 Acc0,Acc1,Acc2 ;
/*
* Initialize the LPC vector
*/
for ( i = 0 ; i < LpcOrder ; i ++ )
Lpc[i] = (Word16) 0 ;
/*
* Do the recursion. At the ith step, the algorithm computes the
* (i+1)th - order MMSE linear prediction filter.
*/
for ( i = 0 ; i < LpcOrder ; i ++ ) {
/*
* Compute the partial correlation (parcor) coefficient
*/
/* Start parcor computation */
Acc0 = L_deposit_h( Corr[i] ) ;/* 将值扩至32位,对应递推中的R(m)值 */
Acc0 = L_shr( Acc0, (Word16) 2 ) ;/* 右移2位,缩小4倍,因为L_msu扩了2 ??? */
for ( j = 0 ; j < i ; j ++ )
Acc0 = L_msu( Acc0, Lpc[j], Corr[i-j-1] ) ;/* 即ai * R(|i-m|),acc0就是k的分子 */
Acc0 = L_shl( Acc0, (Word16) 2 ) ;
/* Save sign */
Acc1 = Acc0 ;
Acc0 = L_abs( Acc0 ) ;
/* Finish parcor computation */
Acc2 = L_deposit_h( Err ) ;
if ( Acc0 >= Acc2 ) {
*Pk2 = 32767;
break ;
}
Pk = div_l( Acc0, Err ) ;/* lsc k的分子除分母 得出k pk就是k */
if ( Acc1 >= 0 )
Pk = negate(Pk) ;//lsc 负负得正,正负得负
/*
* Sine detector
*/
if ( i == 1 ) *Pk2 = Pk;
/*
* Compute the ith LPC coefficient
*/
Acc0 = L_deposit_h( negate(Pk) ) ;
Acc0 = L_shr( Acc0, (Word16) 2 ) ;
Lpc[i] = round( Acc0 ) ;/* a[m]=k */
/*
* Update the prediction error
*/ //lsc 这段一比较绕,为了节省计算量itu使用了一些技巧,读者需要注意
Acc1 = L_mls( Acc1, Pk ) ;//lsc (-分子原值) * k 得出的是 -k^2 * Err k * Err = 分子
Acc1 = L_add( Acc1, Acc2 ) ;//lsc Err(Acc2) + ???
Err = round( Acc1 ) ;/* lsc 更新了k分母 */
/*
* Compute the remaining LPC coefficients
*/
for ( j = 0 ; j < i ; j ++ )
Temp[j] = Lpc[j] ;
for ( j = 0 ; j < i ; j ++ ) { /* lsc 更新lpc系数 anew(i) = ai + ka(m-i), 这个是由A(i+1)(z) = Ai(z) + kBi(z) */
Acc0 = L_deposit_h( Lpc[j] ) ;
Acc0 = L_mac( Acc0, Pk, Temp[i-j-1] ) ;
Lpc[j] = round( Acc0 ) ;
}
}
return Err ;//lsc 返回的是最后一个分母,即残差信号的能量,这个返回值被用于计算cng时的增益估计
}
ok 到这里已完成了10 lpc系数计算了,
由于一个lpc系数的误差,会影响信号的每个频段(这一点由A(z)的形式很容易想出)
毕竟 A(z) = P1 * Z^(-10) + P2 * z^(-9) + ....
这个系统函数不是以因式分解形式给出的,必然有些问题
我们必须把它转化为因式分解的形式,然后对其求根,而每个根的误差只会影响某频域的能量,这样
就适合进行矢量量化了
于是构造了lsf系数,它实际是A(z)的变形,并且能把所有的根者限定的单位圆上,这样求lsf系数就可以
简单地沿着单位圆进行暴力搜索.(解一元10次方程是很复杂的,所以itu将求根过程简化了)
至此,已经分析完了g723的整个lpc求解过程,
笔者将在后继章节中讲述 lpc转成lsf的过程,
完成了声道系数求解,剩下的就是激励编码了,这些笔者均会在后继章节中进行详细分析
未完待继...
版权木有,笔者不对本文转载产生的任何后果负责(如耳麦炸裂,显示器花屏等)
林绍川
2011-04-15 于杭州