直交射影とウィナーフィルタ

19580 ワード
音声信号処理 音声信号処理 テキストリンク

音源による直交射影

下記によると音源による直交射影を求める方法が書いてある。

Performance measurement in blind audio source separation

\begin{align}
S_{target} &:= \boldsymbol P_{s_j} \hat s_j \\
\boldsymbol P_{s_j} &:= \Pi \{(s^\tau_j)_{0\leq\tau\leq\mathit L - 1}\} \\
s^\tau_j &= s_j(t-\tau)\\
\\
\\
\text{特に$L$=1の場合}\\

S_{target} &= \langle \hat s_j, s_j \rangle s_j / \|s_j\|^2 \\
\langle a, b \rangle &:= \sum^{T-1}_{t=0} a(t) \overline b(t)
\end{align}
\\

これにより，音源分離時の元の音源の成分を抽出することで，どの程度分離がうまくいっているか評価しているのだけれども，
直交射影でそんなにうまくいくのかなと疑問に思っていました。

この間，ウィナーフィルタを勉強したのでこっちの方がいいのではと，導出方法を見直していたら，

J = \sum_{k=1}^{L_s}|\varepsilon_k|^2
  = \sum_{k=1}^{L_s}|d_k-\boldsymbol w^H \boldsymbol u_k|^2 \\

簡単にするため$\boldsymbol w$の係数を１つにして，

J = \sum_{k=1}^{L_s}|d_k- w u_k|^2 \\

Jを最小化する$w$は

\begin{align}
\frac{\partial J}{\partial w} &= -\sum_{k=1}^{L_s} u_k(d_k- w u_k) \\
                              &= -\sum_{k=1}^{L_s} (u_k d_k - w u_k^2) \\
                              &= w \sum_{k=1}^{L_s} u_k^2 -\sum_{k=1}^{L_s} u_k d_k = 0 \\
\end{align}

整理して

\begin{align}
  w &= \frac{\sum_{k=1}^{L_s} u_k d_k}{\sum_{k=1}^{L_s} u_k^2} \\
    &= \frac{\langle d, u \rangle}{\| u_k \|^2} \\
  y &= w u = \langle d, u \rangle u / \| u_k \|^2 
\end{align}

で，実は同じことを示していた。。。

複数音源のウィナーフィルタ

これだけだと，悔しいのでもう少し考察。

mir_eval/separation.py
この中の_project が直交射影を求める計算をしているが，自分には何をしているかわからなかった。

なので，$n$個の音源の場合のウィナーフィルタで考えてみた。

\begin{align}
  \boldsymbol J & = 
      (\boldsymbol d - \sum_{j=1}^{n} \boldsymbol w_j^H \boldsymbol U_j)
      (\boldsymbol d - \sum_{j=1}^{n} \boldsymbol w_j^H \boldsymbol U_j)^H \\
                & = 
      (\boldsymbol d   - \sum_{j=1}^{n} \boldsymbol w_j^H \boldsymbol U_j)
      (\boldsymbol d^H - \sum_{j=1}^{n} \boldsymbol U_j^H \boldsymbol w_j) \\
\\
 \frac{\partial \boldsymbol J}{\partial \boldsymbol w_i^*}
    &= -\boldsymbol U_i(\boldsymbol d^H - \sum_{j=1}^{n} \boldsymbol U_j^H \boldsymbol w_j) \\
    &= -\boldsymbol U_i \boldsymbol d^H + \sum_{j=1}^{n} \boldsymbol U_i \boldsymbol U_j^H \boldsymbol w_j \\
    &= \boldsymbol 0\\
\\
\\
\text{ただし}\\
  {\boldsymbol U}_i &=  


\begin{bmatrix}
   u_{i,1} & u_{i,2} & \cdots & u_{i,N}   & \cdots & u_{i,L-1} & u_{i,L}     & 0           & \cdots & 0       \\
   0       & u_{i,1} & \cdots & u_{i,N-1} & \cdots & u_{i,L-2} & u_{i,L-1}   & u_{i,L}     & \cdots & 0       \\
   \vdots  &         & \ddots & \vdots    &        & \vdots    & \vdots      &             & \ddots & \vdots  \\
   0       & 0       & \cdots & u_{i,1}   & \cdots & u_{i,L-N} & u_{i,L-N+1} & u_{i,L-N+2} & \cdots & u_{i,L} \\
\end{bmatrix} \\
  \\

  {\boldsymbol d} &= [d_1, \cdots, d_{L}, 0, \cdots , 0] \\

\end{align}

並べなおすと

\begin{array}{c}
  \boldsymbol U_1 \boldsymbol U_1^H \boldsymbol w_1 & + & \boldsymbol U_1 \boldsymbol U_2^H \boldsymbol w_2 & + & \cdots & + & \boldsymbol U_1 \boldsymbol U_n^H \boldsymbol w_n & = & \boldsymbol U_1 \boldsymbol d^H \\
  \boldsymbol U_2 \boldsymbol U_1^H \boldsymbol w_1 & + & \boldsymbol U_2 \boldsymbol U_2^H \boldsymbol w_2 & + & \cdots & + & \boldsymbol U_2 \boldsymbol U_n^H \boldsymbol w_n & = & \boldsymbol U_2 \boldsymbol d^H \\
  \vdots                                            &   & \vdots                                            &   & \ddots &   & \vdots                                            & = & \vdots \\
  \boldsymbol U_n \boldsymbol U_1^H \boldsymbol w_1 & + & \boldsymbol U_n \boldsymbol U_2^H \boldsymbol w_2 & + & \cdots & + & \boldsymbol U_n \boldsymbol U_n^H \boldsymbol w_n & = & \boldsymbol U_n \boldsymbol d^H \\
\end{array}

これを行列で書くと

\left(
\begin{array}{c}
  \boldsymbol U_1 \boldsymbol U_1^H & \boldsymbol U_1 \boldsymbol U_2^H & \cdots & \boldsymbol U_1 \boldsymbol U_n^H \\
  \boldsymbol U_2 \boldsymbol U_1^H & \boldsymbol U_2 \boldsymbol U_2^H & \cdots & \boldsymbol U_2 \boldsymbol U_n^H \\
  \vdots                            & \vdots                            & \ddots & \vdots                            \\
  \boldsymbol U_n \boldsymbol U_1^H & \boldsymbol U_n \boldsymbol U_2^H & \cdots & \boldsymbol U_n \boldsymbol U_n^H \\
\end{array}
\right)
\left(
\begin{array}{c}
  \boldsymbol w_1 \\
  \boldsymbol w_2 \\
  \vdots          \\
  \boldsymbol w_n \\
\end{array}
\right)
=
\left(
\begin{array}{c}
  \boldsymbol U_1 \boldsymbol d^H \\
  \boldsymbol U_2 \boldsymbol d^H \\
  \vdots \\
  \boldsymbol U_n \boldsymbol d^H \\
\end{array}
\right)

なので，_project 関数内の下記$G$

    # inner products between delayed versions of reference_sources
    G = np.zeros((nsrc * flen, nsrc * flen))
    for i in range(nsrc):
        for j in range(nsrc):
            ssf = sf[i] * np.conj(sf[j])
            ssf = np.real(scipy.fftpack.ifft(ssf))
            ss = toeplitz(np.hstack((ssf[0], ssf[-1:-flen:-1])),
                          r=ssf[:flen])
            G[i * flen: (i+1) * flen, j * flen: (j+1) * flen] = ss
            G[j * flen: (j+1) * flen, i * flen: (i+1) * flen] = ss.T

は

\boldsymbol G = 
\left(
\begin{array}{c}
  \boldsymbol U_1 \boldsymbol U_1^H & \boldsymbol U_1 \boldsymbol U_2^H & \cdots & \boldsymbol U_1 \boldsymbol U_n^H \\
  \boldsymbol U_2 \boldsymbol U_1^H & \boldsymbol U_2 \boldsymbol U_2^H & \cdots & \boldsymbol U_2 \boldsymbol U_n^H \\
  \vdots                            & \vdots                            & \ddots & \vdots                            \\
  \boldsymbol U_n \boldsymbol U_1^H & \boldsymbol U_n \boldsymbol U_2^H & \cdots & \boldsymbol U_n \boldsymbol U_n^H \\
\end{array}
\right) \\

を計算していて，下記$D$

   # inner products between estimated_source and delayed versions of
    # reference_sources
    D = np.zeros(nsrc * flen)
    for i in range(nsrc):
        ssef = sf[i] * np.conj(sef)
        ssef = np.real(scipy.fftpack.ifft(ssef))
        D[i * flen: (i+1) * flen] = np.hstack((ssef[0], ssef[-1:-flen:-1]))

は

\boldsymbol D =
\left(
\begin{array}{c}
  \boldsymbol U_1 \boldsymbol d^H \\
  \boldsymbol U_2 \boldsymbol d^H \\
  \vdots \\
  \boldsymbol U_n \boldsymbol d^H \\
\end{array}
\right)

を計算し，

\boldsymbol C = 
\left(
\begin{array}{c}
  \boldsymbol w_1 \\
  \boldsymbol w_2 \\
  \vdots          \\
  \boldsymbol w_n \\
\end{array}
\right)

を，

    # Computing projection
    # Distortion filters
    try:
        C = np.linalg.solve(G, D).reshape(flen, nsrc, order='F')
    except np.linalg.linalg.LinAlgError:
        C = np.linalg.lstsq(G, D)[0].reshape(flen, nsrc, order='F')

で計算して求めているみたい。
ただ，$(flen*nsrc \times 1)$だと使いにくいのでreshapeで$(flen \times nsrc)$に変換しているみたい。

直交射影とフーリエ変換

なんとなく，$G,D$が何を計算しているかわかってきたけど，"ss = toeplitz(...)" や "np.hstack((ssef[0], ssef[-1:-flen:-1]))"が気持ち悪い。
$U_iU_j^H$や$U_id^H$が何をしているかもう少し見てみる



\begin{align}
&\text{$\boldsymbol A$と$\boldsymbol B$を下記のように定義。}\\


\boldsymbol A &= 

\begin{bmatrix}
   a_{1}       & a_{2}       & \cdots & a_{N}       & \cdots & a_{L-1} & a_{L}     & 0         & \cdots  & 0      \\
   0           & a_{1}       & \cdots & a_{N-1}     & \cdots & a_{L-2} & a_{L-1}   & a_{L}     & \cdots  & 0      \\
  \vdots       &             & \ddots &  \vdots     &        & \vdots  & \vdots　  &           & \ddots  & \vdots \\
   0           & 0           & \cdots & a_{1}       & \cdots & a_{L-N} & a_{L-N+1} & a_{L-N+2} & \cdots  & a_{L}  \\
\end{bmatrix} \\
\\
\boldsymbol B &= 

\begin{bmatrix}
   b_{1}       & b_{2}       & \cdots & b_{M}       & \cdots & b_{L-1} & b_{L}     & 0         & \cdots  & 0      \\
   0           & b_{1}       & \cdots & b_{M-1}     & \cdots & b_{L-2} & b_{L-1}   & b_{L}     & \cdots  & 0      \\
  \vdots       &             & \ddots &  \vdots     &        & \vdots  & \vdots    &           & \ddots  & \vdots \\
   0           & 0           & \cdots & b_{1}       & \cdots & b_{L-M} & b_{L-M+1} & b_{L-M+2} & \cdots  & b_{L}  \\
\end{bmatrix} \\
\\
\text{ただし，} & M \geqq N \\
                & a_i = b_i = 0 \  (i \leq 0, \  i>L) \\
                & \hat L = L + M - 1 \\

\text{とする。} & \text{この時，$\boldsymbol A \boldsymbol B^H$ は} \\ 
\\
\boldsymbol A \boldsymbol B^H &= 
\begin{bmatrix}
   \sum^{\hat L}_{i=1} a_{i}b_{i}       & \sum^{\hat L}_{i=1} a_{i}b_{i-1}        & \cdots & \sum^{\hat L}_{i=1} a_{i}b_{i-M+1}     \\
   \sum^{\hat L}_{i=1} a_{i-1}b_{i}     & \sum^{\hat L}_{i=1} a_{i-1}b_{i-1}      & \cdots & \sum^{\hat L}_{i=1} a_{i-1}b_{i-M+1}   \\
   \vdots                               &                                         & \ddots & \vdots                                 \\
   \sum^{\hat L}_{i=1} a_{i-N+1}b_{i}   & \sum^{\hat L}_{i=1} a_{i-N+1}b_{i-1}    & \cdots & \sum^{\hat L}_{i=1} a_{i-N+1}b_{i-M+1} \\
\end{bmatrix} \\
\text{となる。}& \\
\\
\text{特に}&\text{$M=1$の場合}\\
\boldsymbol A \boldsymbol B^H &= 
\begin{bmatrix}
   \sum^{\hat L}_{i=1} a_{i}b_{i}     \\
   \sum^{\hat L}_{i=1} a_{i-1}b_{i}   \\
   \vdots                             \\
   \sum^{\hat L}_{i=1} a_{i-N+1}b_{i} \\
\end{bmatrix} \\
\text{となる。}& \\


\end{align}

ここで，$\hat a_k=FFT(a_n, \hat L),$ $\hat b_k=FFT(b_n, \hat L),$ $c(n)=IFFT(\hat a_k \hat b_k^*, \hat L) $とした場合，
$c_n = \sum_{i=1}^{\hat L} a_{i+n} b_{i} = \sum_{i=1}^{\hat L} a_{i} b_{i-n}$となる。(範囲外は$0$なので)

よって

\begin{align}
\boldsymbol A \boldsymbol B^H &= 
\begin{bmatrix}
   c_0       & c_{1}        & \cdots & c_{N-1} & \cdots & c_{M-1}  \\
   c_{-1}    & c_0          & \cdots & c_{N-2} & \cdots & c_{M-2}  \\
   \vdots    &              & \ddots & \vdots  & \ddots & \vdots   \\
   c_{-N+1}  & c_{-N+2}     & \cdots & c_{0}   & \cdots & c_{M-N}  \\
\end{bmatrix} \\
&= toeplitz([c_0, c_{-1}, \dots, c_{-N+1}],~ r=[c_0, c_1, \dots, c_{M-1}])\\
&= toeplitz(np.hstack((c[0], c[-1:-N:-1])), r=c[:M])\\
\text{となる。}& \text{(M>Nの場合)}\\
\\
\text{特に}&\text{$M=1$の場合}\\
\boldsymbol A \boldsymbol B^H &= 
\begin{bmatrix}
   c_0      \\
   c_{-1}   \\
   \vdots   \\
   c_{-N+1} \\
\end{bmatrix} \\
&=np.hstack((c[0], c[-1:-N:-1])) \\
\text{となる。}& \\

\end{align}

あとは，$A,B$に$U,d$を当てはめればOK。しかし，共役にする方を逆にすれば，多少見やすくなる気がするのだが，どうだろう。。。
それと，toeplitzの仕様も若干気に入らない。。。