基於深度學習的單通道語音增強

凌逆戰發表於2020-12-25

本文程式碼請見：https://github.com/Ryuk17/SpeechAlgorithms

部落格地址（轉載請指明出處）：https://www.cnblogs.com/LXP-Never/p/14142108.html

如果你覺得寫得還不錯，點贊?，關注是對我最大的支援，謝謝?

　　傳統的語音增強方法基於一些設定好的先驗假設，但是這些先驗假設存在一定的不合理之處。此外傳統語音增強依賴於引數的設定，人工經驗等。隨著深度學習的發展，越來越多的人開始注意使用深度學習來解決語音增強問題。由於單通道使用場景較多，本文就以單通道語音增強為例。

　　目前基於深度神經網路單通道語音增強方法大致可以分為兩類：

第一種基於對映的方法
第二種基於mask的方法

基於對映的語音增強

　　基於對映的語音增強方法通過訓練神經網路模型直接學習帶噪語音和純淨語音之間的對映關係，有兩種對映方案：

頻譜對映：使用模型預測語音的時頻域表示，之後再將語音的時頻域表示通過波形合成技術恢復成時域訊號。使用最多的時頻域表示特徵是短時傅立葉變換譜，利用人耳對相位不敏感的特性，一般只預測短時傅立葉變換幅度譜，並使用混合語音的相位合成預測語音的波形。

波形對映：直接將帶噪語音波形輸入到模型，模型直接輸出純淨語音波形的方法。

　　我們以頻譜對映舉例說明一下：

訓練階段

輸入：這裡採用較為簡單地特徵，即帶噪聲語音訊號的幅度譜，也可以採用其他的特徵。值得一提的是，如果你的輸入是一幀，對應輸出也是一幀的話效果一般不會很好。因此一般採用擴幀的技術，如下圖所示，即每次輸入除了當前幀外還需要輸入當前幀的前幾幀和後幾幀。這是因為語音具有短時相關性，對輸入幾幀是為了更好的學習這種相關性

Label：資料的label為純淨語音訊號的幅度譜，這裡只需要一幀就夠了。
損失函式：學習噪聲幅度譜與純淨語音訊號的幅度譜類似於一個迴歸問題，因此損失函式採用迴歸常用的損失函式，如均方誤差（MSE）、均方根誤差（RMSE）或平均絕對值誤差（MAE）等....
最後一層的啟用函式：由於是迴歸問題，最後一層採用線性啟用函式
其他：輸入的幅度譜進行歸一化可以加速學習過程和更好的收斂。如果不採用幅度譜可以採用功率譜，要注意的是功率譜如果採用的單位是dB，需要對資料進行預處理，因為log的定義域不能為0，簡單的方法就是在取對數前給等於0的功率譜加上一個非常小的數

增強階段

輸入：輸入為噪聲訊號的幅度譜，這裡同樣需要擴幀。對輸入資料進行處理可以在語音訊號加上值為0的語音幀，或者捨棄首尾的幾幀。如果訓練過程對輸入進行了歸一化，那麼這裡同樣需要進行歸一化
輸出：輸入為估計的純淨語音幅度譜
重構波形：在計算輸入訊號幅度譜的時候需要儲存每一幀的相位資訊，然後用儲存好的相位資訊和模型輸出的幅度譜重構語音波形，程式碼如下所示。

spectrum = magnitude * np.exp(1.0j * phase)

基於Mask的語音增強

Mask這個單詞有的地方翻譯成掩蔽有的地方翻譯成掩膜，我個人傾向於翻譯成“掩蔽”，本文就用掩蔽作為Mask的翻譯。

時頻掩蔽

　　我們都知道語音訊號可以通過時域波形或者頻域的各種頻譜表示，此外語譜圖可以同時展示時域和頻域的資訊，因此被廣泛應用，如下圖所示。語譜圖上的畫素點就可以稱為 時頻單元。

現在我們假設有兩段語音訊號，一段是純淨訊號，另一段是噪聲，他們混合在一起了，時域波形和對應的語譜圖分別如下圖所示：

如果我們想將純淨語音訊號從混合訊號中抽離在時域方面是很難做到的。現在我們從語譜圖（語音的時頻單元）角度入手去解決語音分離問題。首先我們提出兩個假設：

1、我們假設訊號能量稀疏的，即對於大多數時頻區域它的能量為0，如下圖所示，我們可以看到大多數區域的值，即頻域能量為0。

2、我們假設訊號能量不相交的，即它們的時頻區域不重疊或者重疊較少，如下圖所示，我們可以看到時頻區域不為0的地方不重疊或者有較少部分的重疊。

基於以上兩點假設，我們就可以分離我們想要的訊號和噪聲訊號。給可能屬於一個訊號源的區域分配掩碼為1，其餘的分配掩碼0，如下圖所示。

　　我們通過0和1的二值掩碼然後乘以混合訊號的語譜圖就可以得到我們想要喜好的語譜圖了，如下圖所示。

　　神經模型一般直接預測時頻掩蔽$M(t,f)$，之後再通過$M(t,f)$與混合語音$Y(t,f)$相乘得到預測的純淨語音$\hat{S}(t,f)=\hat{M}(t,f)\otimesY(t,y)$，其中$\otimes$代表哈達瑪乘積(Hadamard Product)。在語音增強研究的發展過程中，研究人員提出了一系列的時頻掩蔽作為訓練目標：

理想二值掩蔽（Ideal Binary Mask，IBM）

原理：由於語音在時頻域上是稀疏分佈的，對於一個具體的時頻單元，語音和噪聲的能量差異通常比較大，因此大多數時頻單元上的訊雜比極高或極低。IBM 是對這種現實情況的簡化描述，將連續的時頻單元訊雜比離散化為兩種狀態 1 和0，在一個時頻單元內：如果語音佔主導（高訊雜比），則被標記為 1；反之如果噪聲佔主導（低訊雜比），則標記為 0。最後將 IBM 和帶噪語音相乘，實際上就是將低訊雜比的時頻單元置零，以此達到消除噪聲的目的。

　　因此，IBM 的值由時頻單元上的訊雜比SNR(t,f)和設定的閾值比較之後決定：

$$公式1：I B M(t, f)=\left\{\begin{array}{l}
1, \operatorname{SNR}(t, f)>L C \\
0, \text { else }
\end{array}\right.$$

其中LC為閾值，一般取0，SNR計算公式為：

$$公式2：\operatorname{SNR}(t, f)=10 * \log 10\left(\frac{|S(t, f)|^{2}}{|N(t, f)|^{2}}\right)$$

優點：IBM 作為二值目標，只需要使用簡單的二分類模型進行預測，並且可以有效地提高語音的可懂度。
缺點：IBM 只有 0 和 1 兩種取值，對帶噪語音的處理過於粗暴，處理過程中引入了較大的噪聲，無法有效地改善語音質量。

我看到過很多種寫法

def IBM(clean_speech, noise):
    """計算 ideal binary mask (IBM)
    Erdogan, Hakan, et al. "Phase-sensitive and recognition-boosted speech separation using deep recurrent neural networks." ICASSP, 2015.
    :param clean_speech: 純淨語音 STFT
    :param noise: 噪聲 STFT
    :return: 純淨語音的理想二值掩膜 IBM
    """
    mask = np.zeros(np.shape(clean_speech), dtype=np.float32)
    mask[np.abs(clean_speech) >= np.abs(noise)] = 1.0
    return mask

第二種

def IBM_SNR(clean_speech, noise_speech):
    """計算 ideal binary mask (IBM)
    Erdogan, Hakan, et al. "Phase-sensitive and recognition-boosted speech separation using deep recurrent neural networks." ICASSP, 2015.
    :param clean_speech: 純淨語音 STFT
    :param noise_speech: 帶噪語音 STFT
    :return: 純淨語音的理想二值掩膜 IBM
    """
    _eps = np.finfo(np.float).eps  # 避免除以0
    theta = 0.5  # a majority vote
    alpha = 1  # ratio of magnitudes
    mask = np.divide(np.abs(clean_speech) ** alpha, (_eps + np.abs(noise_speech) ** alpha))
    mask[np.where(mask >= theta)] = 1
    mask[np.where(mask < theta)] = 0
    return mask

第二種

第三種

def IBM_SNR(clean_speech, noise_speech,delta_size):
    """計算 ideal binary mask (IBM)
    Erdogan, Hakan, et al. "Phase-sensitive and recognition-boosted speech separation using deep recurrent neural networks." ICASSP, 2015.
    :param clean_speech: 純淨語音 STFT
    :param noise_speech: 帶噪語音 STFT
    :return: 純淨語音的理想二值掩膜 IBM
    """
    _eps = np.finfo(np.float).eps  # 避免除以0
    local_snr = 0
    ibm  = np.where(10. * np.log10(np.abs(clean_speech) ** 2 / np.abs(noise_speech) ** 2) >= local_snr, 1., 0.)

    if delta_size > 0:
        ibm = ibm[:, delta_size: -delta_size]
    return ibm

View Code

理想浮值掩蔽（Ideal Ratio Mask，IRM）

原理：基於語音和噪聲正交，即不相關的假設下，即$S(t,f) ⋅ N(t,f) = 0$，IRM直接刻畫了時頻單元內純淨語音能量和帶噪語音能量的比值，是目前使用非常廣泛的一種掩蔽方法。

　　在這個假設下帶噪語音的能量可以表示為：

$$公式2：|\boldsymbol{Y}(t, f)|^{2}=|\boldsymbol{S}(t, f)+\boldsymbol{N}(t, f)|^{2}=|\boldsymbol{S}(t, f)|^{2}+|\boldsymbol{N}(t, f)|^{2}$$

因此得到 IRM 為：

$$公式3：I R M(t, f)=\left(\frac{|S(t, f)|^{2}}{|Y(t, f)|^{2}}\right)^{\beta} =\left(\frac{|S(t, f)|^{2}}{|S(t, f)|^{2}+|N(t, f)|^{2}}\right)^{\beta}$$

其中$\beta$IRM取值在 0 到 1 之間，值越大代表時頻單元內語音佔的比重越高。另外，IRM 的平方形式就是經典的維納濾波器（Wiener Filter），它是均方誤差意義上的最優濾波器。

優點：IRM 是分佈在 0 到 1 的連續值，因此 IRM 是對目標語音更加準確的刻畫，這使得 IRM 可以有效地同時提升語音的質量和可懂度。
缺點：使用未處理的相位資訊進行語音重構（相位對於感知質量也很重要）

def IRM(clean_speech, noise):
    """計算Compute ideal ratio mask (IRM)
    "Phase-sensitive and recognition-boosted speech separation using deep recurrent neural networks," in ICASSP 2015, Brisbane, April, 2015.
    :param clean_speech: 純淨語音 STFT
    :param noise: 噪音 STFT
    :return: 在原始音訊域中分離(恢復)的語音
    """
    _eps = np.finfo(np.float).eps  # 防止分母出現0
    mask = np.abs(clean_speech) / (np.abs(clean_speech) + np.abs(noise) + _eps)
    return mask


def Wiener_like(clean_speech, noise):
    """計算Wiener-like Mask
    "Phase-sensitive and recognition-boosted speech separation using deep recurrent neural networks," in ICASSP 2015, Brisbane, April, 2015.
    :param clean_speech: 純淨語音 STFT
    :param noise: 噪音 STFT
    :return: 在原始音訊域中分離(恢復)的語音
    """
    _eps = np.finfo(np.float).eps  # 防止分母出現0
    mask = np.divide((np.abs(clean_speech) ** 2 + _eps),
                     (np.abs(clean_speech) ** 2 + np.abs(noise) ** 2) + _eps)
    return mask

理想幅度掩蔽（Ideal Amplitude Mask，IAM）

原理：IAM也稱為Spectral Magnitude Mask（SMM），不對噪聲和語音做出正交假設，IAM刻畫的也是純淨語音和帶噪語音的能量比值

$$公式4：\operatorname{IAM}(t, f)=\frac{|S(t, f)|}{|Y(t, f)|}$$

由於在語音和噪聲疊加的過程中，存在反相相消的情況，因此並不能保證帶噪語音的幅值總是大於純淨語音的幅值，因此 IAM 的範圍是$[0,+\infty ]$。如果目標中出現非常大的數值，會導致訓練過程出現異常。為了穩定訓練，一般會將 IAM 進行截斷到一定的範圍內。為了確定合適的截斷範圍，我們可以在訓練資料上取樣 100 句語音並計算 IAM，就可以對IAM 的數值範圍得到一個近似的估計，得到如圖 3.4 的結果。一般將 IAM 截斷到[0, 1]或者[0, 2]即可，因為只有非常少部分的 IAM 落在了$[2,+\infty ]$的區間內。

圖* IAM數值分佈直方圖

def IAM(clean_speech, noise_speech):
    """計算ideal amplitude mask (IAM)
    "Phase-sensitive and recognition-boosted speech separation using deep recurrent neural networks," in ICASSP 2015, Brisbane, April, 2015.
    :param clean_speech: 純淨語音 STFT
    :param noise_speech: 帶噪語音 STFT
    :return:
    """
    _eps = np.finfo(np.float).eps  # 避免除以0
    mask = np.abs(clean_speech) / (np.abs(noise_speech) + _eps)
    return mask

相位敏感掩蔽（Phase Sensitive Mask，PSM）

原理：PSM考慮到相位誤差的時頻掩蔽

　　PSM在形式上是 IAM 乘上純淨語音和帶噪語音之間的餘弦相似度

$$公式5：P S M(t, f)=\frac{|S(t, f)|}{|Y(t, f)|} \cos \left(\theta^{S}-\theta^{Y}\right)$$

式中$\theta^{S}-\theta^{Y}$表示純淨語音和帶噪語音的相位差，不難看出，PSM 的取值範圍是$[-\infty,+\infty]$，因此也需要截斷，我們同樣使用直方圖統計PSM的數值分佈範圍，從下圖可以看出在0 和 1 附近出現兩個明顯的峰值，這也再次說明了 IBM 目標設計的合理性。為了方便，一般將 PSM 截斷到[0, 1]，或者是適當將截斷的區間放大到[-1, 2]。

PSM數值分佈直方圖

優點：純淨語音相位和帶噪語音相位的差異，加入相位資訊之後，PSM方法可以獲得更高的SNR，因而降噪效果比IAM更好。

def PSM(clean_speech, noise_speech):
    """計算ideal phase-sensitive mask (PSM)
    :param clean_speech: 純淨語音 STFT
    :param noise_speech:帶噪語音 STFT
    :return:
    """
    _eps = np.finfo(np.float).eps  # 防止分母出現0
    clean_speech_phase = np.angle(clean_speech)
    noise_speech_phase = np.angle(noise_speech)
    mask = np.abs(clean_speech) / np.abs(noise_speech) * np.cos(clean_speech_phase - noise_speech_phase)
    # Truncated Phase Sensitive Masking
    # Theta = np.clip(np.cos(clean_speech_phase-noise_speech_phase), a_min=0., a_max=1.)
    # mask = np.divide(np.abs(clean_speech), _eps + np.abs(noise_speech)) * Theta
    return mask

複數理想浮值掩蔽（Complex Ideal Ratio Mask，cIRM）

參考文獻：2015_Complex ratio masking for monaural speech separation

原理：在複數域的理想浮值掩膜，同時增強幅度譜和相位譜

$條件：\left\{ \begin{array}{l}Y = {Y_r} + i{Y_i}\\M = {M_r} + i{M_i}\\S = {S_r} + i{S_i}\\{S_{t,f}} = {M_{t,f}}*{Y_{t,f}}\end{array} \right.$==>${S_r} + i{S_i} = ({M_r} + i{M_i})*({Y_r} + i{Y_i}) = ({M_r}{Y_r} - {M_i}{Y_i}) + i({M_r}{Y_i} + {M_i}{Y_r})$，

那麼：$\left\{ \begin{array}{l}{S_r} = {M_r}{Y_r} - {M_i}{Y_i}\\{S_i} = {M_r}{Y_i} + {M_i}{Y_r}\end{array} \right.$　　解方程得：$\left\{ \begin{array}{l}{M_r} = \frac{{{Y_r}{S_r} + {Y_i}{S_i}}}{{Y_r^2 + Y_i^2}}\\{M_i} = \frac{{{Y_r}{S_i} - {Y_i}{S_r}}}{{Y_r^2 + Y_i^2}}\end{array} \right.$

最終：$M_{cIRM} = {M_r} + i{M_i} = {\frac{{{Y_r}{S_r} + {Y_i}{S_i}}}{{Y_r^2 + Y_i^2}}} + i\frac{{{Y_r}{S_i} - {Y_i}{S_r}}}{{Y_r^2 + Y_i^2}}$

式中，$Y$是帶噪語音，$S$是純淨語音。

優點：cIRM能夠同時增強嘈雜語音的幅度和相位響應，cIRM是加性噪聲假設下的最優掩蔽，可以從帶噪語音中完美重構純淨語音訊號

def cIRM(clean_speech, noise_speech):
    """使用復理想比率掩碼將語音從源訊號的短時傅立葉變換和混合源訊號的短時傅立葉變換中分離出來
    :param clean_speech:純淨語音
    :param noise_speech:帶噪語音
    :return:
    """
    cIRM_r = (np.real(noise_speech) * np.real(clean_speech) + np.imag(noise_speech) * np.imag(clean_speech)) / \
             (np.real(noise_speech) ** 2 + np.imag(noise_speech) ** 2)

    cIRM_i = (np.real(noise_speech) * np.imag(clean_speech) - np.imag(noise_speech) * np.real(clean_speech)) / \
             (np.real(noise_speech) ** 2 + np.imag(noise_speech) ** 2)

    mask = cIRM_r + cIRM_i * 1j
    return mask

總結

　　語音增強中的大部分掩蔽類方法，都可以看成在特定的假設條件下cIRM 的近似。如果將 cIRM 在直角座標系下分解，cIRM 在實數軸上的投影就是 PSM。如果再將 cIRM在極座標系下分解，cIRM 的模值就是 IAM。而 IRM 又是 IAM 在噪聲語音不相關假設下的簡化形式，IBM 則可以認為是 IRM 的二值版本。

各種理想掩蔽的效能比較

度量	IBM	IRM	IAM	PSM	cIRM
PESQ	2.47	3.33	3.45	3.71	4.49
STOI	0.91	0.94	0.97	0.97	1

　　從上表中我們可以看到 cIRM 可以實現對純淨語音幾乎無損地重構，其他掩蔽由於進行了某些特定的假設，所以都會在一定程度上造成效能損失。雖然 cIRM 是最優掩蔽，但是使用其他簡化的掩蔽方法可以降低預測的難度。這也是早期的語音增強研究選擇使用 IBM 或者是 IRM 等簡單掩蔽目標的原因。在模型容量有限的情況下，cIRM 經常並不是最好的選擇，選擇和模型建模能力匹配的目標才能獲得最優的增強效能。

題外話

　　但是，這裡存在一個問題，我們無法從語譜圖中還原語音訊號。為了解決這一問題，我們首先還原所有的頻率分量，即對二值掩碼做個映象後拼接。假設我們計算語譜圖時使用的是512點SFTF，我們一般去前257點進行分析和處理，在這裡我們將前257點的後255做映象，然後拼接在一起得到512點頻率分量，如下圖所示。

然後根據這個還原語音訊號。這裡指的一提的是，在進行STFT後的相位資訊要儲存，用於還原語音訊號。

基於掩蔽的語音增強和基於對映的語音增強模型訓練和增強過程類似，這裡只提幾個重要的地方，其餘地方參考上面內容。

Label：資料的label為根據訊雜比計算的IBM或者IRM，這裡只需要一幀就夠了
損失函式：IBM的損失函式可以用交叉熵，IRM的損失函式還是用均方差
最後一層的啟用函式：IBM只有0和1兩個值，IRM範圍為[0,1]，因此採用sigmoid啟用函式就可以了
重構波形：首先用噪聲幅度譜與計算的Mask值對應位置相乘，程式碼如下，然後根據相位資訊重構語音波形。

enhance_magnitude = np.multiply(magnitude, mask)

Demo效果以及程式碼

首先看下實驗效果，首先是基於對映語音增強的結果：

基於IBM語音增強的結果：

基於IRM語音增強的結果：

訓練程式碼：

"""
@FileName: IBM.py
@Description: Implement IBM
@Author: Ryuk
@CreateDate: 2020/05/08
@LastEditTime: 2020/05/08
@LastEditors: Please set LastEditors
@Version: v0.1
"""

import numpy as np
import librosa
from sklearn.preprocessing import StandardScaler
from keras.layers import *
from keras.models import Sequential


def generateDataset():
    mix, sr = librosa.load("./mix.wav", sr=8000)
    clean,sr = librosa.load("./clean.wav",  sr=8000)

    win_length = 256
    hop_length = 128
    nfft = 512

    mix_spectrum = librosa.stft(mix, win_length=win_length, hop_length=hop_length, n_fft=nfft)
    clean_spectrum = librosa.stft(clean, win_length=win_length, hop_length=hop_length, n_fft=nfft)

    mix_mag = np.abs(mix_spectrum).T
    clean_mag = np.abs(clean_spectrum).T


    frame_num = mix_mag.shape[0] - 4
    feature = np.zeros([frame_num, 257*5])
    k = 0
    for i in range(frame_num - 4):
        frame = mix_mag[k:k+5]
        feature[i] = np.reshape(frame, 257*5)
        k += 1

    snr = np.divide(clean_mag, mix_mag)
    mask = np.around(snr, 0)
    mask[np.isnan(mask)] = 1
    mask[mask > 1] = 1

    label = mask[2:-2]

    ss = StandardScaler()
    feature = ss.fit_transform(feature)
    return feature, label


def getModel():
    model = Sequential()
    model.add(Dense(2048, input_dim=1285))
    model.add(BatchNormalization())

    model.add(LeakyReLU(alpha=0.1))
    model.add(Dropout(0.1))

    model.add(Dense(2048))
    model.add(BatchNormalization())
    model.add(LeakyReLU(alpha=0.1))
    model.add(Dropout(0.1))

    model.add(Dense(2048))
    model.add(BatchNormalization())
    model.add(LeakyReLU(alpha=0.1))
    model.add(Dropout(0.1))

    model.add(Dense(257))
    model.add(BatchNormalization())
    model.add(Activation('sigmoid'))
    return model

def train(feature, label, model):
    model.compile(optimizer='adam',
                  loss='mse',
                  metrics=['mse'])
    model.fit(feature, label, batch_size=128, epochs=20, validation_split=0.1)
    model.save("./model.h5")

def main():
    feature, label = generateDataset()
    model = getModel()
    train(feature, label, model)


if __name__ == "__main__":
    main()

增強程式碼：

"""
@FileName: Inference.py
@Description: Implement Inference
@Author: Ryuk
@CreateDate: 2020/05/08
@LastEditTime: 2020/05/08
@LastEditors: Please set LastEditors
@Version: v0.1
"""

import librosa
import numpy as np
from basic_functions import *
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from keras.models import load_model

def show(data, s):
    plt.figure(1)
    ax1 = plt.subplot(2, 1, 1)
    ax2 = plt.subplot(2, 1, 2)
    plt.sca(ax1)
    plt.plot(data)
    plt.sca(ax2)
    plt.plot(s)
    plt.show()


model = load_model("./model.h5")
data, fs = librosa.load("./test.wav", sr=8000)

win_length = 256
hop_length = 128
nfft = 512

spectrum = librosa.stft(data, win_length=win_length, hop_length=hop_length, n_fft=nfft)
magnitude = np.abs(spectrum).T
phase = np.angle(spectrum).T

frame_num = magnitude.shape[0] - 4
feature = np.zeros([frame_num, 257 * 5])
k = 0
for i in range(frame_num - 4):
    frame = magnitude[k:k + 5]
    feature[i] = np.reshape(frame, 257 * 5)
    k += 1

ss = StandardScaler()
feature = ss.fit_transform(feature)
mask = model.predict(feature)
mask[mask > 0.5] = 1
mask[mask <= 0.5] = 0

fig = plt.figure()
plt.imshow(mask, cmap='Greys', interpolation='none')
plt.show()
plt.close(fig)

magnitude = magnitude[2:-2]
en_magnitude = np.multiply(magnitude, mask)
phase = phase[2:-2]

en_spectrum = en_magnitude.T * np.exp(1.0j * phase.T)
frame = librosa.istft(en_spectrum, win_length=win_length, hop_length=hop_length)

show(data, frame)
librosa.output.write_wav("./output.wav",frame, sr=8000)