Fourier Transform

フーリエ変換

周波数という視点 (ドメイン) から, 音や光の特徴を調べる方法がフーリエ変換です. 音の場合であれば, 基本周波数と倍音の比率などを分析することができます.

フーリエ変換

フーリエ変換は, $-\infty$ ~ $\infty$ の時刻まで観測したアナログ信号の周波数特性を調べるための数学的処理で, フーリエ変換と, その逆変換である, 逆フーリエ変換は以下のように定義されています.

$ \begin{eqnarray} X(f)=\int_{-\infty}^{\infty}x(t)\exp(-j2{\pi}ft)dt \quad (-\infty \leqq f \leqq \infty) \end{eqnarray} $

$ \begin{eqnarray} x(t)=\int_{-\infty}^{\infty}X(t)\exp(j2{\pi}ft)df \quad (-\infty \leqq t \leqq \infty) \end{eqnarray} $

$x(t)$ は, 時間 $t$ を変数とするアナログ信号, $X(f)$ は周波数 $f$ を変数とする $x(t)$ の周波数特性を表しています.

離散フーリエ変換

アナログ信号のフーリエ変換では, 無限の長さの連続信号となりますが, コンピュータでは無限大連続した信号をあつかうことができません. したがって, コンピュータでフーリエ変換を実装する場合には, 離散フーリエ変換 (DFT : Discrete Fourier Transform) と, その逆変換である逆離散フーリエ変換 (IDFT : Inverse Discrete Fourier Transform) を利用します.

$ \begin{eqnarray} t=nt_{s} \quad (1 \leqq n \leqq N-1) \end{eqnarray} $

$ \begin{eqnarray} f=\frac{kf_{s}}{N} \quad (f \leqq k \leqq N - 1) \end{eqnarray} $

と定義できるので,

$ \begin{eqnarray} X(k)=\sum_{n=0}^{N-1}x(n)\exp\left(\frac{-j2{\pi}kn}{N}\right) \end{eqnarray} $

$ \begin{eqnarray} x(n)=\frac{1}{N}\sum_{n=0}^{N-1}X(k)\exp\left(\frac{j2{\pi}kn}{N}\right) \end{eqnarray} $

音の場合, $x(n)$ は実数となります. 一方で, $X(k)$ は複素数となります. $X(k)$ を極座標形式で表すと, 以下のようになります.

$ \begin{eqnarray} X(k)=A(k)\exp(j{\theta}(k)) \end{eqnarray} $

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

ここで, $X(k)$ の逆離散フーリエ変換を考えます.

$ \begin{eqnarray} x(n)=\frac{1}{N}\sum_{n=0}^{N-1}A(k)exp\left(\frac{-j2{\pi}kn}{N}+j\theta(k)\right) \quad (0 \leqq n \leqq N - 1) \end{eqnarray} $

オイラーの公式を適用すると,

$ \begin{eqnarray} x(n)=\frac{1}{N}\sum_{n=0}^{N-1}A(k)\left[\cos\left(\frac{2{\pi}kn}{N}+\theta(k)\right)+j\sin\left(\frac{2{\pi}kn}{N}+\theta(k)\right)\right] \quad (0 \leqq n \leqq N - 1) \end{eqnarray} $

音の場合, $x(n)$ は, 虚数になるので, 虚数部を消去すると, 以下のようになります.

$ \begin{eqnarray} x(n)=\frac{1}{N}\sum_{n=0}^{N-1}A(k)\cos\left(\frac{2{\pi}kn}{N}+\theta(k)\right) \quad (0 \leqq n \leqq N - 1) \end{eqnarray} $

この式は, $N$ 個の $\cos$ 波の重ねあわせによって $x(n)$ を合成できることを表しています. したがって, 離散フーリエ変換によって計算される周波数特性 (スペクトル) は, $\cos$ 波の振幅と位相の情報です ($A(k)$ は, 振幅スペクトル, $\theta(k)$は, 位相スペクトル です).

周波数特性は, $k=\frac{N}{2}$ を中心に, 振幅スペクトルは線対称, 位相スペクトルは点対称になります. さらに, $k=\frac{N}{2}$ は, $\frac{f_{s}}{2}$, つまり, サンプリング周波数の $\frac{1}{2}$ に対応し, サンプリング定理にも関係しています.

人間の聴覚は, 位相の違いに鈍感 なので, 音の周波数特性といえば振幅スペクトルを意味していることがほとんどです (ただし, ホワイトノイズとインパルス音のように位相の違いを分析する必要がある場合もあります).

窓関数

離散フーリエ変換では, $x(n)$, $x(k)$ が, 離散フーリエ変換のサイズ $N$ の周期信号という前提があります. つまり, $N$ が $sin$ 波の整数倍となっていれば, 離散フーリエ変換によって求められる周波数特性は, 基本周波数をのみを含んだものになります.

ところが, そうでない場合, 実際の音にはない周波数成分が含まれてしまいます.

この問題を緩和するのが, 窓関数です. 窓関数を乗算することで, 端部の不連続点をあいまいにすることができ, 周波数成分の広がりをある程度おさえることができます.

窓関数にはいくつか種類がありますが, 代表的な窓関数として, ハニング窓があり, 以下のように定義されています(ほかには, ブラックマン窓などがあります).

($N$ が偶数の場合)

$ w(n)=\begin{cases} \begin{eqnarray} &0.5-0.5\cos\left(\frac{2{\pi}n}{N}\right)& \quad (0 \leqq n \leqq N - 1) \\
&0& \quad (otherwise) \\
\end{eqnarray} \end{cases} $

($N$ が奇数の場合)

$ w(n)=\begin{cases} \begin{eqnarray} &0.5-0.5\cos\left(\frac{2{\pi}(n+0.5)}{N}\right)& \quad (0 \leqq n \leqq N - 1) \\
&0& \quad (otherwise) \\
\end{eqnarray} \end{cases} $

実装

C++ で実装した, 離散フーリエ変換の関数です.

#include <vector>
#include <cmath>
#include <complex>

void DFT(std::vector<std::complex<double>> &x, std::vector<std::complex<double>> &X, int N) {
  for (int n = 0; n < N; n++) {
    for (int k = 0; k < N; k++) {
      double real =  std::cos((2 * M_PI * k * n) / N);
      double imag = -std::sin((2 * M_PI * k * n) / N);

      double real_sum = X[k].real() + ((real * x[n].real()) + (imag * x[n].imag()));
      double imag_sum = X[k].imag() + ((real * x[n].imag()) + (imag * x[n].real()));

      X[k] = std::complex<double>(real_sum, imag_sum);
    }
  }
}

ハニング窓

#include <vector>
#include <cmath>

void hanning_window(std::vector<double> &w, int N) {
  if ((N % 2) == 0) {
    for (int n = 0; n < N; n++) {
      w[n] = 0.5 - (0.5 * std::cos((2.0 * M_PI * n) / N));
    }
  } else {
    for (int n = 0; n < N; n++) {
      w[n] = 0.5 - (0.5 * std::cos((2.0 * M_PI * (n + 0.5)) / N));
    }
  }
}

実行例

#include <iostream>
#include <vector>
#include <cmath>
#include <complex>
#include "wave.h"
#include "window_functions.h"
#include "dft.h"

enum {
  N = 8
};

int main(void) {
  MONO_PCM pcm;

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

  wave_read(&pcm, "monoral.wav");

  hanning_window(w, N);

  for (int n = 0; n < N; n++) {
    x[n] = std::complex<double>((w[n] * pcm.s[n]), 0.0);
  }

  DFT(x, X, N);

  for (int k = 0; k < N; k++) {
    std::cout << "[" << k << "]" << " " <<  X[k].real() << " + j" << std::fixed << X[k].imag() << std::endl;
  }
}

離散コサイン変換

フーリエ変換を利用することで, 周波数解析が可能になります. しかし, フーリエ変換では, その演算結果が複素数になってしまうというデメリットがあります. 実数である時間領域の信号を実数の周波数信号に変換できれば, コンピュータで処理することが容易になります. このセクションでは, MPEG などで利用されている離散コサイン変換 (Discrete Cosine Transform) について解説します.

$ \begin{eqnarray} X(k)=\sqrt{\frac{2}{N}}a_{K}\sum_{n=0}^{N-1}x(n)\cos\left(\frac{\left(2n+1\right)k{\pi}}{2N}\right) \end{eqnarray} $

$ \begin{eqnarray} a_{k}= \begin{cases} &1&\quad(k=1,2,\cdots,N-1) \\
&\frac{1}{\sqrt{2}}&\quad(k=0)\\
\end{cases} \end{eqnarray} $

また, その逆変換である, 逆離散コサイン変換は以下のように定義されます.

$ \begin{eqnarray} x(n)=\sqrt{\frac{2}{N}}\sum_{k=0}^{N-1}a_{K}X(k)\cos\left(\frac{\left(2n+1\right)k{\pi}}{2N}\right) \end{eqnarray} $

$ \begin{eqnarray} a_{k}= \begin{cases} &1&\quad(k=1,2,\cdots,N-1) \\
&\frac{1}{\sqrt{2}}&\quad(k=0)\\
\end{cases} \end{eqnarray} $

変形離散コサイン変換

離散コサイン変換では, $N$ 個ずつデータを取得して, 周波数変換, 信号処理, そして, 逆変換を実行します. しかしながら, それでは, データのまとまりに不連続な点が発生してしまいます.

そこで, 信号を窓関数を利用してオーバラップ (隣接するデータのまとまりを重ねて連結していく処理) をしながら, 信号を取得して, 周波数解析をする処理として, 変形離散コサイン変換 (修正離散コサイン変換 (Modified Discrete Cosine Transform)) が利用されます.

オーバーラップ (赤い矩形の重なった部分がオーバーラップしている部分

$ \begin{eqnarray} X(k)&=&\sum_{n=0}^{2N-1}x(n)\cos\left(\frac{\left(2n+N+1\right)\left(2k+1\right){\pi}}{4N}\right) \\
&=&\sum_{n=0}^{2N-1}x(n)\cos\left(\frac{4\left(n+\frac{N}{2}+\frac{1}{2}\right)\left(k+\frac{1}{2}\right){\pi}}{4N}\right) \\
&=&\sum_{n=0}^{2N-1}x(n)\cos\left(\frac{{\pi}}{N}\left(n+\frac{N}{2}+\frac{1}{2}\right)\left(k+\frac{1}{2}\right)\right) \quad (k=0,1,\cdots,N-1) \\
\end{eqnarray} $

また, 逆変形離散コサイン変換は, 以下のように定義されます.

$ \begin{eqnarray} x(n)&=&\frac{1}{N}\sum_{k=0}^{N-1}X(k)\cos\left(\frac{\left(2n+N+1\right)\left(2k+1\right){\pi}}{4N}\right) \\
&=&\frac{1}{N}\sum_{k=0}^{N-1}X(k)\cos\left(\frac{4\left(n+\frac{N}{2}+\frac{1}{2}\right)\left(k+\frac{1}{2}\right){\pi}}{4N}\right) \\
&=&\frac{1}{N}\sum_{k=0}^{N-1}X(k)\cos\left(\frac{{\pi}}{N}\left(n+\frac{N}{2}+\frac{1}{2}\right)\left(k+\frac{1}{2}\right)\right) \quad (n=0,1,\cdots,2N-1) \\
\end{eqnarray} $

実装

C++ による, 変形離散コサイン変換の実装例です.

#include <vector>
#include <cmath>

void MDCT(std::vector<double> &x, std::vector<double> &X, int N) {
  for (int k = 0; k < N; k++) {
    for (int n = 0; n <= ((2 * N) - 1); n++) {
      X[k] += x[n] * std::cos((M_PI / N) * (n + (N / 2.0) + (1.0 / 2.0)) * (k + (1.0 / 2.0)));
    }
  }
}

void IMDCT(std::vector<double> &x, std::vector<double> &X, int N) {
  for (int n = 0; n <= ((2 * N) - 1); n++) {
    for (int k = 0; k < N; k++) {
      x[n] += X[k] * std::cos((M_PI / N) * (n + (N / 2.0) + (1.0 / 2.0)) * (k + (1.0 / 2.0)));
    }

    x[n] = x[n] / N;
  }
}

実行例

#include <iostream>
#include <vector>
#include <cmath>
#include "../wave/wave.h"
#include "mdct.h"

enum {
  N = 8
};

int main() {
  MONO_PCM pcm;

  std::vector<double> x(2 * N);
  std::vector<double> X(N);

  WAVE::wave_read(&pcm, "sample.wav");

  std::cout << "x(n)" << std::endl;

  for (int n = 0; n < (2 * N); n++) {
    x[n] = pcm.s[n];
    std::cout << "[" << n << "]" << " " << x[n] << std::endl;
  }

  MDCT(x, X, N);

  std::cout << "X(k)" << std::endl;

  for (int k = 0; k < N; k++) {
    std::cout << "[" << k << "]" << " " << X[k] << std::endl;
  }

  x.assign((2 * N), 0);

  IMDCT(x, X, N);

  std::cout << "x(n)" << std::endl;

  for (int n = 0; n < (2 * N); n++) {
    std::cout << "[" << n << "]" << " " << x[n] << std::endl;
  }
}

リファレンス

Share Comments
comments powered by Disqus