LPC语谱图绘制

利用matlab自带的lpc函数可以很方便的求出指定阶数的线性预测器系数。再由全极点模型下的传输函数表达式：

LPC谱就是线性预测系数的傅里叶变换的倒数再乘上模型增益G。（由matlab求出的线性预测系数包含了a0=1，如下图：)

传递函数写成差分方程式为：

而使用自相关函数法的预测误差的表达式为：

移向后会发现，两个式子中的G的求法相同，由此我们就可以利用自相关函数预测的误差来表示模型增益G。使用matlab自带的自相关法解线性预测系数的函数levinson函数。f这个函数的形式如下：
[Lpc,e]=levinson(r,p);输入r自相关函数序列，预测器阶数p，输出线性预测器系数Lpc,和预测的误差e最小预测均方误差（）为：
这里需要注意的是，matlab里边变量的索引是从1开始的，线性预测系数从数组中第二个开始（第一个为a0=1），因此在写程序的时候具体为：lpcError=r(1)-Lpc(2:end)r(2:15);
r是自相关函数，Lpc是线性预测器系数根据预测误差估计模型增益G=squr（lpcError）.
由此求出频率响应为：Y=sqrt(lpcError)./fft(Lpc,N);%求频率响应
U=10log10(abs(Y(1:end/2)).^2+eps);%转换单位为dB

下边附上完整程序代码：

% Experiment 3 线性预测谱估计% Apr.11 2020%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~clearclose allclc %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~% LPC谱load('ah.mat','data');fs=10000;[frame_m,frame_w,frame_length,frame_shift,frame_number]=enframe(data,10000,20e-3,10e-3,'hamming');% 调用分帧加窗函数，获取分帧矩阵、分帧加窗矩阵、帧长、帧数frame22=frame_w(:,22);       %取出第22帧(加窗后)f22=frame_m(:,22);       %取出第22帧(加窗前)N=1024;               %设置fft运算点数f = (0:N/2-1)/N*fs;  %计算频率序列fframe22=fft(frame22,N); %将加汉明窗第22帧数据由时域转换到频域ff22=fft(f22,N); %将加矩形窗第22帧数据由时域转换到频域w=10*log10(abs(fframe22(1:end/2)).^2+eps);%计算功率谱w2=10*log10(abs(ff22(1:end/2)).^2+eps);%计算功率谱plot(f,w,'r');      %画出第22帧短时功率谱图hold onp=14;%设定 LPC 预测器阶[r,lg]=xcorr(frame22,'biased');%求自相关函数序列r(lg<0)=[];[Lpc,e]=levinson(r,p);%求出预测器系数和预测误差lpcError=r(1)-Lpc(2:end)*r(2:15);%最小均方误差Y=sqrt(lpcError)./fft(Lpc,N);%求出频率响应U=10*log10(abs(Y(1:end/2)).^2+eps);%求功率谱衰减plot(f,U)xlabel('频率/Hz  （hamming窗）');ylabel('衰减/dB'); title('LPC谱图');legend('短时功率谱','LPC谱','location','best') [r,lg]=xcorr(f22,'biased');%求自相关函数序列r(lg<0)=[];[Lpc,e]=levinson(r,p);%求出预测器系数和预测误差lpcError=r(1)-Lpc(2:end)*r(2:15);%最小均方误差Y=sqrt(lpcError)./fft(Lpc,N);%求出频率响应U=10*log10(abs(Y(1:end/2)).^2+eps);%求功率谱衰减figure();plot(f,w2);hold onplot(f,U)xlabel('频率/Hz  （矩形窗）');ylabel('衰减/dB'); title('LPC谱图');legend('短时功率谱','LPC谱','location','best') %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

其实通过单帧的LPC功率谱计算，使用for循环，对每帧进行相对应的求解，最后整合到一起，就可以得出语音的LPC语谱图，但是由于for循环的计算量比较大，这里借鉴索引分帧矩阵运算的思路，由矩阵运算，求解整个LPC频率响应。在matlab上help下lpc函数：X也可以是一个矩阵，如果是矩阵的化则对矩阵每一个单个的列进行线性预测，并将预测系数返回到每一行中。这里尤其要注意，它得到的预测系数是返回到每一行中，因此我们在利用预测系数fft到频域时就要注意，分帧矩阵按列排布的，所以我们就需要将返回的预测系数转置后再进行fft。a1=lpc(frame_w,14);%求 LPC系数pf1=fft(a1.’,N);%声道模型功率谱响应曲线

关于增益G的计算，由于G只在短时内可以估计成常数，因此需要对每一帧单独来求增益系数。增益系数是根据自相关序列计算，我查了下资料，没能找到按列进行求自相关序列的函数，所以还是选择了使用for循环，一帧一帧的求。求解方法和单帧的LPC谱一样，只是要把增益系数放在一个序列中：for i=1:frame_number
[r,lg]=xcorr(frame_w(:,i),‘biased’);
r(lg<0)=[];
[Lpc,e]=levinson(r,p);
lpcError(i)=sqrt(r(1)-Lpc(2:end)r(2:15));%预测最小均方误差序列，即增益系数的平方序列
end得出的误差序列要转换为增益序列与每一帧的频率响应相乘即频率响应的第x列，要与增益序列的第x个数据相乘。可以将增益序列拓展至帧长列，这样就可以表达为增益矩阵与频率响应矩阵的点乘：G=lpcError(ones(1,N)????;
pf1=G./pf1;
U1=10log10(abs(pf1(1: end/2,: )).^2+ eps);%求功率谱衰减至此LPC谱的频域参数就求完了。下边附上完整代码：% Experiment 3 绘制LPC语谱图
% Apr. 11 2020
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
clear
close all
clc

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
we = audioread(‘we.wav’);%读取.wav文件数据
save(‘we.mat’, ‘we’) ;%转换为.mat文件
we=filter([1 -0.98],1,we);%预加重
fs=16000;
N=1024;
p=14;
[frame_m,frame_w,frame_length,frame_shift,frame_number]=enframe(we,fs,10e-3,5e-3,‘hamming’);
% 调用分帧加窗函数，获取分帧矩阵、分帧加窗矩阵、帧长、帧数
timeaxis=(1:frame_number)10e-3; %将帧数转换为离散时间轴
frameSpectrum = fft(frame_w, N); % 计算FFT
spectrumGram = 10log10(abs(frameSpectrum(1: end/2,: )).^2+ eps);
% 计算功率谱
f = (0: N/2-1)*fs/N; % 设置语谱图频率刻度

a1=lpc(frame_w,14);%求 LPC系数
aa=levinson(frame_w,14);
lpcError=[];
for i=1:frame_number
[r,lg]=xcorr(frame_w(:,i),‘biased’);
r(lg<0)=[];
[Lpc,e]=levinson(r,p);
lpcError(i)=sqrt(r(1)-Lpc(2:end)r(2:15));
end
pf1=fft(a1.’,N);%声道模型功率谱响应曲线
G=lpcError(ones(1,N)????;
pf1=G./pf1;
U1=10log10(abs(pf1(1: end/2,: )).^2+ eps);%求功率谱衰减

imagesc(timeaxis, f, spectrumGram); % 绘制语谱图
axis xy % 设置语谱图的坐标原点为左下角
colormap(jet) % 设置语谱图配色模式
colorbar % 设置语谱图色彩对照条
title(‘we 16K 语谱图’);
xlabel('时间 (s) (宽带) ') % 设置横轴标签
ylabel(‘频率 (Hz)’) % 设置纵轴标签

figure();
imagesc(timeaxis, f, U1); % 绘制语谱图
axis xy % 设置语谱图的坐标原点为左下角
colormap(jet) % 设置语谱图配色模式
colorbar % 设置语谱图色彩对照条
title(‘we 16K LPC语谱’);
xlabel('时间 (s) (宽带) ') % 设置横轴标签
ylabel(‘频率 (Hz)’) % 设置纵轴标签
​