ari23のブログ

ari23のブログ

メーカ勤務エンジニアの技術ブログです

ディジタル信号処理 |周波数解析とFFT(解説)

今回のテーマはディジタル信号処理です🐜
その中でも、周波数解析またはFFTの要点について整理したいと思います。1

もちろん、継続的に加筆修正していきます!

周波数解析

周波数解析(スペクトル解析)とは、信号にどんな周波数の周期信号が含まれているのかを解析することである。

信号にフーリエ変換を施して、時間領域から周波数領域に変換して解析する。

フーリエ変換|連続時間領域と連続周波数領域

フーリエ変換は、信号を連続時間領域→連続周波数領域に変換する。 もちろん、フーリエ逆変換または逆フーリエ変換は、信号を連続周波数領域→連続時間領域に変換する。

フーリエ変換の定義式を以下に示す。

 \displaystyle
F(\Omega) = \int^{\infty}_{-\infty} f(t) e^{-j \Omega t} dt

フーリエ逆変換の定義式を以下に示す。

 \displaystyle
f(t) = \frac{1}{2 \pi} \int^{\infty}_{-\infty} F(\Omega) e^{ j \Omega t} d \Omega

ただし、 T を周期、 f を周波数としたとき、 \Omega = 2 \pi f = \frac{2 \pi}{T} が成り立つ。 この \Omega は(非正規化)角周波数を表す連続変数であり、 F(\Omega )は f(t) に含まれる角周波数 \Omega の振動成分の量を示す。 ちなみに、 |F(\Omega)|  f(t) の振幅スペクトル、 \angle F(\Omega)  f(t) の位相スペクトル、 |F(\Omega)|^2  f(t) パワースペクトルと呼ぶ。

これを線形・時不変システムで図示すると以下の通り。

FT 連続時間領域と連続周波数領域 模式図
FT 連続時間領域と連続周波数領域 模式図

なお、畳み込み式は次式で定義される。

 \displaystyle
y(t) = \int^{\infty}_{-\infty} x(\tau) h(t-\tau) d\tau = x(t) * h(t)

しかし、計算機(ディジタルな領域)で解析するには、離散化が必要であることがわかる。

離散時間フーリエ変換|離散時間領域と連続周波数領域

離散時間フーリエ変換(Dicrete-Time Fourier Transform; DTFT)は、信号を離散時間領域→連続周波数領域に変換する。 もちろん、離散時間フーリエ逆変換は、信号を連続周波数領域→離散時間領域に変換する。

離散時間フーリエ変換の定義式を以下に示す。

 \displaystyle
F(\omega) = \sum^{\infty}_{n=-\infty} f[n] e^{-j \omega n}

離散時間フーリエ逆変換の定義式を以下に示す。

 \displaystyle
f[t] = \frac{1}{2 \pi} \int^{\pi}_{-\pi} F(\omega) e^{ j \omega n} d \omega

この \omega は正規化角周波数を表す連続変数であり、 F(\omega)  f[n] に含まれる正規化角周波数 \omega の振動成分の量を表す。 また、 F(\omega) は周期 2\pi の周期関数で、通常 -\pi \lt \omega \lt \pi の範囲で考える。

なお、サンプリング周期を T_s としたとき、非正規化角周波数 \Omega と正規化角周波数 \omega の間に以下の式が成り立つ。

 \displaystyle
\omega = \Omega T_s

これを線形・時不変システムで図示すると以下の通り。

DTFT 離散時間領域と連続周波数領域 模式図
DTFT 離散時間領域と連続周波数領域 模式図

しかし、計算機で解析するには、無限和の解消と周波数領域の離散化が必要であることがわかる。

離散フーリエ変換|離散時間領域と離散周波数領域

離散フーリエ変換(Discrete Fourier Transform; DFT)は、信号を離散時間領域→離散周波数領域に変換する。 もちろん、離散フーリエ逆変換は、信号を離散周波数領域→離散時間領域に変換する。

離散フーリエ変換の定義式を以下に示す。

 \displaystyle
F[k] = \sum^{N-1}_{n=0} f[n] e^{-j \frac{2\pi}{N} kn}

離散フーリエ逆変換の定義式を以下に示す。

 \displaystyle
f[n] = \frac{1}{N} \sum^{N-1}_{k=0} F[k] e^{j \frac{2\pi}{N} kn}

 f[n] は離散時間 n (整数)で定義された周期 N の周期関数であり、 F[k] も整数 k で定義される周期 N の周期関数である。  k は周波数のインデックスを表す整数で、 F[k]  f[n] に含まれる正規化角周波数 \frac{2\pi k}{N} [rad]の振動成分の量を表す。

これを線形・時不変システムで図示すると以下の通り。

DFT 離散時間領域と離散周波数領域 模式図
DFT 離散時間領域と離散周波数領域 模式図

なお、畳み込み式は次式で定義される。

 \displaystyle
y[n] = \sum^{\infty}_{k=-\infty} x[k] h[n-k] = x[t] * h[t]

高速フーリエ変換(FFT)

結局、周波数解析で求めたいのは F(\Omega) または |F(\Omega)| であるが、計算機で解析するために様々な近似を経て、 F[k] を計算して解析を行う。

しかし、DFTでも計算時間が O(N^{2}) であり、現実的ではない。

そこで一般には、この計算時間が O(N \log_2 N) になる高速フーリエ変換(Fast Fourier Transform; FFT)を利用する。

このFFTでは、そのアルゴリズムや計算機の性質から N=1024 など2の乗数の点数で計算を行うと高速にできる。

アルゴリズムの詳細はここでは議論しない。

フーリエ変換のまとめ

ここまで3種のフーリエ変換について概要を述べたが、これらにフーリエ級数展開を加えて、離散性と周期性に注目して以下に整理する。

フーリエ変換まとめ
フーリエ変換まとめ

なお、フーリエ級数の定義式を以下に示す。ただし、ここでは議論しない。

 \displaystyle
F_k = \frac{1}{T_0} \int^{T_0 /2}_{-T_0 /2} f(t) e^{-j \Omega_0 kt} dt
 \displaystyle
f(t) = \sum^{\infty}_{k=-\infty} F_k e^{j \Omega_0 kt}

実装時のポイント

Pythonの信号処理ライブラリで有名なSciPyを使う場合、主に以下の2つの関数を使います。

  • scipy.fftpack.fft
    中身はDFTのFFTです。戻り値は F[k] ですので、複素数であることに注意しましょう。
    また、一般に周波数解析でy軸にとるのは |F[k]| なので、適宜パワースペクトルに変換します。
  • scipy.fftpack.fftfreq
    scipy.fftpack.fftの戻り値に対応する、離散化された周波数を計算する関数です。
    必要に応じて正規化各周波数→非正規化角周波数に変換します。

NumPyにも numpy.fft.fftとnumpy.fft.fftfreqという同じ関数が用意されていますが、SciPyの方が速度が速いそうです。

付録|複素数 オイラーの公式

複素数の理解のため、オイラーの公式を下記に示す。

 \displaystyle
e^{i\theta} = \cos\theta + i \sin\theta

なお、オイラーの等式は以下の通り2

 \displaystyle
e^{i\pi} + 1 = 0

詳細は下記を参照のこと。

参考文献

おわりに

数学的な議論や定義式の導出は省き、要点だけをまとめました。
(大変だった。。。)

これを見れば、なぜFFTの戻り値が複素数なのか、なぜfreqの計算結果が -\omega から \omega の値なのかがわかると思います。

参考になれば幸いです(^^)


  1. FFTで大事なのはあくまでも近似式であることです。ほんとはフーリエ変換やりたいけど、計算機では実現が難しいからFFTをしているということ。

  2. 恥ずかしながらこのブログをきっかけに初めて知ったんですが、ふつくしい