以MFCC提取为例:
MFCC的计算过程
原始语言——>分帧、预加重、加窗——>FFT——>Mel滤波器组——>DTC——>MFCC
PCM音频数据
PCM(Pulse Code Modulation)是一种编码方式,可以将语音信号离散化
主要过程是将话音、图像等模拟信号每隔一定时间进行取样,使其离散化,同时将抽样值按分层单位四舍五入取整量化,同时将抽样值按一组二进制码来表示抽样脉冲的幅值。
声信号:由声源的激励脉冲与声道冲击响应卷积而成
下面通过同态处理实现声源的激励脉冲与声道冲击响应的分离。设声信号
x
(
n
)
x(n)
x(n)为
x
(
n
)
=
x
1
(
n
)
∗
x
2
(
n
)
x(n)=x_{1}(n)*x_{2}(n)
x(n)=x1(n)∗x2(n)
式中,
x
1
x_{1}
x1为声源激励,
x
2
(
n
)
x_{2}(n)
x2(n)为声道冲击响应,"*"表示卷积运算。
同态处理可实现解卷积的过程,即将卷积信号转换为加性信号的运算,具体步骤为:
(1)将时域信号
x
(
n
)
x(n)
x(n)经过傅里叶变换得到频域信号
X
(
w
)
X(w)
X(w)
F
T
[
x
(
n
)
]
=
F
T
[
x
1
(
n
)
∗
x
2
(
n
)
]
=
X
1
(
w
)
X
2
(
w
)
=
X
(
w
)
FT[x(n)]=FT[x_{1}(n)*x_{2}(n)]=X_{1}(w)X_{2}(w)=X(w)
FT[x(n)]=FT[x1(n)∗x2(n)]=X1(w)X2(w)=X(w)
由于信号在时域以卷积运算的方式存在,经过傅里叶变换到频域后即变为乘积形式,因此可利用这个性质实现线性变换。
(2)两边取对数:
l
n
X
(
w
)
=
l
n
X
1
(
w
)
+
l
n
X
2
(
w
)
=
X
^
1
(
w
)
+
X
2
^
(
w
)
=
X
(
w
)
^
lnX(w)=lnX_{1}(w)+lnX_{2}(w)=X\hat{}_{1} (w)+X_{2}\hat{}(w)=X(w)\hat{}
lnX(w)=lnX1(w)+lnX2(w)=X^1(w)+X2^(w)=X(w)^
(3)两边取傅里叶逆变换:
F
T
−
1
[
X
^
(
w
)
]
=
F
T
−
1
[
X
^
1
(
w
)
+
X
2
^
(
w
)
]
=
x
^
1
(
w
)
+
x
2
^
(
w
)
=
x
(
w
)
^
FT^{-1}[X\hat{}(w)]=FT^{-1}[X\hat{}_{1} (w)+X_{2}\hat{}(w)]=x\hat{}_{1} (w)+x_{2}\hat{}(w)=x(w)\hat{}
FT−1[X^(w)]=FT−1[X^1(w)+X2^(w)]=x^1(w)+x2^(w)=x(w)^
x ( ^ n ) x\hat(n) x(^n)称为x(n)的复倒谱
预处理
分帧
处理的一段数据就是一帧,帧的大小取决于我们使用的是多少点的FFT。这里FFT可以暂且理解为将时域信号变为频域的工具。
加窗
用一个函数或者映射的关系将每一帧数据进行包装。具体就是这一帧的每个数据都乘以一个用窗函数生成的系数。
通常窗函数选择海明窗
公式
w
(
n
)
=
{
0.54
+
0.46
c
o
s
(
2
π
n
N
−
1
)
,
0
≥
n
≥
N
−
1
0
,
其
他
w(n)=\left\{ \begin{aligned} 0.54+0.46cos(\frac{2\pi n}{N-1}),&&0\geq n \geq N-1 \\ 0,&&其他 \\ \end{aligned} \right.
w(n)=⎩⎨⎧0.54+0.46cos(N−12πn),0,0≥n≥N−1其他
FFT
用途:加速多项式乘法
多项式
A
(
x
)
=
∑
a
i
X
i
A(x)=\sum a_{i}X^{i}
A(x)=∑aiXi
单位根
w
n
=
1
w^{n}=1
wn=1
如果n=3,w可以为 1 , − 1 + s q r t ( 3 ) i 2 , − s q r t ( 3 ) i 2 1,\frac{-1+sqrt(3)i}{2},\frac{-sqrt(3)i}{2} 1,2−1+sqrt(3)i,2−sqrt(3)i
即将单位元分为三等分
单位根的性质
w 2 n 2 k = w n k w_{2n}^{2k}=w_{n}^{k} w2n2k=wnk
w n k = − w n k + n 2 w_{n}^{k}=-w_{n}^{k+\frac{n}{2}} wnk=−wnk+2n
w n k = e 2 k i π n = c o s ( 2 k π n ) + i s i n ( 2 k π n ) w_{n}^{k}=e^{\frac{2ki\pi}{n}}=cos(\frac{2k\pi}{n})+isin(\frac{2k\pi}{n}) wnk=en2kiπ=cos(n2kπ)+isin(n2kπ)
FFT的过程
系数表示法——>点值表示法——>系数表示法
A
(
x
)
=
∑
a
i
X
i
A(x)=\sum a_{i}X^{i}
A(x)=∑aiXi
1.分为奇数项和偶数项
分别记为
A
1
(
x
2
)
,
A
2
(
x
2
)
A_{1}(x^{2}),A_{2}(x^{2})
A1(x2),A2(x2)
2.将
x
x
x换为
w
n
k
w_{n}^{k}
wnk
3.所以,只需要知道
A
1
(
w
n
2
k
)
,
A
2
(
w
n
2
k
)
A_{1}(w_{\frac{n}{2}}^{k}),A_{2}(w_{\frac{n}{2}}^{k})
A1(w2nk),A2(w2nk)就可以表达这个多项式了,然后算多项式的乘法的复杂度就降低了。
例子:求 N = 2 3 = 8 N=2^{3}=8 N=23=8点FFT变换,按N=8——>N/2=4,做4点DFT(多项式乘法)
X ( k ) = X 1 ( k ) + W n k X 2 ( k ) X(k)=X_{1}(k)+W_{n}^{k}X_{2}(k) X(k)=X1(k)+WnkX2(k)
X ( k + N / 2 ) = X 1 ( k ) − W n k X 2 ( k ) X(k+N/2)=X_{1}(k)-W_{n}^{k}X_{2}(k) X(k+N/2)=X1(k)−WnkX2(k)
先将N=8点的DFT分解为2个4点DFT
可知:
时域上:
x ( 0 ) , x ( 2 ) , x ( 4 ) , x ( 6 ) x(0),x(2),x(4),x(6) x(0),x(2),x(4),x(6)为偶子序列,
x ( 1 ) , x ( 3 ) , x ( 5 ) , x ( 7 ) x(1),x(3),x(5),x(7) x(1),x(3),x(5),x(7)为奇子序列。
频域上:
X ( 0 ) X ( 3 ) X(0)~X(3) X(0) X(3),由 X ( k ) X(k) X(k)给出
X ( 4 ) X ( 7 ) X(4)~X(7) X(4) X(7),由 X ( k + N / 2 ) X(k+N/2) X(k+N/2)给出
如果直接DFT计算的计算量为:
复乘:
N
2
N^{2}
N2次=64次,复加
N
(
N
−
1
)
N(N-1)
N(N−1)次=8×7=56次
但是用
X
1
(
k
)
X_{1}(k)
X1(k)和
X
2
(
k
)
X_{2}(k)
X2(k)计算需要:
复乘:
(
N
/
2
)
2
+
(
N
/
2
)
2
(N/2)^{2}+(N/2)^{2}
(N/2)2+(N/2)2次=32次,复加
N
/
2
(
N
/
2
−
1
)
+
N
/
2
(
N
/
2
−
1
)
N/2(N/2-1)+N/2(N/2-1)
N/2(N/2−1)+N/2(N/2−1)次=12+12=24次
4.上面还不是计算量,下面介绍蝴蝶变换
原序列::::0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
位逆序置换后:0 8 4 12 2 10 6 14 1 9 5 13 3 11 7 15
回想:1的二进制是0001,8的二进制则是1000
规律:原序列和置换后的序列的对应位转换为二进制之后是刚好翻转过来的
上面的例题有4个蝶形结,每个蝶形结需要1次复乘,2次复加
所以,用分解的方法X(k)需要:
复乘:32+4=36次
复加:24+8=32次
5.因为4点DFT还是比较麻烦,所以再继续分解,这就叫做DIT-FFT
x ( 0 ) , x ( 2 ) , x ( 4 ) , x ( 6 ) x(0),x(2),x(4),x(6) x(0),x(2),x(4),x(6)为偶子序列,
x ( 1 ) , x ( 3 ) , x ( 5 ) , x ( 7 ) x(1),x(3),x(5),x(7) x(1),x(3),x(5),x(7)为奇子序列。再分
x 1 ( r ) : x ( 0 ) , x ( 4 ) x_{1}(r):x(0),x(4) x1(r):x(0),x(4)偶序列 x ( 2 ) , x ( 6 ) x(2),x(6) x(2),x(6)奇序列
x 2 ( r ) : x ( 1 ) , x ( 5 ) x_{2}(r):x(1),x(5) x2(r):x(1),x(5)偶序列 x ( 3 ) , x ( 7 ) x(3),x(7) x(3),x(7)奇序列
整个过程(DIT-FFT运算流图N=8)
在Matlab里如何进行FFT变换
[y,Fs,bits]=wavread(‘q.wav’)
y代表着语音信号的变量名;FS采样率;BITS采样位数
q.wav保存的源语音文件;wavread函数名
sound(y,Fs,bits)
有这三个变量就可以播放
n=length(y)
t=(0:n-1)/Fs
plot(t,y);xlabel(‘s’);title(‘时域波形图’)
绘制时域波形
x=fft(y)
《FFT matlab 算法 实现》
https://blog.csdn.net/qq_36666115/article/details/90045013
Y=abs(x)
W=fftshift(Y)
f=linspace(-Fs/2,Fs/2,n)
%将频域横轴改为以Hz为单位
plot(f,W);xlabel(‘频率(Hz)’);title(‘频域波形图’)
绘制频域图形
声纹时间频率(spectrum)
是将声音信号进行分帧、加窗及离散傅里叶变换之后,按照时间维度堆叠而成的能量谱图。具有声音信号的时域频域特征,以及直观的图像特征。
Mel滤波器组
在语音识别领域,由于人耳对于正常频率标度的感知是非线性的,因此常用Mel滤波将声纹时域谱变换到Mel标度下的Mel时域谱。
1.Mel时频谱可以实现人耳对频率感知的线性化处理,降低干扰频段的权重
2.Mel时频谱在保留声纹特征的前提下,大大降低样本尺寸
转换关系如下
M e l ( k ) = 2595 ∗ l g ( 1 + f / 700 ) Mel(k)=2595*lg(1+f/700) Mel(k)=2595∗lg(1+f/700)
M e l − 1 ( f ) = 700 ∗ ( 1 0 1 + k / 2595 − 1 ) Mel^{-1}(f)=700*(10^{1+k/2595}-1) Mel−1(f)=700∗(101+k/2595−1)
其中f为正常标度的频率,KaTeX parse error: Undefined control sequence: \leg at position 3: 0 \̲l̲e̲g̲ ̲f \leg 22050;k为Mel标度频率
针对样本的频率分布特点,若声纹主要频率成分为2000Hz以下的中、低频,采用Mel滤波器组对低频、中低频成分进行滤波放大处理,对中高频与高频进行滤波降权处理,从而实现样本数据的降维以及声纹特征的准确提取。
该等高三角形滤波器组的函数为:
KaTeX parse error: Undefined control sequence: \leg at position 89: …(m-1)},&&x(m-1)\̲l̲e̲g̲ ̲f \leg x(m)\\ \…
其中,m为滤波器编号,可设置为——;x(m)为三角形滤波器的中心频率,公式为
x
(
m
)
=
(
N
f
s
)
M
e
l
−
1
(
M
e
l
(
f
m
i
n
)
+
m
M
e
l
(
f
m
a
x
)
−
M
e
l
(
f
m
i
n
)
M
+
1
)
x(m)=(\frac{N}{f_{s}})Mel^{-1}(Mel(f_{min})+m\frac{Mel(f_{max})-Mel(f_{min})}{M+1})
x(m)=(fsN)Mel−1(Mel(fmin)+mM+1Mel(fmax)−Mel(fmin))
其中
f
m
a
x
,
f
m
i
n
g
f_{max},f_{ming}
fmax,fming分别为滤波器滤波范围的频率最大值和最小值,
f
s
f_{s}
fs为采样频率,取。N为离散傅里叶时的帧长度。
通过将时域谱矩阵与Mel滤波器组矩阵相乘,
再取对数,
l
o
g
X
[
k
]
=
l
o
g
H
[
k
]
+
l
o
g
E
[
k
]
log X[k]=log H[k]+log E[k]
logX[k]=logH[k]+logE[k]
得到Mel时谱图
以上是Mel时谱图的求取过程
DCT(离散余弦变换)
这一步的目的是将频域变回时域,然后取最低值(离散序列的包络)得到MFCC系数。