语音信号线性预测分析

el/2024/7/24 2:47:49

语音信号线性预测分析

基本思想:一个语音取样的现在值可以用若干个语音取样过去值的加权线性组合来逼近(最小均方误差)。线性预测最重要的优势在于它可以较为精确地估计语音的参数,而这些极少的参数可以正确地表现语音信号的时域和频域特性。

基本原理:线性预测分析的基本原理是把信号用一个模型来表示,即将信号看作某一个模型的输出,这样就可以用该模型的参数来描述信号。

数学公式:
x(n)=∑i=1paix(n−i)+e(n)x(n)=\sum_{i=1}^{p} a_ix(n-i)+e(n) x(n)=i=1paix(ni)+e(n)
x(n)为真实信号,加权项为预测信号,e(n)为预测误差。根据e(n)的最小均方误差准则来计算滤波器系数ai。
e(n)=x(n)−∑i=1paix(n−i)e(n)=x(n)-\sum_{i=1}^{p} a_ix(n-i) e(n)=x(n)i=1paix(ni)

E=∑ne2(n)=∑n[x(n)−∑i=1paix(n−i)]2E=\sum_{n}e^{2}(n)=\sum_{n}[x(n)-\sum_{i=1}^{p}a_ix(n-i)]^{2} E=ne2(n)=n[x(n)i=1paix(ni)]2

为使预测的最小均方误差E最小,对ai求偏导等0:
∂E∂aj=0(1≤j≤p)\frac{\partial E}{\partial a_j} = 0(1\le j\le p) ajE=0(1jp)
最小误差结果为:
∑nx(n)x(n−j)=∑i=1pai∑nx(n−i)x(n−j)\sum_{n}x(n)x(n-j)=\sum_{i=1}^{p}a_i\sum_{n}x(n-i)x(n-j) nx(n)x(nj)=i=1painx(ni)x(nj)

E=∑n[x(n)]2−∑i=1pai∑nx(n)x(n−i)E=\sum_{n}[x(n)]^{2} -\sum_{i=1}^{p}a_i\sum_{n}x(n)x(n-i) E=n[x(n)]2i=1painx(n)x(ni)

通过自相关解法来求解线性预测系数(Yule-Walker方程),根据加窗自相关函数定义,以及最小误差结果可得:
R(j)=−∑i=1paiR(j−i)1≤j≤pR(j)=-\sum_{i=1}^{p}a_iR(j-i) \quad 1\le j\le p R(j)=i=1paiR(ji)1jp
最小均方误差可写为:
E=R(0)+∑i=1paiR(i)E=R(0)+\sum_{i=1}^{p}a_iR(i) E=R(0)+i=1paiR(i)
利用Toeplitz矩阵的性质进行计算,基本思想是递归方法求解,采用Levinson-Durbin算法。

算法过程与步骤:

当i=0时,E0=R(0),a0=1;

对于第i次递归(i=1,2,…,p):
ki=1Ei−1[R(i)−∑j=1i−1aji−1R(j−i)]k_i=\frac{1}{E_{i-1}}\left [ R(i)-\sum_{j=1}^{i-1}a_{j}^{i-1}R(j-i)\right ] ki=Ei11[R(i)j=1i1aji1R(ji)]

iai(i)=kiia_i^{(i)}=k_i iai(i)=ki

aj(i)=aj(i−1)−kiai−j(i−1)a_j^{(i)}=a_j^{(i-1)}-k_ia_{i-j}^{(i-1)} aj(i)=aj(i1)kiaij(i1)

Ei=(1−ki2)Ei−1E_i=(1-k_i^{2})E_{i-1} Ei=(1ki2)Ei1

最终解为:
ai=aj(p)1≤j≤pa_i=a_j^{(p)}\quad 1\le j\le p ai=aj(p)1jp

Ep=R(0)∏i=1p(1−ki2)E_p=R(0)\prod_{i=1}^{p}(1-k_i^2) Ep=R(0)i=1p(1ki2)

由上式可知,最小均方误差Ep一定大于0,且随着预测器阶数的增加而减少。参数ki一定满足
∣ki∣<1,1≤i≤p\left | k_i \right | <1,\quad 1\le i\le p ki<1,1ip
MATLAB工具箱中有lpc函数,使用Levinson-Durbin的自相关方法计算预测系数的。输入一帧的数据以及线性预测阶数,输出预测系数ar和预测误差的方差r。
对一帧语音信号进行实验,预测结果以及误差:
image-20210728164817790
根据Levinson-Durbin自相关法求线性预测系数的原理,编写函数,并于lpc函数进行比较,结果一致。
image-20210729141140143
一帧信号的FFT频谱与LPC预测系数的复频谱比较:
image-20210729145826050

倒谱分析

倒谱(Cepstrum)是表示一帧语音数据特征的一个序列。可以看作和自相关序列类似的东西,也可以理解为自相关序列的对数压缩,因为其携带有和自相关序列类似的信息。信号的倒谱是通过对信号进行离散傅里叶变换,绝对值取平方得到信号的功率谱,然后取对数再做逆离散傅里叶变换得到信号的倒谱。

通过对信号x(n)取离散傅里叶变换,可以得到信号的复数谱:
DFT(x(n))→X(k)DFT(x(n))\to X(k) DFT(x(n))X(k)
通过对复数谱的绝对值取平方得到信号的功率谱:
∣X(k)∣2→P(k)\left | X(k) \right | ^{2}\to P(k) X(k)2P(k)
对功率谱取对数,再做逆离散傅里叶变换
IDFT(log(P(k)))→C(n)IDFT(log(P(k)))\to C(n) IDFT(log(P(k)))C(n)
一段话一帧倒谱实验:
image-20210729184028934

声母a一帧的倒谱实验:

image-20210729191234736

离散余弦变换

离散余弦变换DCT具有信号谱分量丰富,能量集中的特点。对于那些不重要的频域区域和系数就能够直接裁剪掉,因此,DCT变换非常适合于压缩算法的处理。

数学公式:设x(n)是N个有限值的一维实数信号序列,n=0,1,…,N-1,C(k)是正交因子。
X(k)=2N∑n=0N−1C(k)x(n)cos((2n+1)kπ2N)k=0,1,⋯,N−1X(k)=\sqrt{\frac{2}{N} }\sum_{n=0}^{N-1}C(k)x(n)cos\left ( \frac{(2n+1)k\pi }{2N} \right ) \quad k=0,1,\cdots ,N-1 X(k)=N2n=0N1C(k)x(n)cos(2N(2n+1)kπ)k=0,1,,N1

C(k)={2/2,k=01,k∈[1,N−1]C(k)=\begin{cases} \sqrt{2}/2, \quad k=0 \\ 1,\quad \quad \quad k\in [1,N-1] \end{cases} C(k)={2/2,k=01,k[1,N1]

调用DCT函数求序列x(n)的DCT系数,然后仅用幅值大于5的系数进行信号重建,比较差异。
x(n)=cos⁡(2πfn/fs)n∈[0,1000],fs=1000Hz,f=50Hzx(n)=\cos (2\pi fn/f_s)\quad n\in[0,1000],f_s=1000Hz,f=50Hz x(n)=cos(2πfn/fs)n[0,1000],fs=1000Hz,f=50Hz
实验结果:

image-20210730093117344

代码:

%lpc预测结果实验
clear
clc
[x,fs] = audioread('CHMMWJQ.wav');
sound(x,fs);
n = 200; 
xx = enframe(x, n);
m = 60;
y = x((m-1)*n+1:m*n);
p = 12;
ar = lpc(y,p);
est_x = filter([0 -ar(2:end )],1,y);
err = y - est_x;
subplot(411);plot(x);title('原始信号');xlabel('点数');ylabel('幅值');
subplot(412);plot(y);title('一帧数据');xlabel('点数');ylabel('幅值');
subplot(413);plot(est_x);title('预测结果');xlabel('点数');ylabel('幅值');
subplot(414);plot(err);title('差值');xlabel('点数');ylabel('幅值');
%用自相关法求信号s使均方预测误差为最小的预测系数的函数
function [ar,G]=lpc_coeff(s,p)
%算法为Durbin快速递推算法
n=length(s);                                    % 获得信号长度                    
for i=1:pRp(i)=sum(s(i+1:n).*s(1:n-i)); % 计算自相关函数
end
Rp_0=s'*s;      % 即Rn(0)
Ep=zeros(p,1);  % Ep为p阶最佳线性预测反滤波能量
k=zeros(p,1);   % k为偏相关系数
a=zeros(p,p);   % 以上为初始化
%i=1的情况需要特殊处理,也是对p=1进行处理
Ep_0=Rp_0;
k(1)=Rp(1)/Rp_0;
a(1,1)=k(1);
Ep(1)=(1-k(1)^2)*Ep_0;
%i=2起使用递归算法
if p>1for i=2:pk(i)=(Rp(i)-sum( a(1:i-1,i-1).*Rp(i-1:-1:1)'))/Ep(i-1);a(i,i)=k(i);Ep(i)=(1-k(i)^2)*Ep(i-1);for j=1:i-1a(j,i)=a(j,i-1)-k(i)*a(i-j,i-1);endend
end
ar=a(:,p);
ar=[1 -1*ar'];
G=sqrt(Ep(p));
%线性预测系数对比
clear all; clc; close all;
[x,fs]=wavread('CHMMWJQ.wav');
L=200;
y=x(8001:8000+L);
p=12; 
ar1=lpc(y,p);
ar2=lpc_coeff(y,p); 
est_x1=filter([0 -ar1(2:end)],1,y); 
est_x2=filter([0 -ar2(2:end)],1,y); 
err1=y-est_x1; 
err2=y-est_x2; 
subplot 321; plot(x); axis tight;title('原始信号'); ylabel('幅值')
subplot 322; plot(y); xlim([0 L]); title('一帧数据'); ylabel('幅值')
subplot 323; plot(est_x1); xlim([0 L]); title('LPC预测值'); ylabel('幅值')
subplot 324; plot(est_x2,'r'); xlim([0 L]); title('lpc\_coeff预测值'); ylabel('幅值')
subplot 325; plot(err1); xlim([0 L]); title('LPC预测误差'); ylabel('幅值'); xlabel('样点')
subplot 326; plot(err2,'r'); xlim([0 L]); title('lpc\_coeff预测误差'); ylabel('幅值'); xlabel('样点')
%比较LPC预测系数的复频谱与FFT频谱
clear all; clc; close all;
[x,fs]=wavread('CHMMWJQ.wav');
L=200;     
p=12; 
y=x(8001:8000+L); 
ar=lpc(y,p); 
nfft=512;
W2=nfft/2;
m=1:W2+1;
Y=fft(y,nfft); 
ff=(fft(ar.',2*W2+2).').^(-1);
Y1=ff(1:length(ff)/2); 
% 作图
subplot 211; plot(y);title('一帧语音信号的波形'); ylabel('幅值'); xlabel('(a)')
subplot 212; plot(m,20*log10(abs(Y(m))));axis([0 200 -inf inf]);
hold on;
plot(m,20*log10(abs(Y1)),'color','r');axis([0 200 -inf inf]); 
ylabel('幅值/db');legend('FFT频谱','LPC复频谱'); xlabel('样点');title('FFT频谱和LPC谱的比较');
%倒谱计算
clear
clc
[x,fs]=audioread('CHMMWJQ.wav');
len=200;
% nfft=128;
y=x(6001:6000+len);
Xk=fft(y);
Pk=(abs(Xk));
Cn=ifft(log(Pk));
[xh,yh]=rceps(y);
subplot(211);plot(y);xlabel('点数');ylabel('幅值');title('一帧信号');
subplot(212);plot(Cn);xlabel('点数');ylabel('幅值');title('倒谱序列');
%DCT变换
clear
clc
f=50;
fs=1000;
N=1000;
n=(0:N-1);
xn=cos(2*pi*f*n/fs);
xd=dct(xn);
zn=xd;
num=find(abs(xd)<5);
xd(num)=0;
y=idct(xd);
subplot 211;plot(n,xn);title('原始信号');xlabel('样点');ylabel('幅值');
subplot 212;plot(n,y,'r');title('重建信号');xlabel('样点');ylabel('幅值');
%DCT变换
clear
clc
f=50;
fs=1000;
N=1000;
n=(0:N-1);
xn=cos(2*pi*f*n/fs);
xd=dct(xn);
zn=xd;
num=find(abs(xd)<5);
xd(num)=0;
y=idct(xd);
subplot 211;plot(n,xn);title('原始信号');xlabel('样点');ylabel('幅值');
subplot 212;plot(n,y,'r');title('重建信号');xlabel('样点');ylabel('幅值');

http://www.ngui.cc/el/3419325.html

相关文章

在vue中使用fastclick解决移动端300ms延时问题

第一步&#xff1a;安装插件 把fastclick这个包安装到项目的依赖之中&#xff0c;--save表示开发与上线都需要 npm install fastclick --save 第二步:在main.js中引入插件 import fastclick from fastclick 第三步:在body元素上使用插件 fastclick.attach(document.body)

如何实现洗牌算法?

day09 题目描述: 开发一款扑克游戏&#xff0c;需编写一套洗牌算法&#xff0c;公平的洗牌是将洗好的牌存储在一个整型数组里&#xff0c;每张牌被放在任何一个位置的概率是相等的. 解析: 定义一个洗牌函数&#xff0c;函数内用tmp数组存储1~54表示54张牌&#xff0c;然后对5…

球的反弹高度有多高?

day11 题目描述: 一球从100米高度自由落下&#xff0c;每次落地后反弹回原高度的一半&#xff0c;再落下。求它在第十次落地时&#xff0c;共经过多少米&#xff1f;第十次反弹多高&#xff1f; 解析: 设初始总高度为100米&#xff0c;球每次下落高度反弹回的高度为上一次的一…

如何找出1000以内的“完数“

day12 题目描述: 如果一个数恰好等于它的因子之和&#xff0c;这个数就称为"完数"&#xff0c;例如6123.编程找出1000以内的所有完数。 解析&#xff1a; 外层循环1000次&#xff0c;每次循环得到的i传入下个循环内&#xff0c;内部循环求解出符合i整除k等于0的数。…

如何在本地创建分支并推送到码云上

第一步: git branch 查看当前所在的分支&#xff0c;确保在master分支上 第二步: git checkout -b example 在当前分支上(也就是master分支)创建一个example分支&#xff0c;并切换到该分支上 第三步: git push -u or…

如何获取规定的排列组合?

day13 题目描述: 将一组数字&#xff0c;字母或符号进行排列&#xff0c;以得到不同的组合顺序&#xff0c;例如1&#xff0c;2&#xff0c;3这三个数的排列组合有:123,132,213,231,312,321. 解析&#xff1a; 可以使用递归将问题切割为较小的单元进行排列组合&#xff0c;如…

快速排序--C语言

一、快速排序算法&#xff08;Quicksort&#xff09; 1:定义&#xff1a; 快速排序由C. A. R. Hoare在1962年提出。快速排序是对冒泡排序的一种改进&#xff0c;采用了一种分治的策略。 2&#xff1a;基本思想&#xff1a; 通过一趟排序将要排序的数据分割成独立的两部分&#…

浅谈Java多线程基础

一: 线程和进程 进程: 进程是操作系统结构的基础&#xff1b;是一次程序的执行&#xff1b;是一个程序及其数据在处理机上顺序执行时所发生的活动。操作系统中&#xff0c;几乎所有运行中的任务对应一条进程,进程具有并发性&#xff0c;它可以同其他进程一起并发执行&#xff…

HTML+CSS基础----盒模型

在谈盒模型之前,我们先了解下什么是块级元素,什么是行内元素. 块级元素: 如<p>,<li>,<h1>等等; 默认情况下&#xff0c;块级元素会另起一行&#xff0c;并尽可能的充满整个容器。 块级元素可以包含行内元素和其他块级元素&#xff0c;相比于行内元素可以创…

浅谈css几种常用的选择器

首先我们得了解,什么是选择器? 每一条css样式定义由两部分组成&#xff0c;形式如下&#xff1a; [code] 选择器{样式} [/code] 在{}之前的部分就是“选择器”。 “选择器”指明了{}中的“样式”的作用对象&#xff0c;也就是“样式”作用于网页中的哪些元素.要使用css对HTML页…