使用譜減法去除音頻底噪

上一篇文章我主要分享了短時傅立葉變換及其逆變換在python中的實現，有興趣的可以閱讀一下該篇文章，地址如下：
Python音頻信號處理 1.短時傅里葉變換及其逆變換

那么在本篇文章中，我們將利用短時傅立葉變換及其逆變換來實現譜減法。

Part 1 ：譜減法

譜減法的核心思路非常簡單，顧名思義，譜減法是一種頻域上的信號處理方法，其基本思路就是提取出信號本身的頻譜以及噪音的頻譜，通過兩者之差獲取降噪后信號的頻譜，最后利用傅立葉變換逆變換重構初始信號。

Part 2 ：使用譜減法降低或消除信號的底噪

要使用譜減法來進行信號處理，顯然我們首先需要計算出信號的頻譜以及噪音的頻譜。這里我們繼續使用上篇文章所介紹的短時傅立葉變換及其逆變換的實現，Python代碼如下：

# 短時傅立葉變換
def TFCT(trame, Fe, Nfft,fenetre,Nwin,Nhop):L = round((len(trame) - len(fenetre))/Nhop)+1M = Nfftxmat = np.zeros((M,L))print('xmat',xmat.shape)print(Nwin+Nhop)for j in range(L):xmat[:,j] = np.fft.fft(trame[j*Nhop:Nwin+Nhop*j]*window,Nfft)  x_temporel = np.linspace(0,(1/Fe)*len(trame),len(trame))x_frequentiel = np.linspace(0, Fe,Nfft)return xmat,x_temporel,x_frequentiel

# 短時傅立葉變換逆變換
def ITFD(xmat,Fe,Nfft,Nwin,Nhop):window = np.hamming(Nwin)Te = 1/Feyvect = np.zeros(Nfft + (xmat.shape[1]-1)*Nhop)t_vecteur = np.arange(0,Te*len(yvect),Te)index = 0K = 0L = xmat.shape[1]yl = np.zeros((Nfft,L))for j in range(L):yl[:,j] = np.fft.ifft(xmat[:,j])for k in range(L):yvect[Nhop*k:Nfft+Nhop*k] += yl[:,k]for n in range(Nwin-1):K +=  window[n]K /= Nhopyvect /=Kprint(yvect.shape)return t_vecteur, yvect

有關于短時傅立葉變換及其逆變換的實現這里就不過多贅述，有興趣的可以參考上一篇文章，其中詳細解釋了代碼的構成。

首先我們先讀取信號，并觀察其在時域中的圖像：

# 讀取初始信號 mix.wav
Fe, mix = wavfile.read('mix.wav')
mix = mix/2**15
Te = 1 / Fe
xtemp = np.arange(0,Te*len(mix),Te)
plt.figure()
plt.plot(xtemp,mix)
plt.axis([0,1,-0.6,0.6])
plt.xlabel('t(s)')
plt.ylabel('amplitude(V)')
plt.title('representation of the original signal in time domain')
display(Audio(mix,rate=Fe))

Fe為原始信號的采樣頻率，mix為復信號值。
Te = 1/Fe = 6.25e-5

我們不難發現，0-0.4s，時域上只有底噪，因此我們可以利用這段時間來計算噪音的頻譜。
這個例子中，我們選擇頻域的采樣點個數Nfft=1024，窗函數的長度Nwin=1024，每次窗函數滑動的長度Nhop=512，計算該信號的短時傅立葉變換。

Nfft = 1024
Nwin = 1024
Nhop = 512
window = np.hamming(Nwin)
xmat_sound,tvect,fvect = TFCT(mix,Fe,Nfft,window,Nwin,Nhop)

xmat_sound就是我們保存短時傅立葉變換值的矩陣，該矩陣的每一行代表一個在0至采樣頻率范圍內的頻率，單位為Hz，每一列對應該段被窗函數截取的信號的FFT快速傅里葉變換值。
不難算出第一列代表的時域范圍為0-Nwin*（1/Fe）= 0.064s
在該例中，我們就使用第一列的頻譜，近似認為是底噪的頻譜。（顯然我們這里應該取盡可能多的純噪音頻譜，并使用他們幅值的平均值，這里為了簡化僅使用第一個窗函數截取的信號部分作為噪聲頻譜。）

我們可以先使用imshow畫出原始信號的光譜圖，橫軸為時間，縱軸為頻率。

module_tf_xmat = abs(xmat_sound)
plt.figure()
xlim = int(module_tf_xmat.shape[0]/2)
ylim = int(module_tf_xmat.shape[1]/2)
plt.imshow(20*np.log10(module_tf_xmat[0:xlim,:]),extent=[0,Te*len(mix),0,Fe/2],aspect='auto')
plt.colorbar()
plt.xlabel('time(s)')
plt.ylabel('frequence(Hz)')
plt.title('spectrogram of the originql signal')

接下來我們使用譜減法，對信號頻譜所有列減去第一列對應的噪音頻譜，注意這里的全部減法都是針對幅值。

module_tf_xmat = abs(xmat_sound)
angle_tf_xmat = np.angle(xmat_sound)
module_tf_bruit = module_tf_xmat[:,0]module_reconstruit = np.zeros(module_tf_xmat.shape)
for n in range(module_tf_xmat.shape[1]):module_reconstruit[:,n] = module_tf_xmat[:,n] - module_tf_bruit
module_reconstruit[module_reconstruit<0] = 0

譜減法之后的光譜圖如下所示：

最后我們使用短時傅立葉變換逆變換將獲取的降噪后的幅值矩陣，使用原始的信號相位，重構降噪后的復信號。

# 將相位和降噪后的幅值重構復信號的頻域分布
tf_reconstruit = np.zeros(module_tf_xmat.shape,dtype=complex)
for i in range(module_tf_xmat.shape[0]):for j in range(module_tf_xmat.shape[1]):tf_reconstruit[i,j] = module_reconstruit[i,j] * np.exp(angle_tf_xmat[i,j]*1j)
# 使用短時傅立葉變換逆變換重構時域內的信號
t,yvect = ITFD(tf_reconstruit,Fe,Nfft,Nwin,Nhop)

總結

在本例中，我們僅僅使用了短時傅立葉變換矩陣的第一列來計算噪聲的頻譜，實際是非常不準確的，因為此時所取的噪聲長度僅為0.064s。即使是對于平穩的噪音，這段噪聲的頻譜仍然不能很好地代表噪音的真實頻譜。在實際的信號處理過程中，較好的方式是取足夠長度的噪音并計算其幅值的平均值，再使用這個噪音在各個頻率幅值的平均值進行譜減法，可以更加有效地降低底噪。