实时音频编解码之十二Opus编码-SILK编码-基频估计
Posted shichaog
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了实时音频编解码之十二Opus编码-SILK编码-基频估计相关的知识，希望对你有一定的参考价值。
本文谢绝任何形式转载，谢谢。
4.1.9 silk_encode_frame_FLP
在调用该API函数之前，其原始数据会存放在&psEnc->state_Fxx[ 0 ].sCmn.inputBuf指向的地址中，4.1.6小节中silk_memcpy函数拷贝short类型输入，该函数调用若干核心函数完成声道参数和激励信号的编码压缩，调用层次结构清晰，后文若干小节依次分析核心函数。
//encode_frame_FLP.c
84 opus_int silk_encode_frame_FLP( silk_encoder_state_FLP *psEnc, // I/O 编码器状态
85 opus_int32 *pnBytesOut, // O 编码后负载字节数;
86 ec_enc *psRangeEnc, // I/O 区间编码状态
87 opus_int condCoding, // I 条件编码标志，对于40ms帧，该标志没有置位
88 opus_int maxBits, // I 如果大于零，则指示输出比特最大值
89 opus_int useCBR // I 强制使用定比特率标志
124 //当输入信号带宽发生变化时，使用低通滤波（ARMA方法）平滑带宽变化
//第一个参数是低通滤波器状态，第二个参数是输入信号，滤波后信号也存在该指针指向地址中，第三个参数/帧长，对于16kHz，20ms帧长，其长度为320
126 silk_LP_variable_cutoff( &psEnc->sCmn.sLP, psEnc->sCmn.inputBuf + 1, psEnc->sCmn.frame_length );
127
130 //将低通滤波后的信号拷贝到x_frame，
//x_frame是指向psEnc->x_buf + psEnc->sCmn.ltp_mem_length的指针,
//LA_SHAPE_MS是噪声形状分析的lookahead毫秒数，该宏为5，即向前看5ms，这会引入5ms的编码延迟
//这里将short类型数据转换为float类型数据
131 silk_short2float_array( x_frame + LA_SHAPE_MS * psEnc->sCmn.fs_kHz, psEnc->sCmn.inputBuf + 1, psEnc->sCmn.frame_length );
//基频搜索, 初始的LPC分析，层级关系见图4-3。
142 silk_find_pitch_lags_FLP( psEnc, &sEncCtrl, res_pitch, x_frame, psEnc->sCmn.arch );
//噪声形状分析
147 silk_noise_shape_analysis_FLP( psEnc, &sEncCtrl, res_pitch_frame, x_frame );
//LPC 和 LTP系数查找
152 silk_find_pred_coefs_FLP( psEnc, &sEncCtrl, res_pitch_frame, x_frame, condCoding );
//增益处理
157 silk_process_gains_FLP( psEnc, &sEncCtrl, condCoding );
//Low Bitrate Redundant编码
162 silk_LBRR_encode_FLP( psEnc, &sEncCtrl, x_frame, condCoding );
// 循环量化和熵编码以控制比特率
//噪声形状量化
196 silk_NSQ_wrapper_FLP( psEnc, &sEncCtrl, &psEnc->sCmn.indices, &psEnc->sCmn.sNSQ, psEnc->sCmn.pulses, x_frame );
//发声模型参数编码
205 silk_encode_indices( &psEnc->sCmn, psRangeEnc, psEnc->sCmn.nFramesEncoded, 0, condCoding );
//激励信号编码
210 silk_encode_pulses( psRangeEnc, psEnc->sCmn.indices.signalType, psEnc->sCmn.indices.quantOffsetType,
211 psEnc->sCmn.pulses, psEnc->sCmn.frame_length );
//增益量化
343 silk_gains_quant( psEnc->sCmn.indices.GainsIndices, pGains_Q16,
344 &psEnc->sShape.LastGainIndex, condCoding == CODE_CONDITIONALLY, psEnc->sCmn.nb_subfr );
//增益向量索引
347 gainsID = silk_gains_ID( psEnc->sCmn.indices.GainsIndices, psEnc->sCmn.nb_subfr );
4.1.10 pitch估计
在基于模型的发声方法中，pitch/基频是语音的重要特征，激励信号的周期性和pitch息息相关，pitch的准确性因而也是非常重要的，由于语音逐帧计算的，每帧语音都可以计算得到一个pitch，且帧间间隔非常短，在opus算法中，pitch是按照5ms一次计算的，因此pitch之间的相关性非常强， 从下图可以看出，短时间内pitch的变动通常不大，其连成的线被称为Pitch Contour（pitch包络），可以对包络进行压缩，opus采用的方法是矢量量化的方式，对于20ms编码包因其包含4个5ms子帧，因而其包络是4个值组成的向量，通过大量语音基频数据发现，这4个值的大小关系分布并不是完全随机的（正如下图显示的pitch关系，通常相差不大），通过聚类等算法，对于20ms 宽带语音帧，opus使用了34个类（34个长度为4的向量）描述20ms帧间的基频相对关系，这样在通过网络传输的时候，可以传输类（向量）索引值，通过一个索引值解码端可以从有34个码字的码本中索引基频包络对应的向量，这节约了网络传输带宽，另外这34个类出现的概率并不是一样的，有些概率大一些，有些概率小一些，通过纯数学的区间编码再一次编码以进一步压缩传输数据量，10ms语音帧和其类似，其包络是2个值组成的向量。Lag是以采样点计数的基音周期，Fs=8kHz时，Lag=160~20 (50Hz~400Hz)。
图：使用Aubio Pitch Detector绘制的pitch图
图中两段分别是笔者发中文5和中文6的时域波形，其下方数字指示的是pitch延迟，可见基频延迟具有相关性，比如以[130, 127, 122, 119, 117, 116, 114]这一段，如果以100为基准参考则，可以变为100 + [30, 27, 22, 19, 17, 16, 14]，显然在传输100再传输偏移量的方法占用的比特位更少（均小于32，可以用8比特传输偏移量），如果对中括号内的偏移进行矢量码本量化映射，则偏移可以用一个索引值替代，这进一步压缩了传输数据量，由于中括号内的偏移量相对集中，同等量化误差下，对偏移量使用矢量量化比直接对基频延迟进行矢量量化可以使得码本更小，搜索更为快捷。频谱图中基频延迟如下图所示，通常成年男性的基频比成年女性的基频更低，因而基频延迟值也会更小一些，此外，频谱图传递的另一个信息是基频和谐波是语音频谱的重要组成部分，因而正确计算其值对语音的质量有着重要影响。
图 男女声频谱图基频信息
Pitch lag的分析基于LPC滤波之后残差信号，如下从84行起的代码片段计算LPC系数并进行了滤波得到残差信号res。
//find_pitch_lags_FLP.c
36 void silk_find_pitch_lags_FLP(
37 silk_encoder_state_FLP *psEnc, // I/O Encoder state FLP
38 silk_encoder_control_FLP *psEncCtrl, // I/O Encoder control FLP
39 silk_float res[], // 输出 残差信号
40 const silk_float x[], // 输入语音信号 Speech signal
// 输入 运行时架构，pitch/LPC 计算较为耗时，实现了SIMD指令加速和纯C版本，根据arch值执行相应代码
41 int arch
42 )
43
// 根据LPC 阶数计算自相关向量 ，结果存放在auto_corr里，Wsig是加窗语音信号,
// 为了计算准确性，计算pitch的数据向前和向后各看了2ms，对于16kHz采样率而言，就是32个采样点，
// 因而Wsig存放信号长度为32（向前） + 320（当前帧） + 32（向后）组成，只对向前和向后采样点使用正弦窗
85 silk_autocorrelation_FLP( auto_corr, Wsig, psEnc->sCmn.pitch_LPC_win_length, psEnc->sCmn.pitchEstimationLPCOrder + 1 );
// 增加白噪声,作为能量的一部分
88 auto_corr[ 0 ] += auto_corr[ 0 ] * FIND_PITCH_WHITE_NOISE_FRACTION + 1;
// 用Schur算法计算反射系数，根据反射系数计算的LPC系数是稳定的，使用Schur算法的另一个原因是其递归求解特性会减少计算复杂度
// 由于命令行参数设置的复杂度是10，因而阶数为16
// 因而通过反射系数计算LPC系数，而非直接根据自相关计算LPC系数
91 res_nrg = silk_schur_FLP( refl_coef, auto_corr, psEnc->sCmn.pitchEstimationLPCOrder );
// 预测增益计算，该值越大，则值影响合成滤波器对带宽扩展形状的控制
94 psEncCtrl->predGain = auto_corr[ 0 ] / silk_max_float( res_nrg, 1.0f );
// 将反射系数转换为LPC系数，存放在数组A中
97 silk_k2a_FLP( A, refl_coef, psEnc->sCmn.pitchEstimationLPCOrder );
// 根据LPC系数，使用自回归滤波器对带宽扩展，其目的是增加LPC滤波器的稳定性，见1.4.8小节
100 silk_bwexpander_FLP( A, psEnc->sCmn.pitchEstimationLPCOrder, FIND_PITCH_BANDWIDTH_EXPANSION );
//使用LPC系数估计原始信号，原始信号和残差信号的差值存放于res变量中，A是100行计算得到的LPC系数。
105 silk_LPC_analysis_filter_FLP( res, A, x_buf, buf_len, psEnc->sCmn.pitchEstimationLPCOrder );
LPC系数计算基于自相关算法，得到自相关auto_corr之后计算反射系数refl_coef，然后根据反射系数递归得到LPC系数A（滤波器系数），然后用这组系数对原始信号进行滤波得到LPC方法估计的信号，后将真实信号x_buf和估计信号的残差（白化）存放在res变量用于计算pitch，白化的作用使得pitch分析对所有频谱的敏感性是一样的，这避免了谐波的干扰，这也提高了后面基于自相关判断语音/非语言的准确性，这里的LPC计算原理在第二章LPC一节已经涉及，这里不再进行原理分析。
在设置好pitch估计门限之后，调用silk_pitch_analysis_core_FLP进行pitch的核心估计计算。
// pitch 估计其门限设置，语音可能性越大，则门限越低，以防止漏掉语音
109 thrhld = 0.6f;
110 thrhld -= 0.004f * psEnc->sCmn.pitchEstimationLPCOrder;
111 thrhld -= 0.1f * psEnc->sCmn.speech_activity_Q8 * ( 1.0f / 256.0f );
112 thrhld -= 0.15f * (psEnc->sCmn.prevSignalType >> 1);
113 thrhld -= 0.1f * psEnc->sCmn.input_tilt_Q15 * ( 1.0f / 32768.0f );
pitch的范围是[2ms, 18ms），对应于（56Hz, 500Hz]，对应到采样点数在代码和文献中该值称为lag，pitch lag即是用采样点表示的pitch值，在pitch计算时opus先是从粗粒度到细粒度的过程，即先计算4kHz信号的pitch，然后根据4kHz信号计算的pitch范围估计8kHz信号pitch所处的范围，直接在该估计的范围内计算，这样做的目的是为了降低自相关计算的复杂度，这样省去了不必要的计算同时使得pitch估计更为准确，最后计算原始16kHz信号的pitch和其包络。下图中的输入信号1指的是LPC滤波后残差信号res（白化），高复杂度（complexity参数）时LPC的阶数是16，中等复杂度LPC阶数是12，低复杂度LPC阶数是8，阶数越高白化效果越好（相减之后去相关），下图中标号2是代码中stage 1计算得到的Lag候选值数组，3是代码中stage 2计算得到的Lag候选值数组，4是自相关门限，5是语音/非语言标志，6是pitch自相关，7是pitch lag结果。
图：pitch_estimator处理框图
pitch分析首先确定是语音还是非语音，只有判为语音的时候才会做pitch分析，对于语音情况，按帧长为5ms逐个计算lag值，每次会得到4个lag（20ms对应4个5ms），对白化之后的信号，时域自相关越大的情况则越可能是lag存在的时候，这是因为并不能完全意义上去相关，为了降低复杂度，使用三个stage计算得到pitch延迟，第一个阶段是对4kHz信号（图中2倍下采样）根据lag可能区间值（2ms-18ms），对于4kHz即8~72点，做自相关计算。两次下采样的代码如下：
//pitch_analysis_core_FLP.c
//由于命令行参数指定了编码采样率为16kHz，即 Fs_kHz=16，因而进入第一个if将信号下采样到8 kHz
67 opus_int silk_pitch_analysis_core_FLP( // O 语音估计: 0 语音, 1 非语音
68 const silk_float *frame, // I 信号长度 PE_FRAME_LENGTH_MS*Fs_kHz
69 opus_int *pitch_out, // O Pitch lag 值 [nb_subfr]
70 opus_int16 *lagIndex, // O Lag 索引
71 opus_int8 *contourIndex, // O Pitch 包络索引
72 silk_float *LTPCorr, // I/O 归一化互相关; 输入值: 计算自前一帧
73 opus_int prevLag, // I 前一帧延迟值; 非语音设置为0
74 const silk_float search_thres1, // I 第一阶段lag延迟值门限值 0 - 1
75 const silk_float search_thres2, // I 第二阶段lag延迟值门限值 0 - 1
76 const opus_int Fs_kHz, // I 采样频率 (kHz)
77 const opus_int complexity, // I 复杂度设置, 0-2, 2 复杂度最高
78 const opus_int nb_subfr, // I 5 ms 子帧的个数
79 int arch // I 运行时架构，优化用
80 )
81
136 if( Fs_kHz == 16 )
137 //重采样16 -> 8 kHz
138 opus_int16 frame_16_FIX[ 16 * PE_MAX_FRAME_LENGTH_MS ];
139 silk_float2short_array( frame_16_FIX, frame, frame_length );
140 silk_memset( filt_state, 0, 2 * sizeof( opus_int32 ) );
141 silk_resampler_down2( filt_state, frame_8_FIX, frame_16_FIX, frame_length );
142 silk_short2float_array( frame_8kHz, frame_8_FIX, frame_length_8kHz );
143 else if( Fs_kHz == 12 )
//silk_resampler_down2_3是2/3下采样信号，即从 12 -> 8 kHz
145 opus_int16 frame_12_FIX[ 12 * PE_MAX_FRAME_LENGTH_MS ];
146 silk_float2short_array( frame_12_FIX, frame, frame_length );
147 silk_memset( filt_state, 0, 6 * sizeof( opus_int32 ) );
148 silk_resampler_down2_3( filt_state, frame_8_FIX, frame_12_FIX, frame_length );
149 silk_short2float_array( frame_8kHz, frame_8_FIX, frame_length_8kHz );
150 else
151 celt_assert( Fs_kHz == 8 );
152 silk_float2short_array( frame_8_FIX, frame, frame_length_8kHz );
153
//再将8kHz信号下采样到 4 kHz，这里采用的是和重采样小节一样的方法
156 silk_memset( filt_state, 0, 2 * sizeof( opus_int32 ) );
157 silk_resampler_down2( filt_state, frame_4_FIX, frame_8_FIX, frame_length_8kHz );
158 silk_short2float_array( frame_4kHz, frame_4_FIX, frame_length_4kHz );
接下来根据4kHz信号估计pitch lag，
//pitch_analysis_core_FLP.c silk_pitch_analysis_core_FLP函数
//frame_4kHz是下采样之后4kHz信号的起始地址，sf_length_4kHz是4kHz 信号sub-frame=20，左移两位之后为80。
//则target_ptr指向了4kHz信号的中间一半，4kHz下采样之后的帧长为160.
169 target_ptr = &frame_4kHz[ silk_LSHIFT( sf_length_4kHz, 2 ) ];
170 for( k = 0; k < nb_subfr >> 1; k++ )
175 basis_ptr = target_ptr - min_lag_4kHz;
//计算max_lag_4kHz - min_lag_4kHz + 1点的自相关，计算自相关的用到的点数为sf_length_8kHz，结果存放于xcorr。
//xcorr(0)=r_xx(max_lag_4kHz)
//xcorr(1)=r_xx(max_lag_4kHz-1)
181 celt_pitch_xcorr( target_ptr, target_ptr-max_lag_4kHz, xcorr, sf_length_8kHz, max_lag_4kHz - min_lag_4kHz + 1, arch );
// xcorr最后一个值，是间隔为min_lag_4kHz两个序列的自相关，因而这里取了最后一个值，存放在了min_lag_4kHz位置，即归一化后存放于C[ 0 ][ min_lag_4kHz ]
184 cross_corr = xcorr[ max_lag_4kHz - min_lag_4kHz ];
185 normalizer = silk_energy_FLP( target_ptr, sf_length_8kHz ) +
186 silk_energy_FLP( basis_ptr, sf_length_8kHz ) +
187 sf_length_8kHz * 4000.0f;
189 C[ 0 ][ min_lag_4kHz ] += (silk_float)( 2 * cross_corr / normalizer );
// 以下是获取max_lag_4kHz-min_lag_4kHz个长度归一化自相关序列，归一化因子采用了递归迭代
192 for( d = min_lag_4kHz + 1; d <= max_lag_4kHz; d++ )
193 basis_ptr--;
199 cross_corr = xcorr[ max_lag_4kHz - d ];
// 增加新值，剔除旧值，更新归一化因子
202 normalizer +=
203 basis_ptr[ 0 ] * (double)basis_ptr[ 0 ] -
204 basis_ptr[ sf_length_8kHz ] * (double)basis_ptr[ sf_length_8kHz ];
205 C[ 0 ][ d ] += (silk_float)( 2 * cross_corr / normalizer );
206
// 将数据后移40个点，8kHz计数的子帧长度，再次计算归一化自相关，并和之前的累加，得到帧内累积自相关
208 target_ptr += sf_length_8kHz;
209
归一化采用了算术均值而非几何均值，即采用了 C o p u s = x T y 1 2 ( x T x + y T y ) C_opus=\\frac\\mathbfx^T\\mathbf y\\frac12(\\mathbfx^T\\mathbfx + \\mathbfy^T\\mathbfy) Copus=21(xTx+yTy)xTy，在得到归一化累积自相关之后，对其进行排序，根据门限得到候选pitch lag所在的位置，映射到8kHz之后再次进行估算。
//pitch_analysis_core_FLP.c silk_pitch_analysis_core_FLP函数
//根据复杂度确定search的长度，这里的情况是2。
//length_d_srch指示的是8kHz是pitch lag的候选数量，对于这里的情况是8，
//即4kHz前8个最大的归一化累积自相关索引作为8kHz pitch查找指示
217 length_d_srch = 4 + 2 * complexity;
//降序排列归一化累积自相关值，d_srch存放的是由大到小排序后该值对应的原始位置索引，因为是原地址排序（对C数组排序）。
219 silk_insertion_sort_decreasing_FLP( &C[ 0 ][ min_lag_4kHz ], d_srch, max_lag_4kHz - min_lag_4kHz + 1, length_d_srch );
得到排序之后根据归一化累积自相关判断4kHz情况下的pitch 所在位置的候选值，8kHz在4kHz指示的位置附近再进行细粒度的查找。
//pitch_analysis_core_FLP.c silk_pitch_analysis_core_FLP函数
// 如果最大的归一化累积自相关小于0.2，则认为不存在pitch，所有lag相关字段都清零
222 Cmax = C[ 0 ][ min_lag_4kHz ];
223 if( Cmax < 0.2f )
224 silk_memset( pitch_out, 0, nb_subfr * sizeof( opus_int ) );
225 *LTPCorr = 0.0f;
226 *lagIndex = 0;
227 *contourIndex = 0;
228 return 1;
229
//选择前8个归一化累积自相关最大值作为候选pitch位置
232 for( i = 0; i < length_d_srch; i++ )
// 将大于门限的归一化累积自相关所在的索引映射到8kHz情况下的索引
// 新的索引长度记录在length_d_srch变量
234 if( C[ 0 ][ min_lag_4kHz + i ] > threshold )
235 d_srch[ i ] = silk_LSHIFT( d_srch[ i ] + min_lag_4kHz, 1 );
236 else
237 length_d_srch = i;
238 break;
239
240
8kHz的处理方式和4kHz类似，但是不再全局意义上做自相关后查找pitch lag可能的值，而是对4kHz情况的pitch lag位置进行细粒度的滤波之后的可能位置估算互相关，8kHz情况对pitch lag做了个偏移以避免pitch doubling，计算的结果存放在C数值中，第一维是子帧的索引，第二维是自相关值存储位置索引。
最后的估算部分直接基于输入信号（8kHz，12kHz或者16kHz），对于命令行的测试情况是基于16kHz信号估算。如果前一帧存在pitch lag的话，则基于pitch lag的短时相关性在前一帧附近搜索（范围限制在55.6Hz~500Hz）若干pitch lag，对于估算的每个pitch lag，会遍历码本中对应的pitch 包络相对偏移以计算自相关累积值，取最大累积值所在的索引为pitch lag值，pitch包络基于5ms帧计算。比如表34的silk_CB_lags_stage2是20ms帧的pitch包络，其第一个维度是帧数，第二个维度是码本的数量，第一列表示四个子帧的pitch相对偏移为零，第二列表示四个子帧的平均pitch lag和第三个子帧相等，第一个子帧比平均值大两个索引值，其它的以此类推，这样可以计算出自相关累积值，对于8kHz计算出的候选pitch所在的索引，根据此码本计算自相关累积值，然后将最大自相关累积值出现的位置定位pitch lag，这样根据一个索引值和一个pitch lag绝对值就可以得到四个子帧（20ms帧长）pitch lag，如果绝对值pitch lag使用前一帧的，则可以用一个索引指示四个子帧的pitch lag。
//pitch_est_tables.c
//表34: Codebook Vectors for Subframe Pitch Contour:NB, 20 ms Frames 8kHz情况
53 const opus_int8 silk_CB_lags_stage2[ PE_MAX_NB_SUBFR ][ PE_NB_CBKS_STAGE2_EXT ] =
54
55 0, 2,-1,-1,-1, 0, 0, 1, 1, 0, 1,
56 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0,
57 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0,
58 0,-1, 2, 1, 0, 1, 1, 0, 0,-1,-以上是关于实时音频编解码之十二Opus编码-SILK编码-基频估计的主要内容，如果未能解决你的问题，请参考以下文章