1 原理分析

在获取语谱图数据之前，我们需要先了解短时傅里叶变换。语音信号是典型的非平稳信号，但是由于其非平稳性由发声器官的物理运动过程而产生，这种过程是相对变换缓慢的，在10~30ms以内可以认为是平稳的。傅里叶分析时分析线性系统和平稳信号稳态特征的手段，而短时傅里叶分析，是用稳态分析方法处理非平稳信号的一种方法。

假设语音波形时域信号为x ( l ) x(l)x(l),加窗分帧处理后得到的第n nn帧语音信号为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 Nxn​(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∑N​xn​(m)e−jwm

其中，定义角频率w = 2 π k / N w=2\pi k/Nw=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 NXn​(ejN2πk​)=Xn​(k)=m=1∑N​xn​(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, :]12345

输入数据首先分帧处理，使用分帧函数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)|^2P(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()1234

进行绘制语谱图，自己使用短时傅里叶变化得到的结果来做，那么首先输出的结果是一个复数矩阵，所以先求模np.abs(y)的平方，那么用plt.pcolormesh可以得到结果。在绘制之前需要转化为dB单位，就是以10取对数，不然啥也看不见，黑乎乎一片。完整代码如下：

CommandLineSpectrogram.py1

'''手动求求语谱图'''import librosaimport numpy as npimport 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.wavdata, fs = librosa.load(path, sr=None, mono=False)#sr=None声音保持原采样频率， mono=False声音保持原通道数wlen = 256nfft = wlenwin = hanning_window(wlen)inc = 128y = 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()1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071

自己实现的图像如下：

为了验证正确性，直接使用plt.specgram函数绘制出语谱图，结果对比是一致的，如下图所示。

感谢阅读！！！

多说一句，很多人学Python过程中会遇到各种烦恼问题，没有人解答容易放弃。小编是一名python开发工程师，这里有我自己整理了一套最新的python系统学习教程，包括从基础的python脚本到web开发、爬虫、数据分析、数据可视化、机器学习等。想要这些资料的可以关注小编，并在后台私信小编：“01”即可领取。