参量编码LPC:原理分析与仿真实践
在早期通信系统中,带宽资源有限,而波形编码要精确重现语音波形,这就需要较高的码率来传输大量数据,这在带宽不足的情况下就成了阻碍语音传输的大难题。随着通信技术不断进步,人们迫切希望能在低带宽环境下实现语音传输。
研究人员经过深入研究发现,语音信号有着独特的结构和统计特性。简单来说,语音信号并非杂乱无章,它内部存在一定规律。比如,在一段时间内,语音的频率变化、幅度变化等都有迹可循。这就意味着,我们并不需要像波形编码那样,把语音信号的每一个细节波形都完整传输过去,而是可以从语音信号里提取出一些关键的、能够代表其特征的 “参数”(这就是参量)。
1. 编码原理
语音信号具有一定的相关性,即当前样本与过去的样本存在某种线性关系。LPC 利用这一特性,通过线性预测分析来提取语音信号的参量。假设语音信号 s ( n ) s(n) s(n) 可以由其过去的 p p p 个样本的线性组合来预测:
s ^ ( n ) = ∑ k = 1 p a k s ( n − k ) \hat{s}(n)=\sum_{k = 1}^{p}a_k s(n - k) s^(n)=k=1∑paks(n−k)
其中, s ^ ( n ) \hat{s}(n) s^(n) 是预测值, a k a_k ak 是预测系数, p p p 是预测阶数。预测误差 e ( n ) e(n) e(n) 为:
e ( n ) = s ( n ) − s ^ ( n ) = s ( n ) − ∑ k = 1 p a k s ( n − k ) e(n)=s(n)-\hat{s}(n)=s(n)-\sum_{k = 1}^{p}a_k s(n - k) e(n)=s(n)−s^(n)=s(n)−k=1∑paks(n−k)
预测误差的均方值 ( E [ e 2 ( n ) ] ) (E[e^2(n)]) (E[e2(n)])就是对所有样本的预测误差的平方求平均,数学表达式为:
E [ e 2 ( n ) ] = 1 N ∑ n = 1 N e 2 ( n ) = 1 N ∑ n = 1 N [ s ( n ) − ∑ i = 1 p a i s ( n − i ) ] 2 E[e^2(n)]=\frac{1}{N}\sum_{n = 1}^{N}e^2(n)=\frac{1}{N}\sum_{n = 1}^{N}\left[s(n)-\sum_{i = 1}^{p}a_{i}s(n - i)\right]^2 E[e2(n)]=N1n=1∑Ne2(n)=N1n=1∑N[s(n)−i=1∑pais(n−i)]2
通过最小化 E [ e 2 ( n ) ] E[e^2(n)] E[e2(n)],可以找到一组最佳的预测系数 a i {a_{i}} ai,使得预测值 s ^ ( n ) \hat{s}(n) s^(n)尽可能接近真实值 s ( n ) s(n) s(n),从而实现对语音信号的有效预测和编码。 这需要用到求偏导数并令其为0的方法。
对 E [ e 2 ( n ) ] E[e^2(n)] E[e2(n)]关于 a j a_{j} aj ( j j j= 1, 2, ……, p)求偏导数:
-
根据复合函数求导法则以及求和运算的求导性质,我们有:
∂ E [ e 2 ( n ) ] ∂ a j = ∂ ∂ a j ( 1 N ∑ n = 1 N [ s ( n ) − ∑ i = 1 p a i s ( n − i ) ] 2 ) = 1 N ∑ n = 1 N ∂ ∂ a j [ s ( n ) − ∑ i = 1 p a i s ( n − i ) ] 2 \begin{align*} \frac{\partial E[e^2(n)]}{\partial a_{j}}&=\frac{\partial}{\partial a_{j}}\left(\frac{1}{N}\sum_{n = 1}^{N}\left[s(n)-\sum_{i = 1}^{p}a_{i}s(n - i)\right]^2\right)\\ &=\frac{1}{N}\sum_{n = 1}^{N}\frac{\partial}{\partial a_{j}}\left[s(n)-\sum_{i = 1}^{p}a_{i}s(n - i)\right]^2 \end{align*}\ ∂aj∂E[e2(n)]=∂aj∂ N1n=1∑N[s(n)−i=1∑pais(n−i)]2 =N1n=1∑N∂aj∂[s(n)−i=1∑pais(n−i)]2 -
令 u = s ( n ) − ∑ i = 1 p a i s ( n − i ) u = s(n)-\sum_{i = 1}^{p}a_{i}s(n - i) u=s(n)−∑i=1pais(n−i),根据复合函数求导公式 ∂ u 2 ∂ a j = 2 u ∂ u ∂ a j \frac{\partial u^2}{\partial a_{j}} = 2u\frac{\partial u}{\partial a_{j}} ∂aj∂u2=2u∂aj∂u,则:
∂ E [ e 2 ( n ) ] ∂ a j = 1 N ∑ n = 1 N 2 [ s ( n ) − ∑ i = 1 p a i s ( n − i ) ] ∂ ∂ a j [ s ( n ) − ∑ i = 1 p a i s ( n − i ) ] = 2 N ∑ n = 1 N [ s ( n ) − ∑ i = 1 p a i s ( n − i ) ] ( − s ( n − j ) ) = − 2 N ∑ n = 1 N s ( n − j ) [ s ( n ) − ∑ i = 1 p a i s ( n − i ) ] \begin{align*} \frac{\partial E[e^2(n)]}{\partial a_{j}}&=\frac{1}{N}\sum_{n = 1}^{N}2\left[s(n)-\sum_{i = 1}^{p}a_{i}s(n - i)\right]\frac{\partial}{\partial a_{j}}\left[s(n)-\sum_{i = 1}^{p}a_{i}s(n - i)\right]\\ &=\frac{2}{N}\sum_{n = 1}^{N}\left[s(n)-\sum_{i = 1}^{p}a_{i}s(n - i)\right]\left(-s(n - j)\right)\\ &=-\frac{2}{N}\sum_{n = 1}^{N}s(n - j)\left[s(n)-\sum_{i = 1}^{p}a_{i}s(n - i)\right] \end{align*} ∂aj∂E[e2(n)]=N1n=1∑N2[s(n)−i=1∑pais(n−i)]∂aj∂[s(n)−i=1∑pais(n−i)]=N2n=1∑N[s(n)−i=1∑pais(n−i)](−s(n−j))=−N2n=1∑Ns(n−j)[s(n)−i=1∑pais(n−i)] -
令 ∂ E [ e 2 ( n ) ] ∂ a j = 0 \frac{\partial E[e^2(n)]}{\partial a_{j}} = 0 ∂aj∂E[e2(n)]=0:
− 2 N ∑ n = 1 N s ( n − j ) [ s ( n ) − ∑ i = 1 p a i s ( n − i ) ] = 0 ∑ n = 1 N s ( n − j ) [ s ( n ) − ∑ i = 1 p a i s ( n − i ) ] = 0 ∑ n = 1 N s ( n − j ) s ( n ) = ∑ n = 1 N s ( n − j ) ∑ i = 1 p a i s ( n − i ) ∑ n = 1 N s ( n − j ) s ( n ) = ∑ i = 1 p a i ∑ n = 1 N s ( n − j ) s ( n − i ) \begin{align*} -\frac{2}{N}\sum_{n = 1}^{N}s(n - j)\left[s(n)-\sum_{i = 1}^{p}a_{i}s(n - i)\right]&=0\\ \sum_{n = 1}^{N}s(n - j)\left[s(n)-\sum_{i = 1}^{p}a_{i}s(n - i)\right]&=0\\ \sum_{n = 1}^{N}s(n - j)s(n)&=\sum_{n = 1}^{N}s(n - j)\sum_{i = 1}^{p}a_{i}s(n - i)\\ \sum_{n = 1}^{N}s(n - j)s(n)&=\sum_{i = 1}^{p}a_{i}\sum_{n = 1}^{N}s(n - j)s(n - i) \end{align*} −N2n=1∑Ns(n−j)[s(n)−i=1∑pais(n−i)]n=1∑Ns(n−j)[s(n)−i=1∑pais(n−i)]n=1∑Ns(n−j)s(n)n=1∑Ns(n−j)s(n)=0=0=n=1∑Ns(n−j)i=1∑pais(n−i)=i=1∑pain=1∑Ns(n−j)s(n−i) -
定义自相关函数:
令 r ( i , j ) = 1 N ∑ n = 1 N s ( n − i ) s ( n − j ) r(i, j)=\frac{1}{N}\sum_{n = 1}^{N}s(n - i)s(n - j) r(i,j)=N1∑n=1Ns(n−i)s(n−j),则上式可化为:
r ( 0 , j ) = ∑ i = 1 p a i r ( i , j ) r(0, j)=\sum_{i = 1}^{p}a_{i}r(i, j) r(0,j)=i=1∑pair(i,j) ( j = 1 , 2 , … , p ) (j = 1, 2, …, p) (j=1,2,…,p)
这是一个关于 a i a_{i} ai(i = 1, 2, …, p)的线性方程组,通过求解这个方程组,就可以得到使得预测误差均方值最小的预测系数 a i a_{i} ai(i = 1, 2, …, p)。一般可以将其写成矩阵形式 R a = b \mathbf{R}\mathbf{a}=\mathbf{b} Ra=b,其中 R \mathbf{R} R是由自相关函数 r ( i , j ) r(i, j) r(i,j)组成的矩阵, a = [ a 1 , a 2 , ⋯ , a p ] T \mathbf{a} = [a_{1}, a_{2}, \cdots, a_{p}]^T a=[a1,a2,⋯,ap]T是预测系数向量, b = [ r ( 0 , 1 ) , r ( 0 , 2 ) , ⋯ , r ( 0 , p ) ] T \mathbf{b} = [r(0, 1), r(0, 2), \cdots, r(0, p)]^T b=[r(0,1),r(0,2),⋯,r(0,p)]T,然后通过矩阵求逆等方法求解 a \mathbf{a} a。可以得到预测系数 a k a_k ak 。 -
此外,还需要确定增益 G G G ,通常定义为:
G = ∑ n = 1 N s 2 ( n ) ∑ n = 1 N s ^ 2 ( n ) G=\sqrt{\frac{\sum_{n = 1}^{N}s^2(n)}{\sum_{n = 1}^{N}\hat{s}^2(n)}} G=∑n=1Ns^2(n)∑n=1Ns2(n)
从公式可以看出,G是原始信号 s ( n ) s(n) s(n)的能量( ∑ n = 1 N s 2 ( n ) \sum_{n = 1}^{N}s^2(n) ∑n=1Ns2(n))与预测信号 s ^ ( n ) \hat{s}(n) s^(n)的能量( ∑ n = 1 N s ^ 2 ( n ) \sum_{n = 1}^{N}\hat{s}^2(n) ∑n=1Ns^2(n))之比的平方根。它反映了原始信号能量相对于预测信号能量的缩放程度。如果 G > 1 G > 1 G>1,说明原始信号的能量大于预测信号的能量;如果 G = 1 G = 1 G=1,表示原始信号和预测信号的能量相等;若 G < 1 G < 1 G<1,则原始信号能量小于预测信号能量。
2 LPC编码算法Matlab仿真
使用线性预测编码(LPC)方法对模拟话音信号进行处理,计算预测信号,并将原始信号与预测信号进行可视化对比,同时分析预测误差。
2.1 信号生成
- 参数设定:首先设定了模拟话音信号的相关参数,包括采样频率
fs
(8000 Hz)和信号持续时间duration
(1 秒)。采样频率决定了信号在时间上的离散程度,而持续时间则确定了生成信号的长度。 - 时间向量生成:根据采样频率和持续时间,生成时间向量
t
,用于后续信号的时间索引。 - 多频率正弦信号叠加:为了模拟话音信号的复杂特性,使用多个不同频率(300 Hz、800 Hz 和 1200 Hz)的正弦信号进行叠加。通过循环遍历每个频率,将对应频率的正弦信号累加到初始为零的信号
s
上,得到模拟的话音信号。
2.2 信号分帧与处理
- 帧长确定:设定帧长为 20 毫秒,在 8000 Hz 的采样频率下,对应的样本数为
frame_length = 0.02 * fs
,即 160 个样本。这是因为在语音处理中,通常将信号分割成较短的帧进行处理,以保证信号在局部范围内的平稳性。 - 提取话音帧:从模拟信号
s
中提取第一帧信号first_frame
,后续的 LPC 处理将基于这一帧信号进行。
2.3 线性预测编码(LPC)
- 预测阶数设定:设定预测阶数
p
为 12,表示使用过去的 12 个样本值来预测当前样本值。预测阶数的选择会影响预测的准确性和计算复杂度。 - 计算预测系数和增益:使用 MATLAB 的
lpc
函数对提取的话音帧first_frame
进行处理,计算得到 LPC 预测系数a
和增益g
。lpc
函数通过最小化预测误差的均方值来确定最佳的预测系数。 - 计算预测信号:利用
filter
函数根据预测系数a
计算预测信号predicted_frame
。filter
函数实现了一个数字滤波器,根据给定的滤波器系数对输入信号进行滤波处理,从而得到预测信号。
2.4 误差计算
- 计算误差信号:通过将原始信号
first_frame
减去预测信号predicted_frame
,得到误差信号error_signal
。误差信号反映了预测的准确性,其值越小表示预测越接近原始信号。
2.5 可视化
- 原始信号与预测信号对比图:创建一个图形窗口,使用
hold on
命令允许在同一图形上绘制多个对象。分别用蓝色实线和蓝色实心圆点绘制原始信号,用红色实线和红色实心圆点绘制预测信号。添加标题、坐标轴标签、网格线和图例,方便直观地比较原始信号和预测信号的差异。 - 误差信号图:创建另一个图形窗口,绘制误差信号
error_signal
,并添加标题、坐标轴标签和网格线,以便观察误差信号的变化情况。
2.6 误差分析
- 均方误差(MSE)计算:计算误差信号的均方误差(MSE),即误差信号平方的平均值。MSE 是衡量预测误差大小的常用指标,其值越小表示预测效果越好。
- 信噪比(SNR)计算:计算信噪比(SNR),通过比较原始信号的能量和误差信号的能量得到。SNR 表示信号中有用信号与噪声(误差)的比例,其值越大表示信号质量越好。最后将 MSE 和 SNR 的结果输出,方便对预测效果进行量化评估。
综上所述,代码通过模拟话音信号、进行 LPC 处理、计算预测信号和误差信号、可视化结果以及进行误差分析。
2.7 代码分析
% 采样频率
fs = 8000; % Hz
% 信号持续时间
duration = 1; % 秒
% 生成时间向量
t = 0:1/fs:duration - 1/fs;
% 生成多频率正弦信号叠加来模拟话音
frequencies = [300, 800, 1200]; % 正弦波频率
s = zeros(size(t));
for freq = frequenciess = s + sin(2*pi*freq*t);
end% 帧长,20 毫秒的帧,在 8000 Hz 采样率下是 160 个样本
frame_length = 0.02 * fs;
% 提取一个话音帧
first_frame = s(1:frame_length);% 预测阶数
p = 12;
% 计算 LPC 预测系数和增益
[a, g] = lpc(first_frame, p);% 计算预测信号
predicted_frame = filter([0 -a(2:end)], 1, first_frame);% 计算误差信号
error_signal = first_frame - predicted_frame;% 绘制原始信号和预测信号在同一张图中
figure;
hold on;
plot(t(1:frame_length), first_frame, 'b-', 'DisplayName', '原始信号');
plot(t(1:frame_length), first_frame, 'bo', 'MarkerFaceColor', 'b', 'MarkerSize', 3);
plot(t(1:frame_length), predicted_frame, 'r-', 'DisplayName', '预测信号');
plot(t(1:frame_length), predicted_frame, 'ro', 'MarkerFaceColor', 'r', 'MarkerSize', 3);
title('原始信号与预测信号对比');
xlabel('时间 (秒)');
ylabel('幅度');
grid on;
legend;% 绘制误差信号
figure;
plot(t(1:frame_length), error_signal);
title('误差信号');
xlabel('时间 (秒)');
ylabel('幅度');
grid on;% 分析误差,计算均方误差(MSE)
mse = mean(error_signal.^2);
fprintf('均方误差 (MSE): %.6f\n', mse);% 计算信噪比(SNR)
snr_value = 10 * log10(sum(first_frame.^2) / sum(error_signal.^2));
fprintf('信噪比 (SNR): %.2f dB\n', snr_value);
2.8 运行结果
2.8.1 预测系数a
a =1 至 10 列1.0000 -1.8904 1.3882 -0.1940 -0.1323 -0.0557 -0.0015 0.0228 0.0354 0.058111 至 13 列0.0931 -0.4149 0.3676
- 含义:在 LPC(线性预测编码)中,预测系数 a 用于描述当前语音样本与过去若干个样本之间的线性关系。具体来说,当前样本的预测值 s ^ ( n ) \hat{s}(n) s^(n) 是过去 p 个样本值 s ( n − k ) s(n - k) s(n−k) ( k = 1 , 2 , ⋯ , p ) (k = 1, 2, \cdots, p) (k=1,2,⋯,p)的线性组合,即 s ^ ( n ) = ∑ k = 1 p a k s ( n − k ) \hat{s}(n)=\sum_{k = 1}^{p}a_{k}s(n - k) s^(n)=∑k=1paks(n−k)。这里的预测阶数 p 为 12 ,所以 a 包含 13 个元素(第一个元素固定为 1 )。
- 数值分析:系数的正负和大小反映了过去样本对当前样本预测的贡献程度和方向。例如,系数为负表示该过去样本与当前样本的变化趋势相反,系数绝对值越大表示该过去样本对当前样本的影响越大。在这个结果中,a(2) 为 -1.8904 ,绝对值较大,说明前一个样本对当前样本的预测有比较重要的影响,且影响方向相反。
2.8.2 均方误差(MSE)
- 含义:均方误差(MSE)是衡量预测准确性的一个重要指标,它表示预测误差平方的平均值。计算公式为 M S E = 1 N ∑ n = 1 N e 2 ( n ) MSE=\frac{1}{N}\sum_{n = 1}^{N}e^{2}(n) MSE=N1∑n=1Ne2(n),其中 e ( n ) = s ( n ) − s ^ ( n ) e(n)=s(n)-\hat{s}(n) e(n)=s(n)−s^(n) 是预测误差,N 是样本数量。
- 分析:MSE 的值越小,说明预测信号与原始信号越接近,预测效果越好。在这个例子中,MSE 为 0.0296 ,可以根据具体的应用场景和经验来判断这个值是否合理。如果与其他预测方法或不同参数设置下的 MSE 进行比较,能更直观地评估当前预测的优劣。
2.8.3 信噪比(SNR)
snr_value = 17.0505
- 含义:信噪比(SNR)表示信号中有用信号与噪声(这里的噪声可以理解为预测误差)的比例,通常用分贝(dB)表示。计算公式为: S N R = 10 log 10 ∑ n = 1 N s 2 ( n ) ∑ n = 1 N e 2 ( n ) SNR = 10\log_{10}\frac{\sum_{n = 1}^{N}s^{2}(n)}{\sum_{n = 1}^{N}e^{2}(n)} SNR=10log10∑n=1Ne2(n)∑n=1Ns2(n)
- 分析:SNR 值越大,说明信号质量越好,预测误差相对较小。一般来说,在语音处理中,SNR 大于 10 dB 可以认为预测效果尚可,但具体的评判标准还需要结合实际应用和需求。在这个例子中,SNR 为 17.0505 dB ,表明预测信号与原始信号有一定的相似性,但仍存在一定的误差。
3 小结
LPC的核心是线性预测算法,通过理论推导,代码仿真,使我们看到了算法的有效性。那么线性预测编码的效率和波形编码相比有什么优势呢?