概要

フーリエ変換を離散化したものが離散フーリエ変換(DFT)であり、
計算機ではこのDFTを高速化した高速フーリエ変換(FFT)を用いて、
フーリエ変換をしています。

DFTの練習問題としてよく出てくるのが，
(1)方形波
(2)正弦波
の2つです．
正直フーリ変換と違ってDFTは手計算できるものがかなり限られてくるので
この2つを押さえておくだけでもだいぶ違うかなと思います．

今回は(1)方形波を扱います．(2)正弦波は次回扱おうと思います．

[【信号処理】Pythonを使って三角波をフーリエ変換してみた！]
(https://qiita.com/trami/items/f228b476e12578656560)

DFTを実際に手計算で計算し，それをPythonでFFTしたものと一致するかを確かめてみたいと思います．

動作環境

Windows10(64bit)
Python 3.7.2

DFTの復習

入力信号をx(n)とすると，そのDFTであるX(k)は以下のようになります．ただし，$W_N = exp(-j\frac{2\pi}{N})$.


X(k) = \sum_{n=0}^{N-1}{x(n)W_{N}^{kn}}

問題


x(n) = \begin{cases}
1 & (0 \le n \le 3)\\
0 & (4 \le n \le 7)
\end{cases}

のDFTを求めなさい．ただし，サンプル数$N=8$とする．

解答

それでは上記の問題を，まずは手計算で計算していきます．

まずは，回転子を導入する．


W_8 = exp(-j\frac{2\pi}{8})

DFTの定義より


\begin{align}
X(k) & = \sum_{n=0}^{7}x(n)W_8^{kn}\\
& = \sum_{n=0}^{3}W_8^{kn}\\
& = W_8^0 + W_8^k + W_8^{2k} + W_8^{3k}
\end{align}

$k\neq 0$ のとき


\begin{align}
X(k) & = \frac{1-W_8^{4k}}{1-W_8^k}\\
& = \frac{W_8^{2k}}{W_8^{\frac{k}{2}}}\frac{W_8^{-2k}-W_8^{2k}}{W_8^{-\frac{k}{2}}-W_8^{\frac{k}{2}}}\\
& = exp(-j\frac{3\pi}{8}k)\frac{sin(\frac{k\pi}{2})}{sin(\frac{k\pi}{8})}
\end{align}

$k = 0$ のとき


X(k) = \sum_{n=0}^{3}1 = 4

実装

では，早速実装していきましょう．

square_wave.py


import numpy as np
import matplotlib.pyplot as plt

N = 64            
n = np.array(range(0,N))
f = np.zeros(n.shape)
f[0:4] = 1
F = np.fft.fft(f)

fig = plt.figure()
ax1 = fig.add_subplot(211)
ax2 = fig.add_subplot(212)

ax1.stem(n, f)
ax2.stem(n, np.abs(F))
ax1.set_xlabel('t')
ax2.set_xlabel('k')
ax1.set_ylabel('f')
ax2.set_ylabel('abs(F)')
fig.tight_layout()
plt.show()

解説

特に難しいところはないと思います．
高速フーリエ変換はnp.fft.fft()でできます．
グラフ描画は離散データなのでstemで描画するためnp.stem()となります．

実行結果

以上のプログラムを実行すると以下のような画像が得られます．

$k = 7, 5$の成分は$k = -1, -3$の成分と解釈することもできますね．
ただ$N=8$だと計算は手計算でもできますが，なかなかFFT後の画像から定性的に理解するのが難しいと思います．
サンプル数を増やして$N=64$として実行した結果が以下になります．

これを見ると，(ふつうの)フーリ変換で得られる，方形波をフーリ変換すると，
sinc関数になるということが若干わかるのではないでしょうか？

まとめ

今回は，DFTで頻出の方形波のDFTに挑戦してみました．いかがだったでしょうか？
計算に関してはほぼお決まりのパターンなので，たくさん練習をして慣れていくといいと思います．
サンプル数を変えたり，入力信号を変えたりして遊んでみると面白いのではないかなと思います．
次回は，DFTで方形波と並んで頻出の正弦波について扱っていきたいと思います．
今回のプログラムに関して何か改善点や建設的な意見等あればコメントしてもらえると助かります！