《MATLAB语音信号分析与合成（第二版）》：第4章 语音信号的线性预测分析

• 前言
• 1. 数据与函数路径设置
• 2. MATLAB仿真一：LPC的频谱分析
• 3. MATLAB仿真二：普通预测法和格型预测法功率谱响应的比较
• 4. MATLAB仿真三：FFT频谱和LPC谱的比较
• 5. MATLAB仿真四：LPCC参数匹配
• 6. MATLAB仿真五：LPCC参数匹配
• 小结

前言

《MATLAB语音信号分析与合成（第二版）》是中科院声学所的大佬宋知用老师数十年经验积累下的呕心之作，对于语音信号处理相关感兴趣的同学，日后希望在语音信号分析、处理与合成相关领域进行一定研究的话，可以以此进行入门。

语音信号处理是数字信号处理的一个重要分支。本书含有许多数字信号处理的方法和 MATLAB函数。 全书共10章。第1-4章介绍语音信号处理的一些基本分析方法和手段,以及相应的MATLAB函数;第5-9章介绍语音信号预处理和特征的提取,包括消除趋势项和基本的减噪方法,以及端点检测、基音的提取和共 振峰的提取,并利用语音信号处理的基本方法,给出了多种提取方法和相应的 MATLAB程序;第10章结合 各种参数的检测介绍了语音信号的合成、语音信号的变速和变调处理,还介绍了时域基音同步叠加( TD PSOLA)的语音合成,并给出了相应的MATLAB程序。附录A中给出了调试复杂程序的方法和思路。 本书可作为从事语音信号处理的本科高年级学生、研究生或科研工程技术人员的辅助读物,也可作为从 事信号处理研究与应用的科研工程技术人员的参考用书。

我的研究生导师的主攻方向就是语音信号处理相关，虽然自己研究生期间的大论文方向是数字图像处理，但所谓语音图像不分家，自己在老师的研究生主讲课小波变换上虽然划水，但在后期导师的语音信号处理的课程设计和工程应用上自己在语音上还算入了一点小门道，在结课测试中拿到了小组第一，导师还特地发了三百大洋的伙食经费以资鼓励。

这次重新捡起语音识别，正好入手了宋老师的这本书，算是自己重新复习一遍吧，主要以介绍各章节中源码为主，这是本书的第四章的5个仿真应用实例，话不多说，开始！

1. 数据与函数路径设置

书中经常会调用的一些函数(自编函数或取自其他应用工具箱中的函数)已集中在basic_tbx工具箱中，在运行本书的程序前请把该工具箱设置(用set path设置)在工作路径下；

当要运行EMD处理时，要把emd工具箱设置在工作路径下;

当要运行主体延伸基音检测时，要把Pitch_ztlib工具箱设置在工作路径下;

当要进行时域基音同步叠加语音合成时，要把psola_lib工具箱设置在工作路径下;

当要应用本书提供的语音数据时，最好把speech_signal设置在工作路径下。

本书的所有函数和程序都在MATLAB R2009a版本下调试通过。（我用的是MATLAB2015b，有些函数已经更新，所以我会进行修改，以便调试通过）

路径设置的方法如下：

打开MATLAB，点击“主页”，找到设置路径

将上述文件夹路径全部添加到MATLAB搜索路径中

添加完毕，保存，开始仿真。

2. MATLAB仿真一：LPC的频谱分析

%% pr4_2_1 clear all; clc; close all;filedir=[]; % 设置数据文件的路径filename='aa.wav'; % 设置数据文件的名称fle=[filedir filename]    % 构成路径和文件名的字符串% [x,fs]=wavread(fle);      % 读入语音数据[x,fs]=audioread(fle);      % 读入语音数据L=240;      % 帧长y=x(8001:8000+L);  % 取一帧数据p=12;% LPC的阶数ar=lpc(y,p);% 线性预测变换Y=lpcar2ff(ar,255);% 求LPC的频谱值est_x=filter([0 -ar(2:end)],1,y);% 用LPC求预测估算值err=y-est_x;% 求出预测误差fprintf('LPC:\n');fprintf('%5.4f   %5.4f   %5.4f   %5.4f   %5.4f   %5.4f   %5.4f\n',ar);fprintf('\n');% 作图pos = get(gcf,'Position');set(gcf,'Position',[pos(1), pos(2)-200,pos(3),pos(4)+150]);subplot 311; plot(x,'k'); axis tight;title('元音/a/波形'); ylabel('幅值')subplot 323; plot(y,'k'); xlim([0 L]); title('一帧数据'); ylabel('幅值')subplot 324; plot(est_x,'k'); xlim([0 L]); title('预测值'); ylabel('幅值')subplot 325; plot(abs(Y),'k'); xlim([0 L]); title('LPC频谱'); ylabel('幅值'); xlabel('样点')subplot 326; plot(err,'k'); xlim([0 L]); title('预测误差'); ylabel('幅值'); xlabel('样点')

fle =aa.wavLPC:1.0000   -0.9924   -0.6898   1.1805   0.7087   -1.3034   -0.23030.5580   0.3657   -0.4033   -0.1244   0.0074   0.1376  

3. MATLAB仿真二：普通预测法和格型预测法功率谱响应的比较

%% pr4_3_1 clear all; clc; close all;filedir=[]; % 设置数据文件的路径filename='aa.wav'; % 设置数据文件的名称fle=[filedir filename];   % 构成路径和文件名的字符串% [x,fs]=wavread(fle);      % 读入语音数据[x,fs]=audioread(fle);      % 读入语音数据L=240;      % 帧长p=12;% LPC的阶数y=x(8001:8240+p);  % 取一帧数据[EL,alphal,GL,k]=latticem(y,L,p);% 格型预测法ar=alphal(:,p);a1=lpc(y,p);% 普通预测法Y=lpcar2pf(a1,255);% 把a1转成功率谱Y1=lpcar2pf([1; -ar],255);% 把ar转成功率谱fprintf('AR1系数(格型预测法):\n');fprintf('%5.4f   %5.4f   %5.4f   %5.4f   %5.4f   %5.4f\n',-ar);fprintf('AR2系数(普通预测法):\n');fprintf('%5.4f   %5.4f   %5.4f   %5.4f   %5.4f   %5.4f\n',a1(2:p+1));% 作图m=1:257;freq=(m-1)*fs/512;plot(freq,10*log10(Y),'k'); grid;line(freq,10*log10(Y1),'color',[.6 .6 .6],'linewidth',2);legend('普通预测法','格型预测法'); ylabel('幅值/dB');title('普通预测法和格型预测法功率谱响应的比较'); xlabel('频率/Hz'); 

AR1系数(格型预测法):-1.0017   -0.7454   1.2245   0.7761   -1.3801   -0.31590.6597   0.4127   -0.4792   -0.1459   0.0540   0.1303AR2系数(普通预测法):-1.0284   -0.6204   1.1818   0.5881   -1.2736   -0.09300.4840   0.2729   -0.3449   -0.0752   -0.0308   0.1319

4. MATLAB仿真三：FFT频谱和LPC谱的比较

%% pr4_4_1 clear all; clc; close all;filedir=[]; % 设置数据文件的路径filename='aa.wav'; % 设置数据文件的名称fle=[filedir filename]    % 构成路径和文件名的字符串% [x,fs]=wavread(fle);      % 读入语音数据[x,fs]=audioread(fle);      % 读入语音数据L=240;      % 帧长p=12;% LPC的阶数y=x(8001:8000+L);  % 取一帧数据ar=lpc(y,p);% 线性预测变换nfft=512;   % FFT变换长度W2=nfft/2;m=1:W2+1;   % 正频率部分下标值Y=fft(y,nfft);     % 计算信号y的FFT频谱Y1=lpcar2ff(ar,W2-1);     % 计算预测系数的频谱zz=lpcar2zz(ar);   % 计算预测系数的根值for k=1 : 12    fprintf('%4d   %5.6f   %5.6f\n',k,real(zz(k)),imag(zz(k)));end% 作图subplot 211; plot(y,'k');title('一帧语音信号的波形'); ylabel('幅值'); xlabel('(a)')subplot 212; plot(m,20*log10(abs(Y(m))),'k','linewidth',1.5); line(m,20*log10(abs(Y1)),'color',[.6 .6 .6],'linewidth',2)axis([0 W2+1 -30 25]); ylabel('幅值/db');legend('FFT频谱','LPC谱',3); xlabel(['样点' 10 '(b)'])title('FFT频谱和LPC谱的比较');

fle =aa.wav   1   0.859081   0.140565   2   0.859081   -0.140565   3   0.597665   0.759588   4   0.597665   -0.759588   5   0.763608   0.561950   6   0.763608   -0.561950   7   -0.843682   0.363962   8   -0.843682   -0.363962   9   -0.594580   0.516397  10   -0.594580   -0.516397  11   -0.285869   0.575662  12   -0.285869   -0.575662

5. MATLAB仿真四：LPCC参数匹配

%% pr4_4_2 clear all; clc; close all;% [x1,fs]=wavread('s1.wav');      % 读入信号s1% x2=wavread('s2.wav');    % 读入信号s2% x3=wavread('a1.wav');    % 读入信号a1[x1,fs]=audioread('s1.wav');      % 读入信号s1x2=audioread('s2.wav');    % 读入信号s2x3=audioread('a1.wav');    % 读入信号a1wlen=200;  % 帧长inc=80;    % 帧移x1=x1/max(abs(x1));      % 幅值归一化x2=x2/max(abs(x2));x3=x3/max(abs(x3));p=12;      % LPC阶数[DIST12,y1lpcc,y2lpcc]=lpcc_dist(x1,x2,wlen,inc,p);% 计算x1与x2的LPCC距离[DIST13,y1lpcc,y3lpcc]=lpcc_dist(x1,x3,wlen,inc,p);% 计算x1与x3的LPCC距离% 作图figure(1)plot(y1lpcc(3,:),y2lpcc(3,:),'k+'); hold onplot(y1lpcc(7,:),y2lpcc(7,:),'kx'); plot(y1lpcc(12,:),y2lpcc(12,:),'k^');plot(y1lpcc(16,:),y2lpcc(16,:),'kh'); legend('第3帧','第7帧','第12帧','第16帧',2)title('/i1/与/i2/之间的LPCC参数匹配比较')xlabel('信号x1');ylabel('信号x2')axis([-6 6 -6 6]);line([-6 6],[-6 6],'color','k','linestyle','--');figure(2)plot(y1lpcc(3,:),y3lpcc(3,:),'k+'); hold onplot(y1lpcc(7,:),y3lpcc(7,:),'kx'); plot(y1lpcc(12,:),y3lpcc(12,:),'k^');plot(y1lpcc(16,:),y3lpcc(16,:),'kh'); legend('第3帧','第7帧','第12帧','第16帧',2)title('/i1/与/a1/之间的LPCC参数匹配比较')xlabel('信号x1');ylabel('信号x3')axis([-6 6 -6 6]);line([-6 6],[-6 6],'color','k','linestyle','--');

6. MATLAB仿真五：LPCC参数匹配

%% pr4_5_1 clear all; clc; close all;filedir=[];     % 设置数据文件的路径filename='aa.wav';     % 设置数据文件的名称fle=[filedir filename] % 构成路径和文件名的字符串% [x,fs]=wavread(fle);   % 读取语音文件aa.wav[x,fs]=audioread(fle); % 读入语音数据x=x/max(abs(x));% 幅值归一化time=(0:length(x)-1)/fs;      % 求出对应的时间序列N=200;   % 设定帧长M=80;    % 设定帧移的长度  xn=enframe(x,N,M)';    % 按照参数进行分帧s=xn(:,100);    % 取分帧后的笫100帧进行分析p=12;    % 设预测阶次num=257; % 设定频谱的点数a2 =lpc(s,p);   % 利用信号处理工具箱中的函数lpc求预测系数a2Hw=lpcar2ff(a2,num-2); % 调用lpcar2ff函数从预测系数a求出LP谱HwHw_abs=abs(Hw); % 取Hw的模值lsf=ar2lsf(a2); % 调用ar2lsf函数把ar系数转换的lsf参数P_w=lsf(1:2:end);      % 用lsf求出P和Q对应的频率，单位为弧度Q_w=lsf(2:2:end);P_f=P_w*fs/2/pi;% 转换成单位为HzQ_f=Q_w*fs/2/pi;figure(1)pos = get(gcf,'Position');    % 设置绘图框set(gcf,'Position',[pos(1), pos(2)-100,pos(3),(pos(4)-180)]);plot(time,x,'k');      % 画出信号的波形title('语音信号aa.wav的波形图 ');xlabel('时间/s'); ylabel('幅值')    xlim([0 max(time)]);figure(2)subplot 211; plot(s,'k');     % 画出一帧信号的波形title('语音信号aa.wav的一帧波形图 ');xlabel(['样点值' 10 '(a)']); ylabel('幅值')    freq=(0:num-1)*fs/512;   % 计算频域的频率序列m=1:num;K=length(Q_w);ar=lsf2ar(lsf); % 调用lsf2ar函数把lsf转换成预测系数ar Hw1=lpcar2ff(ar,num-2);% 调用lpcar2ff函数,从预测系数ar求出LP谱Hw1Hw1_abs=abs(Hw1);subplot 212;    % 把Hw和Hw1画在一个图中hline1 = plot(freq,20*log10(Hw_abs(m)/max(Hw_abs)),'k','LineWidth',2); hline2 = line(freq+1,20*log10(Hw1_abs(m)/max(Hw1_abs)),...    'LineWidth',5,'Color',[.6 .6 .6]);set(gca,'Children',[hline1 hline2]);axis([0 fs/2 -35 5]);title('语音信号的LPC谱和线谱对还原LPC的频谱 ');xlabel(['频率/Hz' 10 '(b)']); ylabel('幅值')    for k=1 : K     % 把P_f和Q_f也在图中用垂直线标出    line([Q_f(k) Q_f(k)],[-35 5],'color','k','Linestyle','--');    line([P_f(k) P_f(k)],[-35 5],'color','k','Linestyle','-');endfor k= 1 : p+1  % 显示预测系数a2和ar，对两者进行比较    fprintf('%4d   %5.6f   %5.6f\n',k,a2(k),ar(k));end

fle =aa.wav   1   1.000000   1.000000   2   -0.904775   -0.904775   3   -0.412674   -0.412674   4   0.909116   0.909116   5   0.261812   0.261812   6   -0.740536   -0.740536   7   0.131235   0.131235   8   0.007822   0.007822   9   0.128686   0.128686  10   -0.027247   -0.027247  11   -0.049160   -0.049160  12   -0.160520   -0.160520  13   0.109142   0.109142

小结

线性预测系数LPCC是语音信号的分析处理中一个最常见的参数，常与MFCC参数一起拿来进行对比分析，因为它们属于是语音信号匹配识别过程中的基本参数了，希望深入学习的同学，建议去研习书中第四章节的内容。后期会对其中一些知识点在自己理解的基础上进行讨论补充，欢迎大家一起学习交流。

关于宋老师：宋知用——默默传授MATLAB与信号处理知识的老人家

本系列文章列表如下：
《MATLAB语音信号分析与合成（第二版）》：第2章 语音信号的时域、频域特性和短时分析技术
《MATLAB语音信号分析与合成（第二版）》：第3章 语音信号在其他变换域中的分析技术和特性
《MATLAB语音信号分析与合成（第二版）》：第4章 语音信号的线性预测分析
《MATLAB语音信号分析与合成（第二版）》：第5章 带噪语音和预处理
《MATLAB语音信号分析与合成（第二版）》：第6章 语音端点的检测（1）
《MATLAB语音信号分析与合成（第二版）》：第6章 语音端点的检测（2）
《MATLAB语音信号分析与合成（第二版）》：第7章 语音信号的减噪
《MATLAB语音信号分析与合成（第二版）》：第8章 基音周期的估算方法
《MATLAB语音信号分析与合成（第二版）》：第9章 共振峰的估算方法
《MATLAB语音信号分析与合成（第二版）》：第10章 语音信号的合成算法