LPC线性预测谱及其语谱图绘制

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

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

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

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

移向后会发现，两个式子中的G与预测误差kesai的求法相同，由此我们就可以利用自相关函数预测的误差来表示模型增益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

下边附上完整程序代码：

%exp1_enframe：分帧加窗函数(for循环)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%输入参数：
% data              语音数据序列
% fs                采样率
% frame_time        帧时长
% frame_shift_time  帧移时长
% win               加窗类型

%输出参数：
% frame_m       分帧矩阵
% frame_w       分帧加窗后矩阵
% frame_length  帧长（每帧点数）
% frame_shift   帧移
% frame_number  帧数

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

function[frame_m,frame_w,frame_length,frame_shift,frame_number]= enframe(data,fs,frame_time,frame_shift_time,win)
if size(data, 1) ~= 1      % 确保data矩阵为单行矩阵
        data = data';
end
N=length(data);%获取数据长度
ts=1/fs;%获取采样间隔
frame_length=ceil((frame_time/ts+1));%每帧的点数（帧长）
frame_shift=ceil((frame_shift_time/ts+1));%相邻两帧起始点相差的点数（帧移）
frame_number=ceil((N+frame_shift-frame_length)/frame_shift);%分帧总数
N1=frame_number*frame_length;%分帧后数据总长度
N0=N1-N;%补零数
data=[data,zeros(1,N0)];%补零加到原数据序列后
frame_m=zeros(frame_length,frame_number);%初始化分帧矩阵
for i=1:frame_number
    frame_m(:,i)=data((i-1)*frame_shift+1:(i-1)*frame_shift+frame_length);
end
%将分帧矩阵按列赋值
w=str2func(win);%字符串转换窗函数
frame_w=zeros(frame_length,frame_number);%初始化分帧加窗后的矩阵
w=w(frame_length);      %产生对应窗函数数据序列 
for i=1:frame_number
    frame_w(:,i)=frame_m(:,i).*w;%加指定窗
end
%加窗运算赋值到矩阵

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

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% 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帧
N=1024;               %设置fft运算点数
f = (0:N/2-1)/N*fs;  %计算频率序列
fframe22=fft(frame22,N); %将加矩形窗第22帧数据由时域转换到频域
w=10*log10(abs(fframe22(1:end/2)).^2+eps);%计算功率谱
plot(f,w,'r');      %画出第22帧短时功率谱图
hold on
p=14;%设定 LPC 预测器阶
a4=lpc(frame22,p);%求 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');
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个数据相乘。可以将增益序列拓展至fft点数（N）列，这样就可以表达为增益矩阵与频率响应矩阵的点乘：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 = 10*log10(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=10*log10(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)')        %  设置纵轴标签

​