位相鋭敏マスク（PSM）で音声強調 with Python

所望信号が既知の場合に観測信号から所望信号を取り出すフィルタ（マスク）を作る位相鋭敏マスク（PSM）をPythonで実装してみる。音声強調の超基本のSTFT領域の音声マスクの設計というやつです。

音源強調の定式化

$x = s + n$

$x$ 観測信号
$s$ 所望信号
$n$ その他信号

なんとかして$s$が欲しい。

音声強調の基本的な考え方としては、
$\hat{S} = G \odot X$

$X$ : $x$をSTFTしたもの
$S$ : $s$をSTFTしたもの
$\hat{S}$ : 推測された所望信号$s$

$G$ : 時間周波数マスク(行列の要素は $0\leq x \leq 1$の実数値)

観測信号$X$に何らかの手段で作った時間周波数マスク$G$を行列の要素ごとに掛け算（アダマール積）することで所望信号$\hat{S}$を推定する。時間周波数マスク$G$の作り方には、所望信号に含まれる周波数を増幅し、その他の信号を減衰させることで所望信号を取り出すウィナーフィルタなど様々に存在する。

音源強調は$x$から$s$を推定する問題となる。

位相鋭敏マスク(PSM)

$G = T[(|S| \oslash |X|) \odot cos(\angle S - \angle X)] ^1 _0$
$T[x]^1_0=min(max(x,0),1)$

コード

PSM.py

# 計算
import numpy as np
# 音声処理
import librosa
import librosa.display
# データ表示
from IPython.display import display, Audio


# Original Data
# オリジナルの音声ファイル
voice, rate = librosa.load("voice.wav", sr=rate)
# オリジナルの音声ファイルにノイズを付加したもの
noise, _ = librosa.load("noise.wav", sr=rate)


# STFT

# STFTパラメーター
n_fft = 1024*4     #フレーム長 (サンプル数)
win_length = n_fft     #窓関数長 (サンプル数)
hop_length = win_length // 4     #ホップ長 (サンプル数)
sr = rate     #サンプリングレート [Hz]

f_voice = librosa.stft(voice, n_fft=n_fft, win_length=win_length, hop_length=hop_length)
f_noise = librosa.stft(noise, n_fft=n_fft, win_length=win_length, hop_length=hop_length)

#PSM
X, S = f_noise, f_voice #観測信号、所望信号

A = (np.abs(S) / np.abs(X)) * np.cos((np.angle(S)-np.angle(X)))
B = np.where(A < 0, 0, A)
G = np.where(B > 1, 1, B)

# スペクトログラムの表示
# librosa.display.specshow(G, y_axis='log', x_axis='time', sr=sr, hop_length=hop_length)

# マスクの適応
_data = f_cafe * G # アダマール積
data = librosa.istft(_test,win_length=win_length, hop_length=hop_length)

# 完成音声の表示
# display(Audio(test, rate=rate))

まとめ

Qiitaには音声ファイルを掲載できないので各自でやってみてください。
そのうちコラボラトリーかなんかでまとめて公開します。

おまけ

ウィナーフィルタを用いたノイズリダクション

$\hat{X} = E[X|Y] = \frac{\sigma_x ^2}{\sigma_x ^2 + \sigma_d ^2} Y$

ここで$\sigma_x ^2$や$\sigma_d ^2$は所望信号と雑音信号のそれぞれのスペクトルの分散を表す。(多分周波数領域ごとのスペクトル分散)

つまり、誰も会話してない部分をノイズとして選択して、Dとして使ってやれば良い。この分散の分数部分がウィナーフィルタとなる。

win.py

s=voice #所望信号
d=noise #観測信号
axis = 1

#観測信号
_, _, X = stft(d,fs = rate,nperseg=nperseg,noverlap = noverlap)

# 所望信号の分散
_, _, S = stft(s,fs = rate,nperseg=nperseg,noverlap = noverlap)
sigma_s=np.square(np.var(S,axis = axis))

#noiseの分散
noise_start = time2tap(8,rate)
noise_end = time2tap(9,rate)

_noise = d[noise_start:noise_end]
_, _, D = stft(_noise,fs = rate,nperseg=nperseg,noverlap = noverlap)
sigma_d=np.square(np.var(D,axis = axis))

#フィルタの生成
W = sigma_s / (sigma_s + sigma_d) #これだけあれば、リアルタイムな処理もできます
G = np.tile(W,(63,1)).T

#フィルタの適応
ret = G * X
ret = ret / np.amax(ret)
_,test = istft(ret,fs=rate)

display(Audio(test,rate=rate))
display(Audio(d, rate=rate3))

参考文献

深層学習に基づく音源強調と位相制御