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

2 实现代码

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

CommandLineSpectrogram.py
1

'''
手动求求语谱图
'''
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()
1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465666768697071

3 结果显示

自己实现的图像如下：

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

感谢阅读！！！

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