参考文献
「プログラム101付き 音声信号処理」2021/01/01
(著)川村 新
文献元
準備
console
pip install numpy
ソースコード
sample.py
# -*- coding: utf-8 -*-
import numpy as np
def standardization_(x):
return (x-x.mean())/x.std()
def arc_standardization(x,original):
return (x*original.std())+original.mean()
#///////////////////////////////////////////////
#// //
#// 初期設定 //
#// //
#///////////////////////////////////////////////
def init(St,omega,br,Wr,Wi,FFT_SIZE):
St=int(np.log(FFT_SIZE)/np.log(2.0)+0.5) #// FFTアルゴリズムの段数の算出 +0.5はintへの対応
omega=2.0*np.pi/FFT_SIZE
br,St=bitr(br,St,FFT_SIZE)
Wr,Wi=Wnk(Wr,Wi,omega,FFT_SIZE)
return St,omega,br,Wr,Wi
def bitr(br,St,FFT_SIZE): #// ビット反転
br[:,0]=0 #// 0番目は0
br[:,1]=int(FFT_SIZE/2) #// 1番目はN/2
loop=1
for j in range(St-1):
br[:,loop*2]=int(br[:,loop]/2)
loop=loop*2
for i in range(1,loop,1):
br[:,loop+i]=br[:,loop]+br[:,i]
return br,St
def Wnk(Wr,Wi,omega,FFT_SIZE): #// 重みの計算
i=np.arange(0,int(FFT_SIZE/2))
Wr[:,i]=np.cos(omega*i) #// 重み実部
Wi[:,i]=np.sin(omega*i) #// 重み虚部
return Wr,Wi
#///////////////////////////////////////////////
#// //
#// FFT //
#// //
#///////////////////////////////////////////////
def fft(Fr,Fi,br,xin,St,Wr,Wi,Xr,Xi,Xamp,FFT_SIZE):
_2_s=_2_s_1=roop=l=m=p=0
s=j=k=0
Wxmr=Wxmi=0.
j=np.arange(FFT_SIZE) #// FFT入力の設定
Fr[:,br[:,j]]=xin[:,j] #// 入力
Fi[:,br[:,j]]=0.
_2_s=1
for s in range(1,St+1,1): #// ステージ回繰り返し
_2_s_1=_2_s
_2_s=_2_s*2
roop=int(FFT_SIZE/_2_s)
for j in range(roop): #// DFT繰り返し
for k in range(_2_s_1): #// BF演算繰り返し
l=int(_2_s*j+k ) #// BFの上入力番号
m=int(_2_s_1*(2*j+1)+k) #// BFの下入力番号
p=int(roop*k ) #// 下入力への重み番号
Wxmr=Fr[:,m]*Wr[:,p]+Fi[:,m]*Wi[:,p] #// 重み×下入力の実部
Wxmi=Fi[:,m]*Wr[:,p]-Fr[:,m]*Wi[:,p] #// 重み×下入力の虚部
Xr[:,m]=Fr[:,m]=Fr[:,l]-Wxmr #// BFの下出力の実部
Xi[:,m]=Fi[:,m]=Fi[:,l]-Wxmi #// BFの下出力の虚部
Xr[:,l]=Fr[:,l]=Fr[:,l]+Wxmr #// BFの上出力の実部
Xi[:,l]=Fi[:,l]=Fi[:,l]+Wxmi #// BFの上出力の虚部
j=np.arange(FFT_SIZE)
Xamp[:,j]=np.sqrt(Fr[:,j]*Fr[:,j]+Fi[:,j]*Fi[:,j]) #// 振幅スペクトルの算出
return Fr,Fi,Xr,Xi,Xamp
#///////////////////////////////////////////////
#// //
#// IFFT //
#// //
#///////////////////////////////////////////////
def ifft(Fr,Fi,br,z,St,Wr,Wi,Xr,Xi,FFT_SIZE): #// 逆フーリエ変換
_2_s=_2_s_1=roop=l=m=p=0
s=j=k=0
Wxmr=Wxmi=0.
j=np.arange(FFT_SIZE) #// 逆FFT入力の設定
Fr[:,br[:,j]]=Xr[:,j]
Fi[:,br[:,j]]=Xi[:,j]
_2_s=1
for s in range(1,St+1,1): #// ステージ回繰り返し
_2_s_1=_2_s
_2_s=_2_s*2
roop=int(FFT_SIZE/_2_s)
for j in range(roop): #// FFT繰り返し
for k in range(_2_s_1): #// BF演算繰り返し
l=int(_2_s*j+k ) #// BFの上入力番号
m=int(_2_s_1*(2*j+1)+k) #// BFの下入力番号
p=int(roop*k ) #// 下入力への重み番号
Wxmr=Fr[:,m]*Wr[:,p]-Fi[:,m]*Wi[:,p] #// 重み×下入力の実部
Wxmi=Fi[:,m]*Wr[:,p]+Fr[:,m]*Wi[:,p] #// 重み×下入力の虚部
z[:,m]=Fr[:,m]=Fr[:,l]-Wxmr #// BFの下出力の実部
Fi[:,m]=Fi[:,l]-Wxmi #// BFの下出力の虚部
z[:,l]=Fr[:,l]=Fr[:,l]+Wxmr #// BFの上出力の実部
Fi[:,l]=Fi[:,l]+Wxmi #// BFの上出力の虚部
return Fr,Fi,z
#///////////////////////////////////////////////
#// //
#// メイン関数 //
#// //
#///////////////////////////////////////////////
def MAP_(r):#// MAP推定（レイリー分布）
a=standardization_(r)
b=a/np.abs(a).max()
s=b
MEM_SIZE=b.shape[1] #// 音声メモリのサイズ
FFT_SIZE=np.int64(64) #// FFT点数
#///////////////////////////////////////////////
#// //
#// //
#// FFT関数の構築 //
#// //
#// //
#///////////////////////////////////////////////
Wr =np.full((1,FFT_SIZE+1),0.)
Wi =np.full((1,FFT_SIZE+1),0.) #// FFT 重み
Fr =np.full((1,FFT_SIZE+1),0.)
Fi =np.full((1,FFT_SIZE+1),0.)
Xr =np.full((1,FFT_SIZE+1),0.)
Xi =np.full((1,FFT_SIZE+1),0.) #//実部と虚部
xin =np.full((1,FFT_SIZE+1),0.)
z =np.full((1,FFT_SIZE+1),0.) #// FFT入力と出力
Xamp=np.full((1,FFT_SIZE+1),0.) #// 振幅スペクトル
omega=0. #// FFT角周波数
St=0
br=np.full((1,FFT_SIZE+1),0) #// FFTステージ番号，ビット反転
#//************************************************************************//
#///////////////////////////////
#// //
#// 信号処理用変数 //
#// //
#///////////////////////////////
t = 0 #// 時刻の変数
t_out = 0 #// 終了時刻計測用の変数
#s=np.full((1,FFT_SIZE+1),0.) #// 入力データ格納用変数
y=np.full((r.shape[0],r.shape[1]),0.) #// 出力データ格納用変数
l=0
SHIFT = int(FFT_SIZE/2); #// FFTのシフト量
OV = 2.0*SHIFT/FFT_SIZE #// オーバーラップ加算の係数
x=np.full((1,FFT_SIZE+1),0.) #// FFTの入力
yf=np.full((1,FFT_SIZE+1),0.) #// IFFT信号格納用
w=np.full((1,FFT_SIZE+1),0.) #// 窓関数
_lambda=np.full((1,FFT_SIZE+1),0.) #// ノイズ推定値
G=np.full((1,FFT_SIZE+1),0.) #// スペクトルゲイン
NM=0 #// ノイズを推定するフレームの数
cnt = 0 #// フレーム番号のカウンタ
gamma=np.full((1,FFT_SIZE+1),0.)
gamma1=np.full((1,FFT_SIZE+1),0.) #// 事後SNR
xi=np.full((1,FFT_SIZE+1),0.) #// 事前SNR
tmp=0. #// 途中計算で利用
#///////////////////////////////
#// //
#// 変数の初期設定 //
#// //
#///////////////////////////////
St,omega,br,Wr,Wi=init(St,omega,br,Wr,Wi,FFT_SIZE) #// ビット反転，重み係数の計算
l = 0 #// FFT開始時刻管理
NM= 1024*16 #// 初期ノイズ推定のフレーム数
i=np.arange(FFT_SIZE) #// 窓関数の設定
w[:,i]=0.5*(1.0-np.cos(2.0*np.pi*i/FFT_SIZE))
add_len=b.shape[1]-1
#//************************************************************************//
#///////////////////////////////////
#// //
#// メインループ //
#// //
#///////////////////////////////////
while True: #// メインループ
if t_out > add_len:
break # ループ終了判定
#//************************************************************************//
#////////////////////////////////////////////////////
#// //
#// Signal Processing //
#// //
#// 現在時刻tの入力 s[t] から出力 y[t] をつくる //
#// //
#// ※ tは0からMEM_SIZE-1までをループ //
#// //
#////////////////////////////////////////////////////
x[:,l] = s[:,t] #// 入力をx[l]に格納
l=(l+1)%FFT_SIZE #// FFT用の時刻管理
if l%SHIFT==0: #// シフトごとにFFTを実行
for i in range(FFT_SIZE):
xin[:,i] = x[:,(l+i)%FFT_SIZE]*w[:,i] #// 窓関数を掛ける
Fr,Fi,Xr,Xi,Xamp=fft(Fr,Fi,br,xin,St,Wr,Wi,Xr,Xi,Xamp,FFT_SIZE) #// FFT
#print(Fr,Fi,br,xin,St,Wr,Wi,Xr,Xi,Xamp)
#///////////////////////////
#// //
#// 雑音推定 //
#// //
#///////////////////////////
if cnt<NM: #// 初期NMフレームの平均がノイズ推定値
for i in range(FFT_SIZE):
_lambda[:,i]+=Xamp[:,i]*Xamp[:,i]/NM
cnt+=1 #// フレーム番号更新
for i in range(FFT_SIZE):
#///////////////////////////
#// //
#// 事前＆事後SNR推定 //
#// //
#///////////////////////////
gamma1[:,i] = gamma[:,i] #// 過去のγを記録
if _lambda[:,i]!=0:
gamma[:,i]=Xamp[:,i]*Xamp[:,i]/_lambda[:,i] #// 事後SNR γの更新
#// 事前SNR ξの更新．Decision Directed method
xi[:,i] = 0.98*gamma1[:,i]*G[:,i]*G[:,i] + 0.02*( np.fabs(gamma[:,i]-1.0) + (gamma[:,i]-1.0) )*0.5
#////////////////////////////////////////////
#// //
#// MAP推定 //
#// 音声振幅：レイリー分布 //
#// 音声位相：一様分布 //
#// ノイズ ：ガウス分布 //
#////////////////////////////////////////////
if gamma[:,i]!=0: #// ゲインの根号内を計算
tmp=xi[:,i]*xi[:,i]+2.0*(1.0+xi[:,i])*(xi[:,i]/gamma[:,i])
if tmp>0: #// 根号内が正ならMAPゲインを計算
G[:,i]=0.5*(xi[:,i]+np.sqrt(tmp))/(1.0+xi[:,i])
G[:,i]=( np.fabs(G[:,i])+G[:,i] )/2.0 #// 負のゲインは0にする
if G[:,i]>1.0:
G[:,i]=1.0 #// ゲインの最大値は1
Xr[:,i]= G[:,i]*Xr[:,i] #// 実部にゲインを乗じる
Xi[:,i]= G[:,i]*Xi[:,i] #// 虚部にゲインを乗じる
Fr,Fi,z=ifft(Fr,Fi,br,z,St,Wr,Wi,Xr,Xi,FFT_SIZE) #// IFFT
#print(Fr,Fi,br,z,St,Wr,Wi,Xr,Xi)
for i in range(FFT_SIZE): #// 出力信号作成
if i>=FFT_SIZE-SHIFT:
yf[:,(l+i)%FFT_SIZE] = z[:,i]/FFT_SIZE*OV
else:
yf[:,(l+i)%FFT_SIZE] = yf[:,(l+i)%FFT_SIZE]+z[:,i]/FFT_SIZE*OV
y[:,t]=yf[:,l] #// 現在の出力
#//************************************************************************//
y[:,t] = np.arctan(y[:,t])/(np.pi/2.0) #// クリップ防止
y[:,t] = y[:,t]*np.abs(a).max() #// 出力を整数化
t=(t+1)%MEM_SIZE #// 時刻 t の更新
t_out+=1 #// ループ終了時刻の計測
return arc_standardization(y,r) #// メイン関数の終了