简介
FLAC是一种开源的无损压缩编码方法，主要对音频进行压缩，支持对WAV、RAW、RIFF等等格式的音频文件进行压缩编码，以及压缩后解码为原始音频文件。FLAC的使用和Header信息等可在FLAC的官方网站上找到，而对于FLAC中具体的压缩算法的介绍难以找到。本文介绍FLAC的四种编码方式之一——LPC编码。
事实上，对于一段音频，如双声道的，FLAC首先进行分块称为block，分块后编码为frame，然后分块后单个声道的子块称为subblock，编码后成为subframe。对于每个subframe，编码方法可能是verbatim（完整的编码）、constant（如果是常值信号）、fixed或lpc编码四种之一，选择的标准是编码后比特数最少的。
之后可能会更新关于FLAC源代码的其它部分，如对残差部分使用的rice编码等。【PS公式在这用latex敲的，的确latex的比word原先自带的好看很多，word中公式好像也有latex选项】
调用关系
以flac的1.3.2的源代码为例，在FLAC中，lpc编码针对的是一个个的subframe，单声道的音频。因为是第一次看一个具体的开源项目的源代码，没有看代码的经验，找到编码的算法就花了不少时间。
- flac1.3.2\flac\src\main.c中的main()中调用do_it()，do_it()在line 521调用encode_file()，encode_file()在line 1974调用flac__encode_file()。
- flac\src\encode.c中, flac__encode_file()中调用EncoderSession_process()，也就是FLAC__stream_encoder_process()进行编码，FLAC__stream_encoder_process()定义在encode.c中。
- FLAC__stream_encoder_process()调用process_frame_()对每个block处理，process_frame_()调用process_subframes_()来对每个subblock进行具体处理。
- 编码算法体现在process_subframe_()中，这个函数在flac-1.3.2\src\libFLAC\stream_encoder.c中，是重点关注的对象。
- 与lpc算法密切相关的函数有FLAC__lpc_window_data()【加窗】，encoder->private_->local_lpc_compute_autocorrelation()【用函数指针实现的类似C++中的类成员函数，用于计算自相关系数】，FLAC__lpc_compute_lp_coefficients()【计算lpc系数，lpc编码的核心部分】，FLAC__lpc_compute_best_order()【通过比特数判断需要用几阶的lpc，里面涉及一个参数，即精度参数，也是移位参数，称为shift】，FLAC__lpc_compute_expected_bits_per_residual_sample()【编码残差的每个采样点的比特数】。
算法说明
本文主要记录的是encoder->private_->local_lpc_compute_autocorrelation()和FLAC__lpc_compute_lp_coefficients()函数，第一个用于求解自相关系数，第二个用于计算lpc系数。
FLAC中求解lpc系数采用的是自相关法，所以需要求解自相关系数。
假设一段信号s(n)，可以用过去N个采样点的信号来预测当前采样点的信号
s
^
(
n
)
=
−
∑
k
=
1
N
a
N
(
k
)
s
(
n
−
k
)
\hat{s}(n) =-\sum_{k=1}^{N}{a_N(k)s(n-k)}
s^(n)=−k=1∑NaN(k)s(n−k)
其中，
a
N
(
k
)
a_N(k)
aN(k)为lpc采用N阶时的第k个预测系数。预测的误差为
e
(
n
)
=
s
(
n
)
−
s
^
(
n
)
e(n)=s(n)-\hat{s}(n)
e(n)=s(n)−s^(n)。使用最小均方误差准则来求解出每个系数
a
N
(
k
)
a_N(k)
aN(k)，使得总误差最小，总误差
E
N
E_N
EN表示为
E
N
=
∑
n
=
−
∞
+
∞
e
2
(
n
)
=
∑
n
=
−
∞
+
∞
[
s
(
n
)
+
∑
k
=
1
N
a
N
(
k
)
s
(
n
−
k
)
]
2
E_N=\sum_{n=-\infty}^{+\infty}e^2(n)=\sum_{n=-\infty}^{+\infty}[s(n)+ \sum_{k=1}^{N}a_N(k)s(n-k)]^2
EN=n=−∞∑+∞e2(n)=n=−∞∑+∞[s(n)+k=1∑NaN(k)s(n−k)]2
要让总误差最小，总误差最小时，满足
∂
E
n
∂
α
k
=
0
,
1
≤
k
≤
N
\frac{\partial E_n}{\partial \alpha_k} = 0, 1 \leq k \leq N
∂αk∂En=0,1≤k≤N
s
(
n
)
s(n)
s(n)的自相关函数
R
(
k
)
R(k)
R(k)为
R
(
k
)
=
∑
n
=
−
∞
∞
s
(
n
)
s
(
n
−
k
)
,
0
≤
k
≤
N
R(k)=\sum_{n=-\infty}^{\infty}s(n)s(n-k),0\leq k \leq N
R(k)=n=−∞∑∞s(n)s(n−k),0≤k≤N
R
(
k
)
R(k)
R(k)为偶函数，且
E
N
E_N
EN可以用R(k)来表示，表示为
E
n
=
R
(
0
)
+
∑
k
=
1
N
a
N
(
k
)
R
(
k
)
E_n = R(0) + \sum_{k=1}^{N}a_N(k)R(k)
En=R(0)+k=1∑NaN(k)R(k)
偏导数为0，导致了如下方程组的求解，该方程组称为Yule-Walker方程，自相关系数矩阵为
N
×
N
N\times N
N×N阶矩阵，称为Toeplitz矩阵：
R
N
⋅
a
N
=
−
r
N
R_N \cdot a_N=-r_N
RN⋅aN=−rN
其中
a
N
=
[
a
N
(
1
)
,
a
N
(
2
)
,
.
.
.
,
a
N
(
N
)
]
T
a^N=[a_N(1),a_N(2),...,a_N(N)]^T
aN=[aN(1),aN(2),...,aN(N)]T
r
N
=
[
R
(
1
)
,
R
(
2
)
,
.
.
.
,
R
(
N
)
]
T
r^N=[R(1),R(2),...,R(N)]^T
rN=[R(1),R(2),...,R(N)]T
对于该方程，利用Toeplitz矩阵的性质，有高效的Levinson-Durbin递推算法。算法如下：
E
0
=
R
(
0
)
E_0 = R(0)
E0=R(0)
ρ
m
=
a
m
(
m
)
=
−
R
(
m
)
+
r
m
−
1
b
t
a
m
−
1
E
m
−
1
=
−
R
(
m
)
+
r
m
−
1
b
t
a
m
−
1
R
(
0
)
+
r
m
−
1
t
a
m
−
1
\rho_m = a_m(m) = -\frac{R(m)+r^{bt}_{m-1}a_{m-1}}{E_{m-1}} = -\frac{R(m)+r^{bt}_{m-1}a_{m-1}}{R(0)+r^{t}_{m-1}a_{m-1}}
ρm=am(m)=−Em−1R(m)+rm−1btam−1=−R(0)+rm−1tam−1R(m)+rm−1btam−1
a
m
(
k
)
=
a
m
−
1
(
k
)
+
ρ
m
a
m
−
1
(
m
−
k
)
,
k
=
1
,
2
,
.
.
.
,
m
−
1
a_m(k) = a_{m-1}(k)+\rho_ma_{m-1}(m-k) ,\qquad k = 1,2,...,m-1
am(k)=am−1(k)+ρmam−1(m−k),k=1,2,...,m−1
E
m
=
E
m
−
1
[
1
−
a
m
2
(
m
)
]
=
E
m
−
1
[
1
−
ρ
m
2
]
E_m = E_{m-1}[1-a_m^2(m)] = E_{m-1}[1-\rho_m^2]
Em=Em−1[1−am2(m)]=Em−1[1−ρm2]
最后求得均方差
E
N
=
R
(
0
)
∏
i
=
1
p
(
1
−
k
i
2
)
E_N =R(0)\prod_{i=1}^p(1-k_i^2)
EN=R(0)∏i=1p(1−ki2)，和每个系数
a
m
(
k
)
a_m(k)
am(k)。
代码说明
- 求自相关的代码，这一段是弃用了的，但是利于理解。data是输入信号， data_len是长度，lag是滞后阶数，autoc是最后结果。
void FLAC__lpc_compute_autocorrelation(const FLAC__real data[], unsigned data_len, unsigned lag, FLAC__real autoc[])
{
/* a readable, but slower, version */
#if 0
FLAC__real d;
unsigned i;
FLAC__ASSERT(lag > 0);
FLAC__ASSERT(lag <= data_len);
/*
* Technically we should subtract the mean first like so:
* for(i = 0; i < data_len; i++)
* data[i] -= mean;
* but it appears not to make enough of a difference to matter, and
* most signals are already closely centered around zero
*/
while(lag--) {
for(i = lag, d = 0.0; i < data_len; i++)
d += data[i] * data[i - lag];
autoc[lag] = d;
}
- 计算lpc系数。autoc是上一个函数计算的自相关系数，max_order是lpc的阶数，由另一个函数传递进来，lp_coeff是lpc系数，是二维数组，阶数为9也就是max_order=9的lpc系数最终存放在lp_coeff[8][0…8]中，error[j]是对应的阶数为j+1阶的均方误差。
void FLAC__lpc_compute_lp_coefficients(const FLAC__real autoc[], unsigned *max_order, FLAC__real lp_coeff[][FLAC__MAX_LPC_ORDER], double error[])
{
unsigned i, j;
double r, err, lpc[FLAC__MAX_LPC_ORDER];
FLAC__ASSERT(0 != max_order);
FLAC__ASSERT(0 < *max_order);
FLAC__ASSERT(*max_order <= FLAC__MAX_LPC_ORDER);
FLAC__ASSERT(autoc[0] != 0.0);
err = autoc[0];
for(i = 0; i < *max_order; i++) {
/* Sum up this iteration's reflection coefficient. */
r = -autoc[i+1];/*注①*/
for(j = 0; j < i; j++)
r -= lpc[j] * autoc[i-j]; /*注②*/
r /= err; /*注③*/
/* Update LPC coefficients and total error. */
lpc[i]=r;
for(j = 0; j < (i>>1); j++) { /*注④*/
double tmp = lpc[j];
lpc[j] += r * lpc[i-1-j]; /*注⑤*/
lpc[i-1-j] += r * tmp;
}
if(i & 1) /*注⑥*/
lpc[j] += lpc[j] * r;
err *= (1.0 - r * r); /*注⑦*/
/* save this order */
for(j = 0; j <= i; j++)
lp_coeff[i][j] = (FLAC__real)(-lpc[j]); /* negate FIR filter coeff to get predictor coeff */
error[i] = err;
/* see SF bug https://sourceforge.net/p/flac/bugs/234/ */
if(err == 0.0) {
*max_order = i+1;
return;
}
}
}
注释①：
r
=
−
R
(
m
)
r=-R(m)
r=−R(m)，注意这里有负号，和之前公式中是一样的
注释②：通过循环计算
r
m
−
1
b
t
⋅
a
m
−
1
r_{m-1}^{bt}\cdot a_{m-1}
rm−1bt⋅am−1
注释③：计算
a
m
(
m
)
=
−
R
(
m
)
+
r
m
−
1
b
t
a
m
−
1
R
(
0
)
+
r
m
−
1
t
a
m
−
1
a_m(m)= -\frac{R(m)+r^{bt}_{m-1}a_{m-1}}{R(0)+r^{t}_{m-1}a_{m-1}}
am(m)=−R(0)+rm−1tam−1R(m)+rm−1btam−1
注释④：可以这么写，本质上是j<i，如此写，每次计算更新了i阶下的两个系数，分别为lpc[j]和lpc[i-1-j]。
注释⑤：计算
a
m
(
k
)
=
a
m
−
1
(
k
)
+
ρ
m
a
m
−
1
(
m
−
k
)
a_m(k) = a_{m-1}(k)+\rho_ma_{m-1}(m-k)
am(k)=am−1(k)+ρmam−1(m−k)
注释⑥：在④中，如果i是偶数，正好做完，因为第一个系数更新是lpc[i]=r；如果i是奇数，正中的系数尚未更新，需要单独更新。
注释⑦：当前阶数下的均方误差。在计算完系数后，之后在find_best_order中会用到误差。
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
