python绘制语谱图(手动实现)

1 原理分析

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

假设语音波形时域信号为$x(l)$,加窗分帧处理后得到的第$n$帧语音信号为$x_n(m)$,那有：
$$x_n(m)=w(m)x(n+m),1⩽m⩽N$$
对分帧信号进行短时傅里叶变化就是：
$$X_n(e^{jw})=\sum _{m=1}^N{x}_n(m)e^{-jwm}$$
其中，定义角频率$w=2\pi k/N$，得到了离散的短时傅里叶变化(DFT)。实际上就是$X_n(e^{jw})$在频域的取样：
$$X_n(e^{j\frac{2\pi k}{N}})=X_n(k)=\sum _{m=1}^N{x}_n(m)e^{-j\frac{2\pi km}{N}},1⩽k⩽N$$
实际中，可以使用FFT算法代替换成$x_n(m)$到$X_n(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(m)$的短时幅度谱估计，而时间处频谱能量密度函数$P(n,k)$表示为：
$$P(n,k)=|X_n(k){|}^2$$
可以看出$P(n,k)$是一个非负的实数矩阵，以时间n作为横坐标，k作为纵坐标，就可以绘制一张热图（或灰度图），这就是语谱图。如果通过$10\lg P(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函数绘制出语谱图，结果对比是一致的，如下图所示。