ちょっとおもしろいフーリエ変換の使い方

　こんにちは、(株)日立製作所  Lumada Data Science Lab. の中川です。普段は人工知能を制御に適用する研究に従事しています。近年、機械学習が注目される中、機械学習理論および機械学習を使った技術開発環境は急速に進歩すると共に、多くの方がデータサイエンスに関わるようになってきました。すでにデータサイエンスに携わっている方や、これからデータサイエンスに関わってみようと思っている方の中で、理論の大切さをあらためて知りたいあるいは感じたいという方がいらっしゃいましたら、それをできるだけわかりやすく伝えられたら、という思いから基本的な内容で記事を書きます。今回は、ちょっとおもしろいフーリエ解析の性質についてご紹介します。



# 1. はじめに

　（1）式は、離散フーリエ変換を基にした、ある周波数の振幅と比例関係にある値を求める式です。


$$
\begin{align}
&A = \sqrt{(x(j) - x(j+2))^2 +(x(j+3)-x(j+1))^2}{\quad}・・・(1)\\\
&x(j):時系列データのj番目の値\\\
&j:0,1,・・・,N\hspace{0.5cm}データ番号\\
\end{align}
$$

　フーリエ変換あるいは離散フーリエ変換と聞くと、複雑な式を想像する人にとっては、比較的簡単な式に見えないでしょうか。この式で、うまくサンプリング周期を設定できれば、任意の周波数の振幅が検出できることが離散フーリエ変換から導けます。以下、詳細をご紹介します。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/641763/e684cf71-1ade-0148-bb10-5fc710a5aee7.png)


# 2. 離散フーリエ変換

　離散フーリエ変換とは、値がとびとびになっているときに用いるフーリエ変換です。考え方は、フーリエ変換と同じですが、実データは一般に離散値ですので、実データにフーリエ変換を適用する場合は、離散フーリエ変換を用いることになります。（2）式は、離散フーリエ変換の式です。
$$
\begin{align}
F(n)&=\frac{1}{N}\sum_{j=0}^{N-1}f(j)\hspace{0.05cm}e\hspace{0.05cm}^{-i2\pi{n\frac{j}{N}}}\hspace{2.1cm}・・・(2)\\\
n&:0,1,・・・,N\hspace{0.5cm}周波数番号\\\
j&:0,1,・・・,N\hspace{0.5cm}データ番号\\
\end{align}
$$

　（3）式と（4）式は、複素フーリエ級数を表した式ですが、（4）式の$\omega_0t$を離散値である$2{\pi}j/N$に置き換えて、値の個数が無限から有限（$N$個）とすることで（2）式が得られます。

$$
\begin{align}
f(t)&=\sum_{n=-\infty}^{\infty}c_{n}e\hspace{0.05cm}^{in\omega_0t}\hspace{4cm}・・・(3)\\\
c_n&=\frac{1}{T}\int_{-\frac{T}{2}}^{\frac{T}{2}}f(t)\hspace{0.05
cm}e\hspace{0.05cm}^{-in\omega_0t}dt\hspace{0.05cm}\hspace{2.39cm}・・・(4)\\
\end{align}
$$

　$F(n)$は複素数で、複素平面上の円周上の値となり、周波数番号$n$における振幅（|$F(n)$|）と位相（${\angle}F(n)$）を表します。



# 3. ちょっとおもしろい離散フーリエ変換の使い方

　ここで、離散フーリエ変換を用いて、$n$=$k$のときの振幅を求めることを考えます。周波数番号$k$のときの振幅は、|$F(k)$|で得られます。（2）式より（5）式が得られます。
$$
\begin{align}
F(k)&=\frac{1}{N}(f(0)e^{-i2{\pi}k\frac{0}{N}}+f(1)e^{-i2{\pi}k\frac{1}{N}}+・・・+f(N-1)e^{-i2{\pi}k\frac{N-1}{N}})\\\
&=\frac{1}{N}{\Bigl(}{\,}f(0)(cos(2{\pi}k\frac{0}{N})-i\hspace{0.05cm}sin(2{\pi}k\frac{0}{N}))+f(1)(cos(2{\pi}k\frac{1}{N})-i\hspace{0.05cm}sin(2{\pi}k\frac{1}{N}))+\\\
&\hspace{2cm}・・・+f(N-1)(cos(2{\pi}k\frac{N-1}{N})-i\hspace{0.05cm}sin(2{\pi}k\frac{N-
1}{N})){\,}{\Bigl)}{\quad}・・・(5)\\
\end{align}
$$

ここで、

$$
cos2{\pi}k\frac{j}{N}\hspace{15cm}
$$

と

$$
sin2{\pi}k\frac{j}{N}\hspace{15cm}
$$
の値が、-1, 0, 1となるような場合を考えてみます。例えば、$k$=1, $N$=4として、（5）式に代入しますと、（6）式が得られます。

$$
\begin{align}
F(1)&=\frac{1}{4}{\Bigl(}f(0)(cos(2{\pi}\frac{0}{4})-i\hspace{0.05cm}sin(2{\pi}\frac{0}{4}))+f(1)(cos(2{\pi}\frac{1}{4})-i\hspace{0.05cm}sin(2{\pi}\frac{1}{4}))+\\\
&\hspace{2cm}f(2)(cos(2{\pi}\frac{2}{4})-i\hspace{0.05cm}sin(2{\pi}\frac{2}{4}))+f(3)(cos(2{\pi}\frac{3}{4})-i\hspace{0.05cm}sin(2{\pi}\frac{3}{4})){\Bigr)}\\\
&=\frac{1}{4}{\Bigl(}f(0)(cos(0)-i\hspace{0.05cm}sin(0))+f(1)(cos(\frac{\pi}{2})-i\hspace{0.05cm}sin(\frac{\pi}{2}))+\\\
&\hspace{2cm}f(2)(cos(\pi)-i\hspace{0.05cm}sin(\pi))+f(3)(cos(\frac{3\pi}{2})-i\hspace{0.05cm}sin(\frac{3\pi}{2})){\Bigr)}\\\
&=\frac{1}{4}{\Bigl(}f(0)(cos(0)-i\hspace{0.05cm}sin(0))+f(1)(cos(\frac{\pi}{2})-i\hspace{0.05cm}sin(\frac{\pi}{2}))+\hspace{4.8cm}\\\
&\hspace{2cm}f(2)(cos(\pi)-i\hspace{0.05cm}sin(\pi))+f(3)(cos(\frac{3\pi}{2})-i\hspace{0.05cm}sin(\frac{3\pi}{2})){\Bigr)}\\\
&=\frac{1}{4}{\Bigl(}(f(0)-if(1)-f(2)+if(3){\Bigr)}\\\
&=\frac{1}{4}{\Bigl(}(f(0)-f(2))+i(f(1)-f(3)){\Bigr)}{\quad}・・・(6)\\
\end{align}
$$

　（6）式から、$k$=1, $N$=4のとき、（7）式で振幅と比例関係にある値$A$を求められます。

$$
A=\sqrt{(f(0)-f(2))^2+(f(1)-f(3))^2}{\quad}・・・(7)\\
$$

　位相を考えなくて良ければ（振幅だけ求めれば良いとき）、（8）式 が得られます。
$$
A=\sqrt{(f(j)-f(j+2))^2+(f(j+1)-f(j+3))^2}{\quad}・・・(8)\\
$$

　（8）式から、ターゲットの周波数に応じて、サンプリング周期を変えることで、比較的簡単な計算でターゲットの周波数の振幅を求めることがわかります。すなわちターゲットの周波数は、複素平面の単位円を一周することに相当しますので、単位円一周を４等分とした周期（ターゲット周波数の４倍の周波数）でサンプリングし（あるいはリサンプリングし）、（8）式を適用することで、ターゲットの周波数の振幅が得られます。下図は、その一例で、ターゲット周波数を0.5Hzとしたときの結果です。時間領域のグラフの黒丸がサンプリング点を示しています。時間領域で求めた振幅と（8）式から求めた振幅が比例関係にあることがわかります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/641763/13600e24-368b-123c-d651-d8e2340921d8.png)

　下図は、ターゲット周波数を2Hzとしたときの結果です。時間領域のグラフの黒丸がサンプリング点を示しています。ターゲット周波数に応じて、サンプリング周期を適切に選べば、時間領域で求めた振幅と（7）式から求めた振幅が比例関係にあることがわかります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/641763/e65da54a-2615-a77e-a6bb-32b733592114.png)

　下図は、ターゲット周波数を0.5Hzとしたときの結果ですが、先に示した0.5Hzの結果とは位相が異なっています。時間領域のグラフの黒丸がサンプリング点を示しています。ターゲット周波数に応じて、サンプリング周期を適切に選べば、時間領域で求めた振幅と（8）式から求めた振幅が比例関係にあることがわかります。

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/641763/0dfda69d-036b-9591-4ecd-236011868bda.png)

　以上、ターゲット周波数の４倍の周波数でサンプリングし（あるいはリサンプリングし）、（8）式を適用することで、ターゲットの周波数の振幅が得られることがわかります。コンピューティングパワーに制約があるときなど、演算負荷を減らしたいときなどに有用ではないでしょうか。


# 4. おわりに

　今回は、ちょっとおもしろいフーリエ変換の使い方をご紹介しました。フーリエ解析は、三角関数、複素数などが出てきて、数式も一見すると複雑に見えるかもしれませんが、数式の意味が理解できると、とても面白い解析法であることに気づくのではないでしょうか。今回の話もその一端ですが、理論を理解していくと使いこなしの面白さも感じるようになるかもしれません。<br>　今回も、理論の大切さをあらためて知りたいあるいは感じたいという方がいらっしゃいましたら、それをできるだけわかりやすく伝えられたら、という思いから記事を書かせて頂きました。フーリエ解析について、より詳しく知りたいという方は、参考文献などをご参考頂ければと思います。



# 参考文献

大石進一：フーリエ解析

金谷健一：これなら分かる応用数学教室

トランスナショナル カレッジ オブ レックス：フーリエの冒険