Python语音基础操作--5.2谱减法

《语音信号处理试验教程》（梁瑞宇等）的代码主要是Matlab实现的，现在Python比较热门，所以把这个项目大部分内容写成了Python实现，大部分是手动写的。使用CSDN博客查看帮助文件：

Python语音基础操作–2.1语音录制，播放，读取
Python语音基础操作–2.2语音编辑
Python语音基础操作–2.3声强与响度
Python语音基础操作–2.4语音信号生成
Python语音基础操作–3.1语音分帧与加窗
Python语音基础操作–3.2短时时域分析
Python语音基础操作–3.3短时频域分析
Python语音基础操作–3.4倒谱分析与MFCC系数
Python语音基础操作–4.1语音端点检测
Python语音基础操作–4.2基音周期检测
Python语音基础操作–4.3共振峰估计
Python语音基础操作–5.1自适应滤波
Python语音基础操作–5.2谱减法
Python语音基础操作–5.4小波分解
Python语音基础操作–6.1PCM编码
Python语音基础操作–6.2LPC编码
Python语音基础操作–6.3ADPCM编码
Python语音基础操作–7.1帧合并
Python语音基础操作–7.2LPC的语音合成
Python语音基础操作–10.1基于动态时间规整(DTW)的孤立字语音识别试验
Python语音基础操作–10.2隐马尔科夫模型的孤立字识别
Python语音基础操作–11.1矢量量化(VQ)的说话人情感识别
Python语音基础操作–11.2基于GMM的说话人识别模型
Python语音基础操作–12.1基于KNN的情感识别
Python语音基础操作–12.2基于神经网络的情感识别
Python语音基础操作–12.3基于支持向量机SVM的语音情感识别
Python语音基础操作–12.4基于LDA，PCA的语音情感识别

代码可在Github上下载：busyyang/python_sound_open

谱减法

对于任何一帧信号$x_i(m)$做FFT变换后：
$$X_i(k)=\sum _{m=1}^N{x}_i(m)\exp (j\frac{2\pi mk}{N})$$

对于$X_i(k)$的幅值为$|X_i(k)|$,角度为$X_{angle}^i(k)=\arctan [\frac{Im(X_i(k))}{Re(X_i(k))}]$,前导噪声段时长为IS,对应帧数为NIS,可以得到该噪声段的平均能量为：
$$D(k)=\frac{1}{NIS}\sum _{i=1}^{NIS}|X_i(k){|}^2$$

谱减公式为：
$$|{\hat{X}}_i(k){|}^2=\{\begin{array}{ll}|X_i(k){|}^2-a\times D(k)& |X_i(k){|}^2⩾a\times D(k)\\ b\times D(k)& |X_i(k){|}^2<a\times D(k)\end{array}$$

其中，$a,b$是两个常数，$a$为过减因子，$b$为增益补偿因子。

利用谱减后的幅值$|{\hat{X}}_i(k)|$,以及原先的相位角$X_{angle}^i(k)$，可以利用iFFT求出增强后的语音序列${\hat{x}}_i(m)$。

Boll改进谱减法

（一）谱减公式改为：
$$|{\hat{X}}_i(k){|}^{\gamma }=\{\begin{array}{ll}|X_i(k){|}^{\gamma }-a\times D(k)& |X_i(k){|}^{\gamma }⩾a\times D(k)\\ b\times D(k)& |X_i(k){|}^{\gamma }<a\times D(k)\end{array}$$

$$D(k)=\frac{1}{NIS}\sum _{i=1}^{NIS}|X_i(k){|}^{\gamma }$$

当$\gamma =1$，算法相当于用谱幅值做谱减法，当$\gamma =2$，算法相当于用功率谱幅值做谱减法。

（二）计算平均谱值代替
$$Y_i(k)=\frac{1}{2M+1}\sum _{j=-M}^M{X}_{i+j}(k)$$

使用$Y_i(k)$代替$X_i(k)$，可以得到较小的谱估算方差。

（三）减小噪声残留
D i ( k ) = { D i ( k ) D i ( k ) ⩾ max ⁡ ∣ N R ( k ) ∣ min ⁡ { D j ( k ) ∣ j ∈ [ i − 1 , i , i + 1 ] } D i ( k ) < max ⁡ ∣ N R ( k ) ∣ D_i(k)=\left \{\begin{array}{ll} D_i(k)& D_i(k)\geqslant \max|N_R(k)|\\ \min\{D_j(k)|j \in [i-1,i,i+1]\}&D_i(k)< \max|N_R(k)| \end{array} \right. Di​(k)={Di​(k)min{Dj​(k)∣j∈[i−1,i,i+1]}​Di​(k)⩾max∣NR​(k)∣Di​(k)<max∣NR​(k)∣​

其中，$\max |N_R(k)|$为最大残余噪声。

from chapter2_基础.soundBase import *
from chapter5_语音降噪.自适应滤波 import *


def awgn(x, snr):
    snr = 10 ** (snr / 10.0)
    xpower = np.sum(x ** 2) / len(x)
    npower = xpower / snr
    return x + np.random.randn(len(x)) * np.sqrt(npower)


data, fs = soundBase('C5_1_y.wav').audioread()
data -= np.mean(data)
data /= np.max(np.abs(data))
IS = 0.25  # 设置前导无话段长度
wlen = 200  # 设置帧长为25ms
inc = 80  # 设置帧移为10ms
SNR = 5  # 设置信噪比SNR
N = len(data)  # 信号长度
time = [i / fs for i in range(N)]  # 设置时间
r1 = awgn(data, SNR)
NIS = int((IS * fs - wlen) // inc + 1)
# 5.2.1
snr1 = SNR_Calc(r1, r1 - data)
a, b = 4, 0.001
output = SpectralSub(r1, wlen, inc, NIS, a, b)
if len(output) < len(r1):
    filted = np.zeros(len(r1))
    filted[:len(output)] = output
elif len(output) > len(r1):
    filted = output[:len(r1)]
else:
    filted = output

plt.subplot(4, 1, 1)
plt.plot(time, data)
plt.ylabel('原始信号')
plt.subplot(4, 1, 2)
plt.plot(time, r1)
plt.ylabel('加噪声信号')
plt.subplot(4, 1, 3)
plt.ylabel('滤波信号')
plt.plot(time, filted)

plt.show()