Spectrum Subtraction

スペクトルサブトラクションとは

スペクトルサブトラクションとは, その名のとおり, 周波数領域で減算処理をすることで, 背景雑音を除去するアルゴリズムです. 音データを $x(n)$, 背景雑音を ${\omega}(n)$ とすると, 背景雑音が混入した音データ $y(n)$ は以下の式で定義できます.

$ \begin{eqnarray} y(n)=x(n)+{\omega}(n) \end{eqnarray} $

つまり, 背景雑音を除去するには, 以下の式で定義できます.

$ \begin{eqnarray} x(n)=y(n)-{\omega}(n) \end{eqnarray} $

この減算処理を, 周波数領域で実行するのが, スペクトラムサブトラクションです.

$x(n)$, ${\omega}(n)$, $y(n)$ の周波数領域の信号をそれぞれ, $X(k)$, $W(k)$, $Y(k)$ と定義すると,

$ \begin{eqnarray} X(k)=Y(k)-W(k) \end{eqnarray} $

ここで, $Y(k)$, $W(k)$ をそれぞれ極座標形式で定義すると,

$ \begin{eqnarray} Y(k)=A_{Y}(k)\exp(j{\theta}_{Y}(k)) \end{eqnarray} $

$ \begin{cases} \begin{eqnarray} A_{Y}(k) &=&\sqrt{real(Y(k))^{2}+imag(Y(k))^{2}} \\
\theta_{Y}(k)&=&\tan^{-1}\left(\frac{imag(Y(k))}{real(Y(k))}\right) \\
\end{eqnarray} \end{cases} $

$ \begin{eqnarray} W(k)=A_{W}(k)\exp(j{\theta}_{W}(k)) \end{eqnarray} $

$ \begin{cases} \begin{eqnarray} A_{W}(k) &=&\sqrt{real(W(k))^{2}+imag(W(k))^{2}} \\
\theta_{W}(k)&=&\tan^{-1}\left(\frac{imag(W(k))}{real(W(k))}\right) \\
\end{eqnarray} \end{cases} $

したがって, スペクトルサブトラクションの定義式は以下のように定義できます.

$ \begin{eqnarray} X(k)=A_{Y}(k)\exp(j{\theta}_{Y}(k))-A_{W}(k)\exp(j{\theta}_{W}(k)) \end{eqnarray} $

さらに, 人間の聴覚は, 位相スペクトルのちがいに鈍感 という特性を考慮すると, 位相スペクトルを無視することが可能で,

$ \begin{eqnarray} X(k)=(A_{Y}(k)-A_{W}(k))\exp(j{\theta}_{Y}(k)) \end{eqnarray} $

と定義できます.

ところで, 一般的に, 背景雑音の振幅は未知です. そこで, 音データが 0 (無音) の部分から背景雑音の振幅スペクトルを推定した, 一般的な, スペクトラムサブトラクションの定義は以下のようになります.

$ X(k)=\begin{cases} \begin{eqnarray} &(A_{Y}(k)-\overline{A_{W}(k)})\exp(j{\theta}_{Y}(k))& \quad (A_{Y}(K)-\overline{A_{W}(k)} \geqq 0) \\
&0& \quad (A_{Y}(K)-\overline{A_{W}(k)} < 0) \\
\end{eqnarray} \end{cases} $

一般的に, 音データの周波数特性は時間の経過とともに変化するので, スペクトラムサブトラクションではフレーム単位で周波数特性を加工します.

フレーム単位のフィルタリング

ノイズサプレッサ

ノイズサプレッサは, 背景雑音を除去するアルゴリズム (あるいは, エフェクター) の 1 つです. 一般的に, スペクトルサブトラクションを利用して実装されることが多いです.

ノイズゲートの問題点

背景雑音を除去するアルゴリズムには, ノイズサプレッサだけではなく, ノイズゲート もあります. 非常に単純なアルゴリズムで, ある閾値以下の振幅は, 0 (無音) にしてしまうというものです. しかしながら, このアルゴリズムでは, 例えば, 楽器音にかさなった背景雑音を除去することはできません.

ノイズサプレッサは, この問題点を解決するアルゴリズムです.

実装

#include <cmath>
#include <complex>
#include "window_functions.h"
#include "fft.h"

enum {
  N = 1024
};

void NoiseSuppressor(
  double threshold,
  std::vector<double> &in,
  std::vector<double> &out,
  int fs,
  int length
) {
  int number_of_frame = (length - (N / 2)) / (N / 2);

  std::vector<std::complex<double>> x(N);
  std::vector<std::complex<double>> y(N);

  std::vector<double> w(N);
  std::vector<double> x_real(N);
  std::vector<double> x_imag(N);
  std::vector<double> A(N);
  std::vector<double> T(N);
  std::vector<double> y_real(N);
  std::vector<double> y_imag(N);

  hanning_window(w, N);

  for (int frame = 0; frame < number_of_frame; frame++) {
    int offset = (N / 2) * frame;

    for (int n = 0; n < N; n++) {
      x_real[n] = in[offset + n] * w[n];
      x_imag[n] = 0.0;

      x[n] = std::complex<double>(x_real[n], x_imag[n]);
    }

    FFT(x, N);

    for (int k = 0; k < N; k++) {
      A[k] = std::sqrt((x[k].real() * x[k].real()) + (x[k].imag() * x[k].imag()));

      if ((x[k].imag() != 0.0) && (x[k].real() != 0.0)) {
        T[k] = std::atan2(x[k].imag(), x[k].real());
      }
    }

    for (int k = 0; k < N; k++) {
      A[k] -= threshold;

      if (A[k] < 0.0) {
        A[k] = 0.0;
      }
    }

    for (int k = 0; k < N; k++) {
      y_real[k] = A[k] * cos(T[k]);
      y_imag[k] = A[k] * sin(T[k]);

      y[k] = std::complex<double>(y_real[k], y_imag[k]);
    }

    IFFT(y, N);

    for (int n = 0; n < N; n++) {
      out[offset + n] += y[n].real();
    }
  }
}

実行例

#include <iostream>
#include <cstdlib>
#include "wave.h"
#include "noise_suppressor.h"

int main(int argc, char **argv) {
  if (argc != 2) {
    std::cerr << "Require threshold" << std::endl;
    std::exit(EXIT_FAILURE);
  }

  STEREO_PCM pcm0, pcm1;

  double threshold = std::stod(argv[1]);

  wave_read(&pcm0, "stereo.wav");

  pcm1.fs     = pcm0.fs;
  pcm1.bits   = pcm0.bits;
  pcm1.length = pcm0.length;

  pcm1.sL.resize(pcm1.length);
  pcm1.sR.resize(pcm1.length);

  NoiseSuppressor(threshold, pcm0.sL, pcm1.sL, pcm1.fs, pcm1.length);
  NoiseSuppressor(threshold, pcm0.sR, pcm1.sR, pcm1.fs, pcm1.length);

  wave_write(&pcm1, "noise-suppressor.wav");
}

リファレンス

Share Comments
comments powered by Disqus