フーリエ変換(1)　

はじめに

今回は，情報処理の応用の一つとして「信号処理」を扱う． 信号処理とは，実験などで得られた何らかの波形データを（数学に基づいた）方法で加工し，その特徴を抽出したり，雑音を軽減したり，必要な情報を抽出したりする手法の総称である． 広い意味では，2次元データである画像などの解析（画像処理と呼ばれ，顔認識，モザイク処理，色調変換や画像圧縮など）も信号処理に含まれる． ここでは簡単に1次元データの処理を考えよう．

ある波形データがあったとして，この波形から何か情報を抽出したり加工したいとする． 特に，周期的に変動する現象の解析は，工学のさまざまな分野において利用されている．
例えば，自動車や鉄道，船舶などの乗り物や，橋脚，建物，内燃機関など振動をなるべく抑制したい場合の解析に用いられる． 逆に楽器やスピーカー，通信アンテナなど所望の周波数でうまく振動・共振させたい場合も多くあり，その目的はさまざまである．

この例のように，実際の波形には何らかの周期性があることが予想されるが，単純な正弦波などには程遠く，それらが組み合わさった複雑な波形であるため，それ自体を眺めていてもよく分からない． そこで，信号にどのような周波数（または振動数）成分が，どのくらいの大きさ（振幅，またはその2乗のパワー）で含まれているかを数値的に解析することによって，その性質を把握することが出来る．

フーリエ変換とは

周期的な波形は，様々な周波数の正弦波の和で表すことができると考えられる．

いま，ある区間 $[0, 2\pi]$ において時間変動する関数 $f(t)$ を，正弦関数，余弦関数（$\sin, \cos$）の和で表すことを考える．

式で書くと，

\begin{align} f(t) &= \sum_{n=0}^\infty (a_n \cos nt + b_n \sin nt) \\ &= a_0 + \sum_{n=1}^\infty (a_n \cos nt + b_n \sin nt) \label{eqn:exp} \end{align}

とあらわせる．

この式は，いろいろな周波数の $\sin, \cos$ 波をある振幅で足してゆけば，任意の関数 $f$ を表すことができることを示している． 言い換えれば，任意の関数は，様々な周波数 $\omega$ にたいする振幅 $a_n, b_n$ に分解することができる．
では， $a_n, b_n$ の値は，どのような計算をすれば良いだろうか．

一つの方法として，三角関数の直交性を利用して，関数 $f$ と $\sin, \cos$ との内積を計算すれば， $a_n, b_n$ を取り出すことができる．
この性質をうまく利用して，各周波数に対応する振幅を求める方法がフーリエ変換である．

（連続）フーリエ変換

このような考え方を一般化して，任意の関数を，様々な周波数の三角関数 $\sin, \cos$ の和で表すことができることが知られており，フーリエ変換と呼ばれる．
フーリエ変換（Fourier Transform)がどのような変換で，どのような性質かを理解しよう．
（フーリエ変換は，フランスの数学者 Joseph Fourier により，熱伝導を解析するために考案された手法である．）

時間 $t$ の関数 $f (t)$ があるとして， この関数のフーリエ変換 $F(\omega)$ は以下のように定義される．
ここでは，オイラーの公式 $e^{j\omega t}=\cos\omega t + j\sin\omega t$ を使用している．
$j$ は虚数単位で，$j^2=-1$ である．
（工学の分野では，虚数単位に $i$ ではなく $j$ が使われることがある．$i$ は電流で使用されることが多いため？）

\begin{align} F(\omega) &= \int_{-\infty}^\infty f(t)\ \left(\cos(-\omega t) + j\sin(-\omega t)\right)\ dt \\ &= \int_{-\infty}^\infty f(t)\ e^{-j\omega t}\ dt \end{align} （注：複素数の指数関数は周期関数．）

$\omega~\mathrm{[rad/s]}$ は角周波数であり，周波数を $f~\mathrm{[Hz]}$ とすると，$\omega=2\pi f$の関係がある．(この $f$ は関数でなく周波数)

フーリエ変換は，時刻 $t$ の関数である $f$ を，周波数 $\omega$ の関数である $F$ に変換する．
$\omega$ の値を与えれば，$f$ 中の周波数 $\omega$ に相当する $\sin, \cos$ 成分の振幅がそれぞれ実部，虚部として得られる．

離散フーリエ変換

サンプリングについて

圧力や電圧など，連続的に変化する信号を計測により一定の時間間隔で取得し，離散的なデータ列に変換することを「サンプリング」や「標本化」と言う．
例えば，音声データや温度変化などの時系列データは，（時間軸）一次元の配列のデータとなる． また，スマホやディジタルカメラなどで撮影した画像は，（空間軸）二次元の配列のデータとなる．
この測定の時間刻みのことをサンプリング間隔，またはサンプリング周期と呼び，単位は時間(sec, min等)となる．
サンプリング周期の逆数をサンプリング周波数と呼び，単位は $\mathrm{Hz}$ または $\mathrm{s^{-1}}$ である．
サンプリング周期は被測定量の変化よりも短く（つまり，対象が変化するよりも短時間で）する必要がある．

例１：24時間の外気温の変化をサンプリングするには　数分〜30分程度の時間間隔で十分であるが，エアコンや冷蔵庫の制御では 1分またはそれ以下の短い時間間隔で温度をサンプリングする必要がある．
例２：音声録音では一般的に 44.1kHz（=44100 Hz，つまり1秒間に44100回のサンプリング）が使用されており，これ以外にも 48kHz, 96kHz（ハイレゾリューション）などが利用されている．
例３：動画撮影のサンプリング周波数(Frames per second, fpsと呼ばれる)は一般に 24 Hz, 30 Hz 程度であるが，最近のスマートフォンの高速モードでは 240 Hz 程度，さらに産業用の高速度カメラでは数千Hz〜10億Hz以上のものがある．

サンプリング速度について詳しく知りたい場合は，標本化定理，ナイキスト周波数を参照．

いま，サンプリングする時間間隔を $\Delta t$ とすれば，データ配列の 0 番目には $f(0)$，1番目には $f(\Delta t)$，2番目には $f(2\Delta t)$・・・，といった計測データが順次入る．
コンピュータによるデータ処理では，必ず有限区間の離散化されたデータを対象とする．

離散フーリエ変換の計算

フーリエ変換の式は，連続である関数に対して定義されるが，一方で，計算機上で表現できるデータは有限長，かつ，サンプリングされた有限個の離散データである．
離散的なデータに対してのフーリエ変換は 離散フーリエ変換(Discrete Fourier Transform, DFT)と呼ばれ，コンピュータで計算できるのはこのDFTである．

ここでは，連続関数をサンプリングした有限個の離散データをC言語の配列風に d[t] とする．
（f[t] とすると周波数 f とまぎらわしいので，d[t]とした．以降，f は周波数を表すものとする．）

離散データ d[0], d[1], ... d[N-1] にたいする離散フーリエ変換は，以下のようにあらわされる．
$\omega = 2\pi f$ 置き換えているだけで，以下の二つの式は同一である．
\begin{align} F[\omega] &= \sum_{t=0}^{N-1} \mathrm{d[t]}\ e^{-j\ (\omega\frac{t}{N}) }\ \Delta t \\ F[f] &= \sum_{t=0}^{N-1} \mathrm{d[t]}\ e^{-j\ (2\pi f\frac{t}{N}) }\ \Delta t \label{eqn:dft} \end{align}

$F[\omega]，F[f]$が，フーリエ変換後の関数（有限個の複素数配列）である．

右辺の指数関数の部分$e^{-j\ 2\pi f\frac{t}{N}}$は，$f = 1,2,3,...$ ($\omega = 2\pi, 4\pi, 6\pi$, ...) と増加するにつれ，周波数が増加する周期関数となる．
即ち，離散フーリエ変換は，いろいろな周波数　$f$ = 1, 2, 3, ... の余弦波・正弦波と，元の波形との内積により計算できる．

例えば，$f=1$ ($\omega=2\pi$) とすれば，元波形と1周期分の余弦波（実部）＋正弦波（虚部）との内積が計算でき，その結果の実部・虚部がそれぞれ式$\eqref{eqn:exp}$ $a_1$，$b_1$ （の定数倍）となる．
同様に，$f = 2,3$ ($\omega=4\pi，6\pi$) とすれば，それぞれ2周期分，3周期分の余弦波＋正弦波と元波形との内積が計算でき，それぞれ $a_2, b_2$，$a_3, b_3$ が計算できる．

プログラム中での複素数の扱い方

C言語の組み込み型には，複素数型は存在しないため，C++ 標準テンプレートライブラリ(STL)の複素数クラス complex 型を用いる．

複素数型の定義には，実部・虚部を表す変数の型をあわせて指定する．（テンプレート機能）
このクラスライブラリには，実部・虚部を格納するメンバ変数だけでなく，複素数に関連する演算（四則演算，絶対値，指数関数など）が定義されているため，簡単に利用できとても便利である．

複素数型変数complexの使用例

#include <iostream>
#include <complex>
using namespace std;

int main(void)
{
    complex<double> z1, z2;  //  実部と虚部が double 型である複素数型変数の定義

    // 実部・虚部を一括で代入
    z1 = complex<double>(0.2, 1.0);

    // 実部・虚部をそれぞれ個別に代入
    z2.real(1.6);
    z2.imag(-2.3);
  
    // 複素数を表示
    cout << "z1 = " << z1 << endl;
    cout << "z2 = " << z2 << endl;

    // 実部，虚部を別々に表示
    cout << "Re{z1} = " << z1.real() << endl;
    cout << "Im{z1} = " << z1.imag() << endl;

    // 四則演算
    cout << "z1 + z2 = " << z1 + z2 << endl;
    cout << "z1 - z2 = " << z1 - z2 << endl;
    cout << "z1 * z2 = " << z1 * z2 << endl;
    cout << "z1 / z2 = " << z1 / z2 << endl;

    // 共役複素数 conj()
    cout << "conj(z1) = " << conj(z1) << endl;

    // 絶対値 abs()
    cout << "abs(z1)  = " << abs(z1) << endl;

    // ノルム(絶対値の2乗) norm()
    cout << "norm(z1) = " << norm(z1) << endl;

    // 偏角 arg()
    cout << "arg(z1)  = " << arg(z1) << endl;

    // 指数関数 exp()
    cout << "exp(z1)  = " << exp(z1) << endl;

    // 累乗  pow()
    int n = 10;
    cout << "pow(z1, " << n << ")  = " << pow(z1, n) << endl;
}

次に，複素数型の「配列」を使用する練習をしてみよう．

// 複素数型の配列の宣言と使用例

#include <complex>
#include <stdio.h>
#include <math.h>
using namespace std;

const double pi = 2.0*asin(1.0);

int main(void)
{
    const int N = 10;
    complex<double> z;     // 複素数型の変数

    complex<double> y[N];     // 複素数型の配列.
//   static complex<double> y[N];     // 複素数型の配列． Nが数千以上の場合，staticをつけてヒープメモリを使用する必要がある．

    int i;
    for(i=0; i<N; i++) {

        // ヒント：2.0*pi*i/N は，i が 0 から N まで変化すると，0 から 2π まで変化する．

        z = complex<double> (0.0, 2.0*pi*i/N);    // 1周期分の純虚数

        y[i] = exp(z);    // 複素数の指数関数

        printf("%lf, %lf\n", y[i].real(), y[i].imag());   // 表示
    }
    return 0;
}

typedef complex<double> dcmpx;

int main()
{
    dcmpx z1, z2;      // 複素数型変数の宣言
    z1 = dcomx (1.0, 2.0);   
    ...
}

using dcmpx = complex<double>;

int main()
{
    dcmpx z1, z2;      /* 複素数型変数の宣言 */
    z1 = dcomx (1.0, 2.0);   
    ...
}

練習問題

上記のサンプルソースコードを用いて，純虚数 $z = 0j \sim 2\pi j$ （$j$は虚数単位）までの区間を100分割し， その指数関数 $y = \exp(z)$ の値を求め，図示したい．
複素数用のライブラリの exp() 関数を利用し，実部，虚部をそれぞれカンマ区切りファイル(csvファイル)に出力し，

1. 得られた値の実部と虚部を線グラフとして描き，確認しよう．
2. 次に，実部を横軸，虚部を縦軸として複素平面上に 散布図として描き，形状を確認しよう．

実行例（.csvファイルの中身）

real part, imaginary part
1.000000, 0.000000
0.998027, 0.062791
0.992115, 0.125333
...
0.992115, -0.125333
0.998027, -0.062791

#include<complex>

//  内積を計算する関数（複素数ver）
complex<double> inner_product(complex<double> a[], complex<double> b[], int N)
{
    ...
}

#include <stdio.h>
#include <complex>
using namespace std;

const double pi = 2.0*asin(1.0);

//  内積を計算する関数（複素数ver）
complex<double> inner_product(complex<double> a[], complex<double> b[], int N)
{
    ...
}

int main(void)
{
    const int N = 100;
    int i;
    complex<double> a[N], b[N];
  
    //  ここを色々変えてみる．scanfでKBD入力にするとよい．
    double f1 = 3;
    double f2 = 4;

    //  initialize array
    for(i=0; i<N; i++) {
        a[i] = complex<double>( sin( f1 * 2.0*pi*i/N), 0.0);  //  実部にf1周期分の正弦波，虚部は0
        b[i] = complex<double>( sin( f2 * 2.0*pi*i/N), 0.0);  //  実部にf2周期分の正弦波，虚部は0
    }

    complex<double> ret = inner_product(a, b, N);

    printf("%lf, %lf\n", ret.real(), ret.imag());   // 表示
}

実習課題

  実行例．値は正しいとは限らない．
  f=0 : F = 1.32123 + 2.3313 i
  f=1 : F =  0.32123 - 1.0301 i
  f=2 : F = -3.91284 + 0.3910 i
  ...