1 原理分析
在获取语谱图数据之前,我们需要先了解短时傅里叶变换。语音信号是典型的非平稳信号,但是由于其非平稳性由发声器官的物理运动过程而产生,这种过程是相对变换缓慢的,在10~30ms以内可以认为是平稳的。傅里叶分析时分析线性系统和平稳信号稳态特征的手段,而短时傅里叶分析,是用稳态分析方法处理非平稳信号的一种方法。
假设语音波形时域信号为
x
(
l
)
x(l)
x(l),加窗分帧处理后得到的第
n
n
n帧语音信号为
x
n
(
m
)
x_n(m)
xn(m),那有:
x
n
(
m
)
=
w
(
m
)
x
(
n
+
m
)
,
1
⩽
m
⩽
N
x_n(m)=w(m)x(n+m),1\leqslant m\leqslant N
xn(m)=w(m)x(n+m),1⩽m⩽N
对分帧信号进行短时傅里叶变化就是:
X
n
(
e
j
w
)
=
∑
m
=
1
N
x
n
(
m
)
e
−
j
w
m
X_n(e^{jw})=\sum\limits_{m=1}^Nx_n(m)e^{-jwm}
Xn(ejw)=m=1∑Nxn(m)e−jwm
其中,定义角频率
w
=
2
π
k
/
N
w=2\pi k/N
w=2πk/N,得到了离散的短时傅里叶变化(DFT)。实际上就是
X
n
(
e
j
w
)
X_n(e^{jw})
Xn(ejw)在频域的取样:
X
n
(
e
j
2
π
k
N
)
=
X
n
(
k
)
=
∑
m
=
1
N
x
n
(
m
)
e
−
j
2
π
k
m
N
,
1
⩽
k
⩽
N
X_n(e^{j\frac{2\pi k}{N}})=X_n(k)=\sum\limits_{m=1}^Nx_n(m)e^{-j\frac{2\pi km}{N}},1\leqslant k \leqslant N
Xn(ejN2πk)=Xn(k)=m=1∑Nxn(m)e−jN2πkm,1⩽k⩽N
实际中,可以使用FFT算法代替换成
x
n
(
m
)
x_n(m)
xn(m)到
X
n
(
k
)
X_n(k)
Xn(k)的转换。
def STFFT(x, win, nfft, inc):
xn = enframe(x, win, inc)
xn = xn.T
y = np.fft.fft(xn, nfft, axis=0)
return y[:nfft // 2, :]
输入数据首先分帧处理,使用分帧函数enframe(x, win, inc)
。然后直接调用np.fft.fft(xn, nfft, axis=0)
进行fft变化处理,中间有一个转置操作,是为了让时间轴作为横坐标,k作为纵坐标。
一般定义
∣
X
n
(
k
)
∣
|X_n(k)|
∣Xn(k)∣为
x
n
(
m
)
x_n(m)
xn(m)的短时幅度谱估计,而时间处频谱能量密度函数
P
(
n
,
k
)
P(n,k)
P(n,k)表示为:
P
(
n
,
k
)
=
∣
X
n
(
k
)
∣
2
P(n,k)=|X_n(k)|^2
P(n,k)=∣Xn(k)∣2
可以看出
P
(
n
,
k
)
P(n,k)
P(n,k)是一个非负的实数矩阵,以时间n作为横坐标,k作为纵坐标,就可以绘制一张热图(或灰度图),这就是语谱图。如果通过
10
lg
P
(
n
,
k
)
10\lg P(n,k)
10lgP(n,k)处理后,语谱图的单位就是dB,将变换后的矩阵精细图像和色彩映射后,就能得到彩色的语谱图。如下图所示。
语谱图中的横杠表示他们是共振峰,从横杠对应的频率和宽度可以确定相应的共振峰的频率域带宽,在一个语音段中,有没有横杠的出现是判断是不是浊音的重要标志。竖条是语谱图中与时间轴垂直的条纹,每个竖直条表示一个基音,条纹的起点相当于声门脉冲的起点,条纹之间的距离表示基音周期。
在python中,读取到语音信号以后可以直接使用plt.specgram
函数绘制出语谱图。
plt.specgram(data, NFFT=256, Fs=fs, window=np.hanning(256))
plt.ylabel('Frequency')
plt.xlabel('Time(s)')
plt.show()
2 实现代码
进行绘制语谱图,自己使用短时傅里叶变化得到的结果来做,那么首先输出的结果是一个复数矩阵,所以先求模np.abs(y)
的平方,那么用plt.pcolormesh
可以得到结果。在绘制之前需要转化为dB单位,就是以10取对数,不然啥也看不见,黑乎乎一片。完整代码如下:
CommandLineSpectrogram.py
'''
手动求求语谱图
'''
import librosa
import numpy as np
import matplotlib.pyplot as plt
#计算每帧对应的时间
def FrameTimeC(frameNum, frameLen, inc, fs):
ll = np.array([i for i in range(frameNum)])
return ((ll - 1) * inc + frameLen / 2) / fs
#分帧函数
def enframe(x, win, inc=None):
nx = len(x)
if isinstance(win, list) or isinstance(win, np.ndarray):
nwin = len(win)
nlen = nwin # 帧长=窗长
elif isinstance(win, int):
nwin = 1
nlen = win # 设置为帧长
if inc is None:
inc = nlen
nf = (nx - nlen + inc) // inc
frameout = np.zeros((nf, nlen))
indf = np.multiply(inc, np.array([i for i in range(nf)]))
for i in range(nf):
frameout[i, :] = x[indf[i]:indf[i] + nlen]
if isinstance(win, list) or isinstance(win, np.ndarray):
frameout = np.multiply(frameout, np.array(win))
return frameout
#加窗
def hanning_window(N):
nn = [i for i in range(N)]
return 0.5 * (1 - np.cos(np.multiply(nn, 2 * np.pi) / (N - 1)))
#短时傅里叶变换
def STFFT(x, win, nfft, inc):
xn = enframe(x, win, inc)
xn = xn.T
y = np.fft.fft(xn, nfft, axis=0)
return y[:nfft // 2, :]
path=r"F:\python\SimilarityAndSpectrogram\voice\bluesky.wav"#audio002.wav
data, fs = librosa.load(path, sr=None, mono=False)#sr=None声音保持原采样频率, mono=False声音保持原通道数
wlen = 256
nfft = wlen
win = hanning_window(wlen)
inc = 128
y = STFFT(data, win, nfft, inc)
FrequencyScale = [i * fs / wlen for i in range(wlen // 2)] #频率刻度
frameTime = FrameTimeC(y.shape[1], wlen, inc, fs) #每帧对应的时间
LogarithmicSpectrogramData=10*np.log10((np.abs(y)*np.abs(y))) #取对数后的数据
#np.savetxt("SpectrogramData.txt",LogarithmicSpectrogramData)
plt.pcolormesh(frameTime, FrequencyScale,LogarithmicSpectrogramData)
plt.colorbar()
#plt.savefig('语谱图22.png')
plt.show()
#掉包实现
plt.specgram(data, NFFT=256, Fs=fs, window=np.hanning(256))
plt.ylabel('Frequency')
plt.xlabel('Time(s)')
#plt.savefig('语谱图11.png')
plt.show()
3 结果显示
自己实现的图像如下:
为了验证正确性,直接使用plt.specgram
函数绘制出语谱图,结果对比是一致的,如下图所示。