语音信号处理——线性预测编码LPC「建议收藏」

语音信号处理——线性预测编码LPC「建议收藏」语音信号处理二:干净语音的特征提取:今天的信号与系统,DSP知识点参考SpokenLanguageProcessing第5,6章LPC方程的Durbin算法推导:语音信号数字处理(杨行峻,迟惠生)第四章,数字语音处理(Rabiner)第九章作业是自己实现语音信号的LPC预测算法:输入一段语音信号,选定不同阶数p,在最小二乘准则下,用自相关法估计预测系数aia_iai​,对比重建语…

大家好,又见面了,我是你们的朋友全栈君。

语音信号处理二:干净语音的特征提取:

今天的信号与系统,DSP知识点参考 Spoken Language Processing 第5, 6 章

LPC方程的Durbin算法推导:语音信号数字处理(杨行峻,迟惠生)第四章,数字语音处理(Rabiner)第九章

作业是自己实现语音信号的LPC预测算法:输入一段语音信号,选定不同阶数p,在最小二乘准则下,用自相关法估计预测系数 a i a_i ai,对比重建语音和原始语音的时域&短时频谱差别

自相关法可以用普通的矩阵求逆,和Durbin算法做对比。

语音信号的生成模型:激励-滤波模型:

在这里插入图片描述
语音信号的激励部分:声门激励

  • 声门:选择激励源
  • 声带振动:周期性信号
    • 声带打开:产生激励
    • 声带闭合:信号为0
    • 对应:浊音/元音
    • 声带有个开闭的过程
  • 声带松弛:白噪声(与频率无关)
    • 时域:高斯分布
    • 短时频谱:均匀分布
      在这里插入图片描述

语音信号的滤波部分:无损声管模型

  • 声管:声道各个器官的抽象模型
    在这里插入图片描述
    H ( z ) = X ( z ) E ( z ) = 1 1 − ∑ k = 1 p a k z − k = 1 A ( z ) H(z)=\frac{X(z)}{E(z)}=\frac{1}{1-\sum_{k=1}^pa_kz^{-k}}=\frac{1}{A(z)} H(z)=E(z)X(z)=1k=1pakzk1=A(z)1

  • 鼻腔的作用:并行的通道

    • 简化:不考虑口鼻同时打开的情况
      在这里插入图片描述

声源是由声带产生的,声带向声道提供激励信号,这种激励可以是周期性的或非周期性的。当声带处于发声状态(振动)时,会产生有声声音(例如,元音);而当声带处于无声状态时,会产生无声声音(例如,辅音)。声道可以看作是一个滤波器,它可以对来自声带的激励信号频谱进行整形以产生各种声音。

线性预测编码(LPC):Linear Predictive Coding

LPC编码的基本思想:

“一个语音取样的现在值可以用若干个语音取样过去值的加权线性组合来逼近”,用过去p个样本点预测当前值:
x ~ [ n ] = ∑ k = 1 p a k x [ n − k ] \widetilde{x}[n]=\sum_{k=1}^pa_kx[n-k] x
[n]=
k=1pakx[nk]

在线性组合中的加权系数 a k a_k ak称为预测器系数。通过使实际语音抽样和线性预测抽样之间差值的平方和达到最小值,能够决定唯一的一组预测器系数。

预测误差:
e [ n ] = x [ n ] − x ~ [ n ] = x [ n ] − ∑ k = 1 p a k x [ n − k ] e[n]=x[n]-\widetilde{x}[n] = x[n] – \sum_{k=1}^pa_kx[n-k] e[n]=x[n]x
[n]=
x[n]k=1pakx[nk]

m个语音信号样本片段的周期延拓: x m [ n ] = x [ m + n ] x_m[n] = x[m+n] xm[n]=x[m+n]

短时预测误差:
E m = ∑ n e m 2 [ n ] = ∑ n ( x m [ n ] − x ~ m [ n ] ) 2 = ∑ n ( x m [ n ] − ∑ j = 1 p a j x m [ n − j ] ) 2 E_m=\sum_ne_m^2[n] = \sum_n\left(x_m[n]-\widetilde{x}_ m[n]\right)^2=\sum_n\left(x_m[n]-\sum_{j=1}^pa_jx_m[n-j]\right)^2 Em=nem2[n]=n(xm[n]x
m
[n])
2
=
n(xm[n]j=1pajxm[nj])2

线性预测编码通过估计共振峰、剔除它们在语音信号中的作用、估计保留的蜂鸣音强度与频率来分析语音信号。剔除共振峰的过程称为逆滤波,经过这个过程剩余的信号称为残余信号(residue)。

描述峰鸣强度与频率、共鸣峰、残余信号的数字可以保存、发送到其它地方。线性预测编码通过逆向的过程合成语音信号:使用蜂鸣参数与残余信号生成源信号、使用共振峰生成表示声道的滤波器,源信号经过滤波器的处理就得到语音信号。
J = E [ e 2 ( k ) ] = E [ ( s ( k ) − ∑ p = 1 P a p s ( k − p ) ) 2 ] J=E\left[e^{2}(k)\right]=E\left[\left(s(k)-\sum_{p=1}^{P}a_{p}s(k-p)\right)^{2}\right] J=E[e2(k)]=E(s(k)p=1Paps(kp))2

LPC特点:

  • LPC分析/AR模型

  • 把声道抽象为一个全极点模型:
    H ( z ) = X ( z ) E ( z ) = 1 1 − ∑ k = 1 p a k z − k = 1 A ( z ) H(z)=\frac{X(z)}{E(z)}=\frac{1}{1-\sum_{k=1}^pa_kz^{-k}}=\frac{1}{A(z)} H(z)=E(z)X(z)=1k=1pakzk1=A(z)1

  • p:级联声管个数,LPC分析阶数

  • 时域表示:
    x [ n ] = ∑ k = 1 p a k x [ n − k ] + e [ n ] x[n] = \sum_{k=1}^pa_kx[n-k]+e[n] x[n]=k=1pakx[nk]+e[n]

LPC分析的正交性原理:

  • 预测误差与当前样本正交:
    &lt; e m , x m i &gt; = ∑ n e m [ n ] x m [ n − i ] = 0 1 ≤ i ≤ p &lt;e_m,x_m^i&gt;=\sum_ne_m[n]x_m[n-i] =0 \quad\quad1≤i≤p <em,xmi>=nem[n]xm[ni]=01ip ∑ n x m [ n − i ] x m [ n ] = ∑ j = 1 p a j ∑ n x m [ n − i ] x m [ n − j ] i = 1 , 2 , . . . , p \sum_nx_m[n-i]x_m[n]=\sum_{j=1}^pa_j\sum_nx_m[n-i]x_m[n-j]\quad\quad i=1,2,…,p nxm[ni]xm[n]=j=1pajnxm[ni]xm[nj]i=1,2,...,p

  • 相关系数: ϕ m [ i , j ] = ∑ n x m [ n − i ] x m [ n − j ] \phi_m[i,j]=\sum_nx_m[n-i]x_m[n-j] ϕm[i,j]=nxm[ni]xm[nj]

  • Yule-Walker Equations: ∑ j = 1 p a j ϕ m [ i , j ] = ϕ m [ i , 0 ] i = 1 , 2 , . . . , p \sum_{j=1}^pa_j\phi_m[i,j]=\phi_m[i,0]\quad i=1,2,…,p j=1pajϕm[i,j]=ϕm[i,0]i=1,2,...,p

  • 预测误差:
    E m = ∑ n x m 2 [ n ] − ∑ j = 1 p a j ∑ n x m [ n ] x m [ n − j ] = ϕ [ 0 , 0 ] − ∑ j = 1 p a j ϕ [ 0 , j ] E_m=\sum_nx_m^2[n]-\sum_{j=1}^pa_j\sum_nx_m[n]x_m[n-j]=\phi[0,0]-\sum_{j=1}^pa_j\phi[0,j] Em=nxm2[n]j=1pajnxm[n]xm[nj]=ϕ[0,0]j=1pajϕ[0,j]

  • 预测误差的能量归一化:
    e m [ n ] = G u m [ n ] ∑ n u m 2 [ n ] = 1 E m = ∑ n e m 2 [ n ] = G 2 ∑ n u m 2 [ n ] = G 2 e_m[n]=Gu_m[n]\quad\quad \sum_nu_m^2[n]=1\quad\quad E_m=\sum_ne_m^2[n]=G^2\sum_nu_m^2[n]=G^2 em[n]=Gum[n]num2[n]=1Em=nem2[n]=G2num2[n]=G2

LPC方程求解

  • Yule-Walker Equations ∑ j = 1 p a j ϕ m [ i , j ] = ϕ m [ i , 0 ] i = 1 , 2 , . . . , p \sum_{j=1}^pa_j\phi_m[i,j]=\phi_m[i,0]\quad\quad i=1,2,…,p j=1pajϕm[i,j]=ϕm[i,0]i=1,2,...,p
    • 本质:矩阵求逆,但是存在有效解法
  • 方差矩阵法(Spoken Language Processing , Section 6.3.2.1)
    在这里插入图片描述
  • 自相关法
    • 加窗 x m [ n ] = x [ m + n ] w [ n ] x_m[n]=x[m+n]w[n] xm[n]=x[m+n]w[n]
    • 预测误差 E m = ∑ n = 0 N + p − 1 e m 2 [ n ] E_m=\sum_{n=0}^{N+p-1}e_m^2[n] Em=n=0N+p1em2[n]
    • 自相关
      ϕ m [ i , j ] = ∑ n = 0 N + p − 1 x m [ n − i ] x m [ n − j ] = ∑ n = 0 N − 1 − ( i − j ) x m [ n ] x m [ n + i − j ] \phi_m[i,j]=\sum_{n=0}^{N+p-1}x_m[n-i]x_m[n-j]=\sum_{n=0}^{N-1-(i-j)}x_m[n]x_m[n+i-j] ϕm[i,j]=n=0N+p1xm[ni]xm[nj]=n=0N1(ij)xm[n]xm[n+ij] ϕ m [ i , j ] = R m [ i − j ] \phi_m[i,j]=R_m[i-j] ϕm[i,j]=Rm[ij] R m [ k ] = ∑ n = 0 N − 1 − k x m [ n ] x m [ n + k ] R_m[k]=\sum_{n=0}^{N-1-k}x_m[n]x_m[n+k] Rm[k]=n=0N1kxm[n]xm[n+k]

语音生成模型

参考线性预测编码
语音生成模型:在这里插入图片描述

LPC正是基于这个模型的语音生成技术。在该模型中,语音信号是由一个激励信号 e ( k ) e(k) e(k)经过一个时变的全极点滤波器产生。全极点滤波器的系数取决于所产生的特定声音的声道形状。激励信号 e k e_{k} ek要么是浊音语音的脉冲序列,要么是无声声音的随机噪声。生成语音信号 s ( k ) s(k) s(k)可以表示为:
s ( k ) = ∑ p = 1 P a p s ( k − p ) + e ( k ) s(k)=\sum_{p=1}^{P}a_{p}s(k-p)+e(k) s(k)=p=1Paps(kp)+e(k)
其中, P 是滤波器的阶数, a p a_{p} ap 是滤波器的系数。LPC就是在已知 s ( k ) s(k) s(k) 的情况下获取 a p a_{p} ap .

求取 a p a_{p} ap 最常用的一个方法就是最小化真实信号与预测信号之间的均方误差(Mean Squared Error, MSE)。MSE函数可以表示为
J = E [ e 2 ( k ) ] = E [ ( s ( k ) − ∑ p = 1 P a p s ( k − p ) ) 2 ] J=E\left[e^{2}(k)\right]=E\left[\left(s(k)-\sum_{p=1}^{P}a_{p}s(k-p)\right)^{2}\right] J=E[e2(k)]=E(s(k)p=1Paps(kp))2
然后,计算 J 关于每个滤波器系数的偏导,并令其值等于0,可得(3):
∂ J ∂ a p = 0 ( 3 ) \frac{\partial J}{\partial a_{p}}=0\quad\quad\quad(3) apJ=0(3)

通过对(3)计算,可以得到(4):
∑ u = 1 P a u E [ s ( k − p ) s ( k − u ) ] = E [ s ( k ) s ( k − u ) ] ,   1 ≤ u ≤ p ( 4 ) \sum_{u=1}^{P}a_{u}E\left[s(k-p)s(k-u)\right]=E\left[s(k)s(k-u)\right],~1\leq u\leq p\quad\quad\quad(4) u=1PauE[s(kp)s(ku)]=E[s(k)s(ku)], 1up(4)

其中, 1 ≤ p ≤ P 1\leq p \leq P 1pP 。用数值 1,2,…,P 分别替换(4)中的变量 p ,我们可以得到 P 个关于滤波器系数的线性方程组,求解该线性方程组,即可得到滤波器系数的解。求解该方程组最常用高效的方法是Levinson-Durbin算法。

Matlab参考:

  1. https://wenku.baidu.com/view/25a086fc19e8b8f67c1cb938.html
  2. https://blog.csdn.net/kaixinshier/article/details/72142889

上述MSE期望也可以写作: e ( n ) = x ( n ) − x ~ ( n ) = x ( n ) − ∑ i = 1 p a i x ( n − i ) e(n)=x(n)-\widetilde{x}(n)=x(n)-\sum_{i=1}^pa_ix(n-i) e(n)=x(n)x
(n)=
x(n)i=1paix(ni)

a i a_i ai求偏导可得:
∑ n x ( n ) x ( n − j ) = ∑ i = 1 p a i ∑ n x ( n − i ) x ( n − j ) \sum_nx(n)x(n-j)=\sum_{i=1}^pa_i\sum_nx(n-i)x(n-j) nx(n)x(nj)=i=1painx(ni)x(nj) E = ∑ n [ x ( n ) ] 2 − ∑ i = 1 p a i ∑ n x ( n ) x ( n − i ) E=\sum_n[x(n)]^2-\sum_{i=1}^pa_i\sum_nx(n)x(n-i) E=n[x(n)]2i=1painx(n)x(ni)
写成自相关形式(Yule-Walker方程):
R ( j ) = − ∑ i = 1 p a i R ( j − i ) 1 ≤ j ≤ p R(j)=-\sum_{i=1}^pa_iR(j-i)\quad\quad\quad1≤j≤p R(j)=i=1paiR(ji)1jp

拆写加权式子,即为Toeplize矩阵:在这里插入图片描述
使用Durbin算法来求解Toeplize矩阵,即可计算出滤波器系数 a i a_i ai

Matlab中自带lpc函数,数学推导过程看《语音信号数字处理(L.R.Rabiner)》。

Matlab程序:

[x,fs] = audioread('1.wav');   %这里读取的双声道信号 x数据,fs采样率
sound(x,fs);    %播放音频

x1 = x(:,1);	
% figure
% plot(t,x1)
n = 200; % n is the length of a frame.  % n行   
p0 = 50; % p0 is the overlap length.    % 在前面填充p0个0。自相关法 两端都需要加P个零取样值,会造成谱估计失真
xx = buffer(x1, n, p0);  % 列:L/(n-p0),列理解为帧数

m = 8; % (m)th frame of data.
y = xx((m-1)*n+1:m*n); % select a frame of data.
p = 12; % p is the order of the AR model.	阶数
ar = lpc(y,p); % calculate the coefficients of AR model.就是前面的滤波器系数a_i
est_x = filter([0 -ar(2:end )],1,y); % calculate the predicted signal. 建立语音帧的正则方程
err = y - est_x; % calculate the residual signal.

figure
plot(x1)
title('原始信号');
figure
subplot(2,2,1); 
plot(y,'r');
title('原始一帧');
subplot(2,2,2); 
plot(est_x);
title('lpc预测的一帧');
subplot(2,2,3); 
plot(err,'r');
title('残余信号');

效果图:
在这里插入图片描述

参考《语音信号处理》实验3-LPC特征提取

I = audioread('1.wav');   %读入原始语音
I = I(:,1);
plot(I);
title('原始语音波形');%对指定帧位置进行加窗处理
Q = I';
N = 256;    %窗长
Hamm = hamming(N); %加窗
frame = 60;%需要处理的帧位置
M = Q(((frame - 1) * (N / 2) + 1):((frame - 1) * (N / 2) + N));
Frame = M.*Hamm';   %加窗后的语音帧
[B,F,T] = specgram(I,N,N/2,N);
[m,n] = size(B);
for i = 1:m 
    FTframe1(i) = B(i,frame);
end
% P = input('请输入预测器阶数?=?');
P = 5;    % 预测器阶数   改变不同阶数 观察变化
ai = lpc(Frame,P);  %计算lpc系数
LP = filter( [0 - ai(2:end)],1,Frame); %建立语音帧的正则方程
FFTlp = fft(LP);
E = Frame - LP;     % 预测误差
figure
subplot(2,1,1),plot(1:N,Frame,1:N,LP,'-r');grid;
title('原始语音和预测语音波形 ');
subplot(2,1,2)
plot(E);
grid;
title('预测误差');
% pause

fLength(1:2*N) = [M,zeros(1,N)];
Xm = fft(fLength, 2 * N);
X = Xm .* conj(Xm);
Y = fft(X , 2* N);
Rk = Y(1 : N);
PART = sum(ai(2:P+1) .* Rk(1:P));
G = sqrt(sum(Frame.^2) - PART);
A = (FTframe1 - FFTlp(1:length(F')))./FTframe1;
figure
subplot(2,1,1),plot(F',20*log(abs(FTframe1)), F',(20*log(abs(1 ./A))),'-r');
grid;
xlabel('频率/dB');ylabel('幅度');
title('短时谱');
subplot(2,1,2),plot(F',(20*log(abs(G./A))));grid;
xlabel('频率/dB');ylabel('幅度');
title('LPC谱');
% pause

%求出预测误差的倒谱
pitch = fftshift(rceps(E));
M_pitch = fftshift(rceps(Frame));
figure
subplot(2,1,1),plot(M_pitch);grid;
xlabel('语音帧');ylabel('/dB');
title('原始语音帧倒谱');
subplot(2,1,2),plot(pitch);grid;
xlabel('语音帧');ylabel('/dB');
title('预测误差倒谱');
% pause

%画出语谱图
ai1 = lpc(I,P);   %计算原始语音lpc系数
LP1 = filter([0 - ai(2:end)], 1 ,I); % 建立原始语音的正则方程
figure
subplot(2,1,1);
specgram(I,N,N/2,N);
title('原始语音语谱图');
subplot(2,1,2);
specgram(LP1,N,N/2,N);
title('预测语音语谱图')


效果图:
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/141205.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)


相关推荐

  • 零基础学Java难不难?

    零基础学Java难不难?很多同学在学Java前都会问零基础能学Java吗?Java到底难不难学?本文我就和大家唠唠这个事。有74%的人认为不难,说难学的仅占26%。

  • Fragment 的 onResume()

    Fragment 的 onResume()/***1.Fragment第一次创建时调用*2.切换程序(如点了Home键)后恢复Fragment可见时调用*3.切换fragment的hide和visible的时候可能不会调用*/@OverridepublicvoidonResume(){super.onResume();if(isAdded()&amp;&amp;!isHidden…

  • 解决:navicat for mysql连接失败[通俗易懂]

    解决:navicat for mysql连接失败[通俗易懂]1、问题描述:在navicatformysql连接mysql8.0.23时,出现如下错误。2、原因:通过百度翻译,发现是由于navicat版本的问题,出现连接失败的原因。这也就是说需要升级navicat版本。通过搜索,发现navicat是收费的,升级将会面临其他不可控的问题。于是需要寻找其他方法。通过查阅资料以及他人的经历分享。我得知了:mysql8之前的版本中加密规则是mysql_native_password,而在mysql8之后,加密规则是caching_sha2_password

    2022年10月14日
  • 如何利用Javascript发送GET/POST请求「建议收藏」

    如何利用Javascript发送GET/POST请求「建议收藏」如何利用Javascript发送GET/POST请求最近在做基于TWS的分析系统,因为采用Flask+Java的技术架构方案,所以需要开发Web,然而我自己没有做过类似的开发,所以很多工作是从头开始学着做的。因此,在实现表单数据提交的时候,当时就想到个问题,如果一个页面里内容足够多的话,仅用form提交的话,后台就需要做非常复杂的判断,以此确认用户提交的是哪类数据,这样工程不仅难看,而且低效。于是咨

  • 【python】lambda表达式与排序

    【python】lambda表达式与排序lambda表达式简单易用的匿名函数文章目录lambda表达式1.什么是lambda表达式2.lambda表达式语法3.lambda表达式的主要用途3.1list.sort()函数3.2自定义属性排序3.3常见的小问题3.4二维列表的排序1.什么是lambda表达式在学习lambda表达式之前,我们先写一个求圆的面积的函数defget_area(radius):return3.14*radius**2radius=float(input())pri

    2022年10月17日
  • linux 开启allow_url_fopen,如何开启allow_url_fopen函数[通俗易懂]

    linux 开启allow_url_fopen,如何开启allow_url_fopen函数[通俗易懂]有些程序比如dede和discuz,都会有需要打开这个函数,不打开这个函数的甚至无法安装!如何解决这个问题呢?这里给出打开这个函数的终极解决办法:1.首先确保你拥有服务器的操作权限,如果只是虚拟空间客户,那么你就联系主机商帮助操作吧;2.打开PHP的配置文件php.ini,如果你的主机是win2003该文件在C:\WINDOWS目录下,直接用记事本打开就可以,如果是apache那么就是在你的php…

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号