语音信号线性预测分析

《视频语音处理技术》语音信号线性预测分析
学院名称：计算机与信息工程学院
专业名称：计算机科学与技术
年级班级：
姓名：
学号：
计算机与信息技术学院综合性、设计性实验报告
一、实验目的：
综合采用各种线性预测分析的方法，能够达到预测更为准确。

要求掌握各种下列语音信号线性预测分析技术，提高学生数字语音信号处理的能力。

利用MATLAB 编程环境和强大的处理功能，实现语音信号线性预测。

主要训练如下的项目从而获得线性预测的综合能力：
1、LPC 方程的自相关解法。

2、LPC 参数到LSP 参数的转换。

3、LSP 参数到LPC 参数的转换。

4、LPC 参数到ISP 参数的转换。

二、实验仪器或设备：w indowsXP 下的Matlab 编程环境 三、总体设计（设计原理、设计方案及流程等）
线性预测编码原理：利用过去的样值对新样值进行预测，然后将样值的实际值与其
预测值相减得到一个误差信号，显然误差信号的动态范围远小于原始语音信号的动态范围，对误差信号进行量化编码，可大大减少量化所需的比特数，使编码速率降低。

1.LPC 方程的自相关解法
利用对称托普利兹(Toeplitz)矩阵的性质，自相关法求解可用Levinson-Durbin （莱文逊-杜宾）递推算法求解。

该方法是目前广泛采用的一种方法。

利用Levinson-Durbin 算法递推时，从最低阶预测器开始，由低阶到高阶进行逐阶递推计算。

自相关法递推过程如下:
p
i E j i r a
i r k i i j i j
i ≤≤⎥⎦
⎤
⎢⎣
⎡
-=-∑= 1 )(-)( )1(1
-11)
-(，()() 0 0r E =)
1(2)1(--=i i i E k E
联立左面5式可对i=1、2…、p 进行递推求解，其最终解为 对于p 阶预测器，在上述求解预测器系数的过程中，阶数低于p 的各阶预测器系数也同时得到。

2、LPC 参数到LSP 参数的转换。

将P(z)和Q(z)中与LSP 系数无关的两个实根取得到如下两个新的多项式
从LPC 系数到LSP 系数的转换过程，其实就是求解使以上两式等于零时的 的值。

可采用以下几种方式求解。

第一种方法：利用代数方程式求解。

由 令 根据
得
第二种方法：离散傅立叶变换（DFT ）方法
对
和 系数求离散傅立叶变换，得到
各点的值，搜索最小值的位置，即是零点所在。

由于除了0和π之外，总共有p 个零点，
)(i i i
k a =()
11
- 1)1()(-≤≤=---i j a k a a
i j i i i j i j
，() 1 p j
j
a a j p =≤≤，)cos 21()1)(1(1)()(212
1
1111--=---=-+-=--=+='∏∏z z e z e z z z P z P p i i j j p i i
i ωωω)cos 21()1)(1(1)()(212
111211
--=---=-+-=--=-='∏∏z z e z e z z z Q z Q p i i j j p i i
i θθθi i θωcos cos 、]
cos )(5.0[2)
5.0cos 5.0(2cos 21111121i i i z z z z z z z z ωωω-+=+-=+-------ω
j e z =ω
ωωsin cos j e j +=x
z z 2cos 21==+-ω)(z P ')(z Q ' exp(，N
jk z k π-=110-=N k 、、、
而且 和 的根是相互交替出现的，因此只要很少的计算量即可解得，其中N
的取值取64～128就可以
第三种方法：利用切比雪夫(Chebyshev)多项式求解
则 可以写作
其中 是m 阶的Chebyshev 多项式 。

f(i)是由递推关系计算得到的 和 的每个系数。

用下面的递推关系
其中 。

多项式C(x)在x=cos ω时的递推关系是：
其中初始值
第四种方法：
将0～π之间均分为60个点，以这60个点的频率值代入(6-41)、(6-42)式，检查它们的符号变化，在符号变化的两点之间均分为4份，再将这三个点频率值代入方程(6-41)、(6-42)，符号变化的点即为所求的解。

这种方法误差略大，计算量较大，但程序实现容易。

3、LSP 参数到LPC 参数的转换
已知量化和内插的LSP 系数 ，可用下式计算 和 的系数 和
： )(z P ')(z Q 'ω
j e z
=)(z P ')(z Q ')
(2)(2x C e z P p
j ω
-=')
(2)(2x C e z Q p
j θ
-='2
)2()()12
(
)()2()()1()()(12
21
22
p
f x T p f x T f x T f x T x C p p p +-++++=-- mx
x T m cos )(=)(z P ')(z Q '2
1 0)()1()
()1(212111p i i f a a i f i f a a i f i p i i p i ，，， =⎪⎩⎪⎨⎧+-=+-+=+-+-+0.1)0()0(21==f f 2
)2()(end
)
2
(2112
for 2121p
f x x C k p f x to p
k k k k +-=-+-=-=
++λλλλλ12
=p λ01
2
=+p λ)(z P ')(z Q '()p i ')(i q '
以下的递推关系可利用q i ，i =0、1…、p -1，来计算
其中 ，初始值
把上面递推关系中的q2i-1替换为q2i ，就可以得到
根据 ，
和
和 利用
得到LPC 系数为
4、LPC 参数到ISP 参数的转换。

语音编码时，可将LPC 系数转换为ISP 系数以进行量化和内插。

LPC 系数与ISP 系数之间的转换与LSP 类似。

从LPC 转换为ISP 系数时，首先应用求解LSP 参数的方法求解出前p -1个ISP 系数，再给第p 个参数赋上合适的值，即可得到ISP 系数。

解码时，首先根据量化ISP 系数得到p -1个LPC 系数，再根据第p 个ISP 系数得到最后一个LPC 系数。

四、实验步骤（包括主要步骤、代码分析等）
1.下面给出从LPC 参数到LSP 参数转换的MATLAB 程序：
a_lsf_main.m
% 已知语音文件求出其LPC 系数后，调用a_lsfLconversation.m 函数求其对应的LSF % a_lsf_main.m
)cos 21()1)(1(1)()(212
111211
--=---=-+-=--=+='∏∏z z e z e z z z P z P p i i j j p i i
i ωωω)cos 21()1)(1(1)()(2121
1111
--=---=-+-=--=-='∏∏z z e z e z z z Q z Q p i i j j p i i
i θθθ)(i p '2121for 1 to 2
()2(1)2(2) for 1 to 1
()()2(1)(2) end end
i i i p p i q p i p i j i p j p j q p j p j --='''=--+-=-''''=--+-1212cos --=i i q ω1)0(='p 0
)1(=-'p )
(i q '()p i ')(i q ')(z P ')(z Q ')(z P )(z Q 11()()(1) 1 2 2
()()(1) i 1 2 2p i p i p i i p q i q i q i p ''=
+-=⎧⎨
''=--=⎩，，，，，
，，，)]()([2
1
)(z Q z P z A +=11110.5()0.5() 1 2 20.5(1)0.5(1) 2 1 2 2 i p i q i i p a p p i q p i i p p p +=⎧=⎨
+--+-=++⎩，
，，，，
，，，
clear;close all% 将所有变量置为0
clc% 清除命令窗口
fid=fopen('sx86.txt','r');
p1=fscanf(fid,'% f');
fclose(fid);
p=filter([1 -0.68],1,p); % 预加重滤波
x=fra(320,160,p); % 分帧，帧移为160个样点
x=x(60,:); % 取第60帧作为分析帧
N=16% 给线性预测分析的阶次赋值
a1=lpc(x,N);% 调用MATLAB库函数中的lpc函数求解出LPC系数a1
% 此处也可以调用本章赋的函数lpc_coeffi-
% cients（s，p），调用语句为a1=lpc_coefficients(x,N)
a=a1(:);% 将线性预测系数a1赋给矩阵a
lsf=a_lsf_conversion(a);% 调用函数a_lsf_conversion实现从LPC系数到LSF参数 % 转换
% lsf=poly21lsf(a);% 也可调用MATLAB库函数中的poly21lsf(a)函数求解出LSF系数，调用
% 结果为归一化角频率
lsf_abnormalized=lsf.*(6400/3.14);% 将求得的lsf参数反归一化到
% 0～6400Hz
% 使用时可根据实际需要进行更改，如窄带语言编码语音信号频带范围为300～3400Hz，此时
% 就需要将6400Hz改为3400Hz
% 将求得的归一化、反归一化lsf参数输出到文本文件：从lpc系数解得的lsf参数.txt fid= fopen('从lpc系数解得的lsf参数.txt','w');
fprintf(fid,'归一化的lsf：\n');
fprintf(fid,'%6.2f,',lsf);
fprintf(fid,'\n');
fprintf(fid,'反归一化的lsf:\n');
fprintf(fid,'% 8.4f,',lsf_abnormalized);
fclose(fid);
% EOF a_lsf_main.m
函数 a_lsf_conversion 的MATLAB程序见 a_lsf_conversion.m.
% 程序 a_lsf_conversion.m.
function lsf=a_lsf_conversion(a)
%如果a不是实数，输出错误信息：LSF不适用于复多项式的求解
if~isreal(a),
error('Line spectral frequencies are not defined for complex polynomials.');
end
% 如果a(1)不等于1，将其归一化为1 if a(1)~=1.0,
a=a./a(1);% 将矩阵a 的每个元素除以a(1)在赋给矩阵a end
% 判断线性预测多项式的跟是否都在单位园内，如果不在，则输出错误信息 if (max(abs(roots(a)))>=1.0),
error('The polynomil must have all roots inside of the unit circle.'); end
% 求对称和反对称多项式的系数
p=length(a)-1; % 求对称和反对称多项式的阶次 a1=[a;0]; % 给行矩阵a 再增加一个元素为0的行
a2=a1(end:-1:1); % a2的第一行为a1的最后一行，最后一行为a1的第一行 p1=a1+a2; % 求对称多项式的系数 Q1=a1-a2; % 求反对称多项式的系数
% 如果阶次p 为偶数次，从p1取掉实数根z=-1，从Q1去取掉实数根z=1 % 如果阶次p 为奇数次，从Q1去取掉实数根z=1及z=-1
if rem(p,2),% 求解p 除以2的余数，如果p 为奇数次，余数为1，否则为0 Q=deconv(Q1,[1 0 -1]);% 奇数阶次，从Q1取掉实数根z=1 p=p1;
else % p 为偶数阶次执行下面操作 Q=deconv(Q1,[1 -1]);% 从Q1取掉实数根z=1 p=deconv(P1,[1 1]);% 从P1取掉实数根z=-1 end
rp=roots(P);% 求去掉实根后的多项式P 的根 rQ=roots(Q);% 求去掉实根后的多项式Q 的根
aP=angle(rp(1:2:end));% 将多项式P 的根转换为角度(为归一化角频率）赋给ap aQ=angle(rQ(1:2:end));% 将多项式Q 的根转换为角度（为归一化角频率）赋给aQ lsf=sort([aP;aQ]);% 将P 、Q 的根（归一化角频率）按从小到大顺序排序后即为lsf % EOF a_lsf_conversion.m
其中a_lsf_conversion.m 为求解LSF 的函数，a_lsf_main.m 为主程序。

由于MATLAB
程序本身有求多项式根的函数，因此在求解
和 零点时直接调用即可，这极大简化了求解过程。

2.
LPC 转换为ISP 系数的MATLAB 程序和ISP 系数转换为LPC 的MATLAB 程序
a_isf_lpc_conversion.m
%已知语音文件求出其LPC 系数后，求其对应的ISF ，再将ISF 转换为a 。

其中isf 为转换后的ISF 值，a 为isf 转换后的lpc 系数
)(z P ')(z Q '
clear;close all
clc
fid=fopen('sx86.txt','r');
p1=fscanf(fid,'% f')
fclose(fid);
p=filter([1-0.68],1,p1);%预加重滤波
x=fra(320,160,p);
x=x(60,:);
a3=lpc(x,15);
a4=a3(:);%将线性预测系数赋给矩阵a
lsf=a_lsf_conversion(a4);%调用函数a_lsf_conversion实现从LPC系数到lsf参数的转换
%函数a_lsf_conversion()的MATLAB程序见本章6.6.2节所赋程序，此处也可调用MATLAB 自带的
%库函数lsf=poly21sf(a4);
isf=lsf;
isf(16,1)=0.5*acos(a4(16,1));
%isp的最后一个参数取为a的最后一个参数，isf的最后一个参数取为0.5*acos(a4(16,1))
%下面是从isf求a的程序，其中前p-1个a参数根据前p-1个isf参数得到，最后一个a 参数根据isf的第p个参数得到
isf1=isf(1:(size(isf)-1),:);%将isf前p-1个参数赋给isf1
a2=lsf_lpc_conversion(isf1)%调用函数lsf_lpc_conversion实现从lsf参数到LPC系数的转换
a2(1,16)=cos(2*isf(16,1));%最后一个a参数根据isf的第p个参数得到
a=a2;%将转换得到的lpc系数赋给a
%EOFa_isf_lpc_conversion.m
a_isf_lpc_conversion_20frame.m
clear;close all
clc
fid=fopen('sx86.txt','r');
p1=fscanf(fid,'%f')
fclose(fid);
p=filter([1-0.68],1,p1);%预加重滤波
x1=fra(320,160,p);
for i=60:79
x=x1(i,:)
a3=lpc(x,15);
a4=a3(:);%将线性预测系数付给矩阵a4
lsf=a_lsf_conversion(a4)%调用函数a_lsf_conversion实现从LPC系数到lsf参数的转换
isf=lsf;
isf(16)=0.5*acos(a4(16,1));%isp的最后一个参数取为a的最后一个参数，isf的最后一个参数取为0.5*acos(a4(16,1)）
isp=cos(isf);
hold on%让连续20帧isp及isf绘制在一个图形中
figure(1);
for j=1:16
isf2(i-59,j)=isf(j);
end
figure(2);
plot(isp)
end
%EOF a_isf_lpc_conversion_20frame.m
运行结果：
连续20帧语音信号的ISP轨迹图
连续20帧语音信号的ISF轨迹图
五、结果分析与总结
1.线性预测分析的基本思想是：一个语音取样的现在值可以用若干个语音取样过去值的加权线性组合来逼近。

线性预测基本原理是建立在语音的数字模型基础上，为估计数字模型中的参数，线性预测法提供了一种可靠精确而有效的方法。

2.在上机实验的过程中，应注意避免一些小细节处的失误，提高效率。

教师签名：
年月日
河南师范大学计算机与信息技术学院。