FFTを訪ねて［前編］電気屋としての理解を探る

はじめに

FFTと私

高速フーリエ変換(Fast Fourier Transform, FFT)は，通常の離散フーリエ変換(Discrete Fourier Transform, DFT)に比べて非常に高速なアルゴリズムとして知られており，ランダウの記号を用いると，DFTの計算量がO(N²)であるところ，FFTではO(N log N)になるようです。このため，理学・工学の幅広い分野で応用されております。言わば，ディジタル信号処理の中の1つの大きな“柱”となっています。ただし，入力される時間領域のディジタル信号のデータ点数が2のべき乗（例えば64とか1,024とか4,096など）となっている必要があります*1。

図1に，FFTのうち，本記事で取り扱う周波数間引き型Cooley-Tukeyアルゴリズムの計算途中経過を可視化したGIFアニメーションを示します。Python (JupyterLab)とMatplotlibを用いて作成しました*2。青が実部，橙が虚部を表しています*3。

図1: FFT（周波数間引き型Cooley-Tukeyアルゴリズム）の計算途中経過を可視化

ここでは，時間領域波形f(θ)として，

$f(\theta) = 0.8 e^{j \theta} + 0.6 e^{j(7 \theta - \pi / 2)} + 0.4 e^{j(12 \theta - \pi / 2)} + 0.3 e^{j(17 \theta + \pi / 6)} \tag{1}$

を与えています*4。jは虚数単位です(j² = −1)。また，θは図1の横軸全体で0から2πまで変化すると考えて下さい*5。FFTの結果，(1)式の通り，1次（基本波），7次，12次，17次にスペクトルが立っていることが分かります。データ点数はN = 64 (= 2⁶)です。

ところで，筆者はかなり以前から電力・電気系の（自称）技術者（俗にいう「電気屋」）として電力系統の高調波(harmonics)と向き合って参りました。当然ながらフーリエ解析 *6のお世話になるのですが，MATLABにしてもNumPyにしても，fftやnumpy.fft.fftとコマンドを打つだけでフーリエ解析ができてしまうので*7，その中身については理解が及んでおりませんでした。しかし，どうしてデータ点数が2のべき乗でなければならないのか，どうしてDFTよりも高速なのか，どのように時間波形からスペクトルを得ているか，本来は理解しないといけないのだろうけれども，よく分からないまま使っているという後ろめたさを抱え続けておりました*8。

また，以前に加速度センサを用いて地震計を作ろう思い立った際，加速度から震度を計算するアルゴリズムにFFTが必須であり*9，挫折した記憶があります。FFTをマイコン（例えばM5Stack Grayなど加速度センサのある系）に実装できれば，地震計の実験もできます。しかし，これも何年間もの休眠プロジェクトとなっています。

インターフェース誌の記事をきっかけとして

そんな中，CQ出版社の「インターフェース」誌の2021年9月号にて「数学とプログラミング」特集が組まれ，マイコン（記事ではソニーのSpresense）にFFTを実装する際の考え方とC言語によるプログラミングの例が記載されておりました*10。これはFFTへの理解を深めるための機運かもしれない。そこで，これまでに遊んできたPython (JupyterLab & Matplotlib)とM5Stack (Arduino IDE)にFFTを実装しながら，FFTの本質を探訪してみようと思い立ちました。

電気屋としての思考の道具

本記事は，筆者がFFTを理解しようとして~~試行錯誤~~七転八倒した結果をまとめたものです。PythonでのFFTの実装*11やMatplotlibでの計算途中経過のアニメーションを作成しながら考えを巡らせていくと，電気屋がもともと持っている思考の道具である，

α-β変換やd-q変換といった座標変換の概念
零相・正相・逆相といった対称座標法の概念
コモンモードとディファレンシャルモードといったモードの概念

との繋がりが深いことが分かってきました。

FFTの解説記事というのは世にあまたありますので，本記事の独自性として，このような電気屋としての思考の道具との関連性という切り口で書いていこうと思います。

M5StackとRaspberry Piへの実装

インターフェース誌2021年9月号に記載のプログラムを参考にして，M5Stack (Arduino IDE)とRaspberry Pi（PythonとC++)にDFTとFFTを実装し，計算時間を比較してみましたが，本記事だけでかなり長くなってしまったので，実装については次記事に譲ろうと思います。

参考資料

ディジタル信号処理の教科書の他，以下の3つを参考にさせて頂きました。なお，数式のノーテーションなどは，筆者が説明し易いように適宜書き換えております。

星野秋人：「マイコンと音声信号で体験！フーリエ変換プログラムの実装」，インターフェース，2021年9月号，p. 49-57
decafish：「2の基数のFFT」，腰も砕けよ膝も折れよ，2009年9月6日
J. W. Cooley and J. W. Tukey, “An algorithm for the machine calculation of complex Fourier series,” Mathematics of Computation, vol. 19, no. 90, Apr. 1965, pp. 297-301.

特に2のdecafishさんのブログは，以前にRaspberry Piのpigpioライブラリについて勉強した際にもお世話になっており*12，この場を借りて感謝を申し上げます。

3は本記事でご紹介しているCooley-Tukeyアルゴリズムを提案している原論文であり，2のべき乗に限らず，より一般的なフーリエ変換の高速化について述べています。この中では興味深いことに，純理論的には3のべき乗がもっとも効率的であるが，実際には（2進数で動いている）コンピュータにとって，2または4のべき乗にアドバンテージがあると書かれています。

The use of r_i = 3 is formally most efficient, but the gain is only about 6% over the use of 2 or 4, which have other advantages. If necessary, the use of r, up to 10 can increase the number of computations by no more than 50%. Accordingly, we can find “highly composite” values of N within a few percent of any given large number.

Whenever possible, the use of N = r^m with r = 2 or 4 offers important advantages for computers with binary arithmetic, both in addressing and in multiplication economy.

一般的な考え方

ここでは，まず連続時間でのフーリエ級数展開を出発点として，DFT，FFTの一般的な（数式ベースでの）考え方について説明します。

連続時間でのフーリエ級数展開とDFT

連続時間での信号をf(t)とします。f(t)が1周期をTとする周期性を持つとすれば，f(t)をフーリエ級数展開した場合のフーリエ係数F[k]は次式で書くことができます。

$F[k]= \displaystyle \frac{1}{T} \int_{0}^{T} f(t) e^{-j k \omega_{1} t} dt \tag{2}$

ここでkは整数であり*13，要するにF[k]とは「k次高調波成分の振幅と位相をフェーザとして表したもの」です。電気屋としてはこう呼んだ方がよほど分かりやすいですね。また，誤解する恐れがなければ，単に「k次高調波成分」と呼んでも差し支えないでしょう。また，F[0] ~ F[N]までの全体を指して「スペクトル」と呼んでも良いかと考えます。

ここではf(t)は複素数とします。もちろんf(t)は実数の信号としても構いませんが，その際にも「たまたま虚部がいつも零となっている複素数」と見なした方が便利です。この場合，後述しますが正相成分と逆相成分がお互いの複素共役となって打ち消しあっていると考えます。なお，ω₁は1周期がTとなる角周波数，すなわち基本波角周波数であり，

$\omega_{1} = \displaystyle \frac{2 \pi}{T} \tag{3}$

です*14。ここで，周期Tを“1”と見なして時間軸を規格化すれば，

$F[k] = \displaystyle \int_{0}^{1} f(t) e^{-j 2 \pi k t} dt \tag{4}$

と書くことができます。時間軸を規格化したので，本来はf(t)からf'(t')のように記号を変えるべきなのですが，面倒なのでf(t)のままとします。(4)式において，時刻tが0から1まで変化する時間をN分割したとすると，iを0からN − 1までの整数として，それぞれの時刻t[i]は

$t[i] = \displaystyle \frac{i}{N} \quad (i = 0, 1, \dots, N - 1) \tag{5}$

と書けます。そこで，f(t)をiの関数として時間軸上で離散化したものをf[i]と表記することにします*15。これに基づいて(4)式を離散化すると，

$F[k] = \displaystyle \sum_{i = 0}^{N - 1} f[i] e^{-j \frac{2 \pi k i}{N}} \tag{6}$

のように表すことができます*16。(6)式においてeの肩に乗っている部分が小さくて見えにくいので，新たにW_N^k[i]という記号を定義します*17。

$W_{N}^{k}[i] = e^{-j \frac{2 \pi k i}{N}} \tag{7}$

(7)式は回転因子とも呼ばれます。絶対値が1で，偏角が−2πki / Nの複素数であり，複素平面でf[i]を時計方向に回転させる効果があります*18。(7)式を使うと(6)式は

$F[k] = \displaystyle \sum_{i = 0}^{N - 1} f[i] W_{N}^{k}[i] \tag{8}$

と表記できます。(8)式が離散フーリエ変換(Discrete Fourier Transform, DFT)です。

さて，(8)式において，iもkも0からN − 1までの値をとります。このとき，F[0]からF[N − 1]までN個のフーリエ係数を計算しようとすると，個々のF[k]を計算するために，(8)式のシグマの中にあるf[i]とW_N^k[i]の乗算をN回実行する必要があります。すなわち，乗算の回数は単純にN²となります。したがって，Nが大きくなると加速度的に計算時間が増加していくことになってしまいます。

DFTからFFTへ

DFTよりも計算時間を節減するためのアルゴリズムとして考え出されたものの1つが高速フーリエ変換(Fast Fourier Transform, FFT)です。FFTの具体的なアルゴリズムの1つがCooley-Tukeyのアルゴリズムですが，本記事ではそのうちの周波数間引き(decimation-in-frequency)型に絞って述べていきます。ここでは，参考資料やWikipedia，その他のサイトと同じように，数式に基づいた説明を試みます。若干（いや，かなり？）天下り的な説明になってしまうと思います。

偶数次と奇数次に分ける

まず，(8)式において「k」を「2k」に置き換え，偶数次成分F[2k]のみを求める計算に書き換えると，

$F[2 k] = \displaystyle \sum_{i = 0}^{N - 1} f[i] W_{N}^{2 k}[i] \tag{9}$

となります。そう，周波数を「間引く」のです。(9)式のシグマでは，iを0からN − 1まで変化させて和を取っています。これに対して，f[i]の前半と後半のぞれぞれでf × Wを計算してその和をとり，それから半分の回数だけシグマを計算するように書き換えることもできます。このように書き換えても，すべてを足していることには変わりありませんよね。ここで，2のべき乗(= 2^m)であるNの半分をnとおきます。

$n = \displaystyle \frac{N}{2} \tag{10}$

このnを用いて(9)式を書き換えていくと，具体的には，

$\begin{align} F[2 k] &= \displaystyle \sum_{i = 0}^{N - 1} f[i] W_{N}^{2 k}[i] \\ &= \displaystyle \sum_{i = 0}^{n - 1} \left \{ f[i] W_{N}^{2 k}[i] + f[n + i] W_{N}^{2 k}[n + i] \right \} \\ &= \displaystyle \sum_{i = 0}^{n - 1} \left \{ f[i] W_{n}^{k}[i] + f[n + i] W_{n}^{k}[n + i] \right \} \\ &= \displaystyle \sum_{i = 0}^{n - 1} \left \{ f[i] W_{n}^{k}[i] + f[n + i] W_{n}^{k}[i] \right \} \\ &= \displaystyle \sum_{i = 0}^{n - 1} \left \{ f[i] + f[n + i] \right \} W_{n}^{k}[i] \end{align} \tag{10}$

のように変形できます。なお，(10)式中で回転因子Wが自由自在に変化しているように見えますが，次式からすべて等しいことが分かります。

$\begin{align} W_{N}^{2k}[i] &= e^{-j \frac{2 \pi 2 k i}{N}} = e^{-j \frac{2 \pi 2 k i}{2 n}} = e^{-j \frac{2 \pi k i}{n}} = W_{n}^{k}[i] \\ &= e^{-j (2 \pi k + \frac{2 \pi k i}{n})} = e^{-j \frac{2 \pi k (n + i)}{n}} = W_{n}^{k}[n + i] \end{align} \tag{11}$

最後から3番目の等号に関しては，kが整数であることに注意してください。(10)式の最初の行と最後の行を見比べると，fとWを乗算する回数が半分になっていることが分かると思います。言い換えると，f[i]のi = 0からi = NまでWと掛け算をしてから足し算（シグマ）をするのではなく，前半部f[i]と後半部f[n + i]を足し算をしてからWと掛け算をすることによって，掛け算の回数を節減していると言えます。

同様に，(8)式において「k」を「2k + 1」に置き換え，奇数次成分F[2k + 1]のみを求める計算に書き換えると，

$F[2 k + 1] = \displaystyle \sum_{i = 0}^{N - 1} f[i] W_{N}^{2 k + 1}[i] \tag{12}$

(10)式と同じように，n = N / 2を用いて，f[i]を前半と後半に分けるように式を変形していくと，

$\begin{align} F[2 k+ 1 ] &= \displaystyle \sum_{i = 0}^{N - 1} f[i] W_{N}^{2 k + 1}[i] \\ &= \displaystyle \sum_{i = 0}^{n - 1} \left \{ f[i] W_{N}^{2 k + 1}[i] + f[n + i] W_{N}^{2 k + 1}[n + i] \right \} \\ &= \displaystyle \sum_{i = 0}^{n - 1} \left \{ f[i] W_{N}^{2 k}[i] W_{N}^{1}[i] + f[n + i] W_{N}^{2 k}[n + i] W_{N}^{1}[n + i] \right \} \\ &= \displaystyle \sum_{i = 0}^{n - 1} \left \{ f[i] W_{n}^{k}[i] W_{N}^{1}[i] + f[n + i] W_{n}^{k}[n + i] W_{N}^{1}[n + i] \right \} \\ &= \displaystyle \sum_{i = 0}^{n - 1} \left \{ f[i] W_{n}^{k}[i] W_{N}^{1}[i] - f[n + i] W_{n}^{k}[i] W_{N}^{1}[i] \right \} \\ &= \displaystyle \sum_{i = 0}^{n - 1} \left \{ \left ( f[i] - f[n + i] \right ) W_{N}^{1}[i] \right \} W_{n}^{k}[i] \end{align} \tag{13}$

のように変形できます。なお，回転因子Wについて，次の2式が成り立ちます。

$\begin{align} W_{N}^{2 k + 1}[i] &= e^{-j \frac{2 \pi (2 k + 1) i}{N}} = e^{-j \frac{2 \pi 2 k i}{N}} e^{-j \frac{2 \pi i}{N}} \\ &= W_{N}^{2 k}[i] W_{N}^{1}[i] = W_{n}^{k}[i] W_{N}^{1}[i] \end{align} \tag{14}$

$\begin{align} W_{N}^{1}[n + i] &= e^{-j \frac{2 \pi (n + i)}{N}} = e^{-j \frac{2 \pi \left (\frac{N}{2} + i \right )}{N}} = e^{-j \left (\pi + \frac{2 \pi i}{N} \right )} \\ &= e^{-j \pi} e^{-j \frac{2 \pi i}{N}} = - e^{-j \frac{2 \pi i}{N}} = -W_{N}^{1}[i] \end{align} \tag{15}$

が成り立ちます。特に(15)式に着目すると，(13)式の下から2行目のように後半f[n + i]に掛かっている回転因子Wの符号が反転し，前半と後半の差が登場することが分かります。

再帰的な計算

(10)式と(13)式の結果のみを抜き出せば，以下のように書くことができます。

$\begin{align} F[2 k] &= \displaystyle \sum_{i = 0}^{n - 1} \left \{ f[i] + f[n + i] \right \} W_{n}^{k}[i] \\ F[2 k+ 1 ] &= \displaystyle \sum_{i = 0}^{n - 1} \left \{ \left ( f[i] - f[n + i] \right ) W_{N}^{1}[i] \right \} W_{n}^{k}[i] \end{align} \tag{16}$

ここで，次式のようにおいてみましょう（添え字の“e”はeven，“o”はoddの略です）。

$\begin{align} f_{e}[i] &= f[i] + f[n + i] \\[5pt] f_{o}[i] &= \left ( f[i] - f[n + i] \right ) W_{N}^{1}[i] \end{align} \tag{17}$

すると，(16)式は次のように書き換えられます。

$\begin{align} F[2 k] &= \displaystyle \sum_{i = 0}^{n - 1} f_{e}[i] W_{n}^{k}[i] \\[5pt] F[2 k+ 1 ] &= \displaystyle \sum_{i = 0}^{n - 1} f_{o}[i] W_{n}^{k}[i] \end{align} \tag{18}$

(18)式の右辺をよく見ると，それぞれがデータ点数をn = N / 2とするf_e[i]とf_o[i]のDFTとなっています。しかし，DFTは計算時間を要するので，FFTにするのが良いでしょう。すると，f[i]のFFTを計算するために，f_e[i]とf_o[i]のFFTを計算することになりますが，f_e[i]とf_o[i]をそれぞれ前半と後半に分け，言わば「2段階目のf_e[i]とf_o[i]」を計算することになります。この「2段階目のf_e[i]とf_o[i]」のデータ点数はnの半分，つまりN / 4となるはずです。このように続けていき，データ点数が2個となるところまで続けます。すなわち，再帰的な計算です。

なお，後述のように，(17)式の計算を「バタフライ演算(butterfly computation)」と呼びます。

参考資料に挙げたインターフェース誌の2021年9月号の記事は，再帰呼び出しを用いたC言語のプログラムを記載しています。一方，多重のforループを使えば同様の計算は可能であるため，必ずしも再帰呼び出しをする必要はありません。どちらの実装も可能です。しかし，再帰呼び出しを用いた実装の方がアルゴリズムを理解する助けになるというのが筆者の個人的な感想です（実際，筆者は再帰呼び出しのプログラムを見てやっと理解できた気になりました）。

シグナルフローグラフとブロック線図，ビット反転並べ替え

図2は，バタフライ演算を用いた（再帰的な）計算*19をシグナルフローグラフとして表現したものです。ここでは，データ点数がN = 8であった場合を例として描いています。N / 2 = 4点どうしの和と差を計算するところから始まって，最後は2点の演算になります。なお，本記事では取り扱いませんが，時間間引き(decimation-in-time)型のCooley-Tukeyアルゴリズムでは，バタフライ演算の形が（きわめて大雑把に言えば）左右反転します（2点の演算から始まって，N / 2点の演算で終わる）。

ただし，電気屋として生活してきた筆者としては，シグナルフローグラフよりもブロック線図の方が馴染み深いので，ブロック線図としても描いてみました。図3をご覧下さい*20。図3では，データを前半と後半に分けながら，より小さいDFTへと分割していく様子を，Stage 1, 2, 3として色分けしています。

図2のビット反転並べ替え部分を除いたものは，教科書やFFTの解説記事には必ずと言っていいほど登場しますが（図3は「邪道」なのかもしれません），恥ずかしながら筆者は，その意味するところをまったく理解できておりませんでした（シグナルフローグラフであるということすら認識できなかった💦）。図2，3は，(17)式のバタフライ演算で得たf_eとf_oを連結し，f_eを前半に，f_oを後半として元のfを上書きすることを表しています。つまり，図2，3のバタフライ演算による再帰的な計算を実行すると，元のf[i]は破壊されます*21。こう言うとデメリットのように聞こえますが，Cooley-Tukeyの原論文（参考文献3）を見ると，「データの記憶領域を節約できる」としてアドバンテージに挙げています。よく考えると，元のデータが，その記憶されているメモリ上の位置を変えずに，少しずつ変形しながらスペクトルに変わっていくというのは不思議ですよね。

また，本記事では図2，3の右側にビット反転並べ替えの部分を記載しました。前述の再帰的な計算では，最終結果として出てくるk次高調波成分（フーリエ係数）F[k]が順番通りに並んでいません。kを2進数で表した場合に，上位桁(MSB)と下位桁(LSB)が反転した形となります。これについては後述します。

計算量について

以上でFFTの周波数間引き型Cooley-Tukeyアルゴリズムについて，その一般的な考え方を筆者なりに説明してきました。このようなアルゴリズムでどうして高速化が図れるかと言えば，それは，乗算の回数を節減できるという点にあります。前述のように，データ点数がN点の場合，DFTでは乗算の回数がN²となります。一方，N = 8である場合を例示した図2，3を見ると，乗算の回数はN² = 64回ではなく12回です*22。一般に，データ点数N = 2^m（mは整数）である「2の基数のFFT」では，乗算の回数をN log₂ N / 2回に削減できます。

アルゴリズム	乗算の回数	ランダウの記号で表した計算量
DFT	$N^{2}$	$O(N^{2})$
FFT	$\displaystyle \frac{N \log_{2} N}{2}$	$O(N \log N)$

表1に，DFTとFFTの計算量をまとめます*23。どちらのアルゴリズムにも，乗算だけでなく加算やデータの移動があるので，トータルでの計算時間は乗算の回数だけでは厳密には決まりません。しかし，データ点数Nが大きい場合，漸近的にはランダウの記号で書かれるような計算量のオーダーとなります*24。Nが大きくなるとその差は歴然であり，次記事で述べますが，M5StackやRaspberry PiにDFTとFFTを実装して比較してみると，FFTの威力を実感できます。

電気屋としての疑問と考え方

さて，以上で数式をベースとした説明を試みてまいりましたが，やはり「イメージ」がないと筆者のような凡人にはうまく理解できません。そこで以下では，筆者がFFTについての一応の理解を得るまでに抱いた疑問点やそれに対する回答，また，電気屋しての思考の道具を応用した考え方について述べます。上述した数式ベースでの理解と互いに補完しあうことができれば，確固とした理解へと繋がるのではないかと期待します。

さまざまな疑問点

FFTについて筆者が持っていた疑問や，インターフェース誌の2021年9月号の記事を読み進めるうちに抱いた疑問は，主に以下の6点です。

入力される時間領域のディジタル信号のデータ点数はなぜ2のべき乗(2^m)でなければならないのか？
どのようにして計算を高速化しているのか？
入力されるディジタル信号を前半と後半に分けてその和と差を計算することにより，偶数次と奇数次に分けられるのはなぜか？
各次高調波成分はどのように計算されるのか？
バタフライ演算を行った場合，なぜ結果として得られるフーリエ係数の位置がビット反転してしまうのか？
入力されるディジタル信号が実数（虚部が零）である場合，FFT結果の高周波側に鏡像が現れるのはなぜか？

このうち，1と2についてはもう答えが出ています。1については，入力される時間領域波形のデータを「前半」と「後半」の「2つ」に分け続けるからであり*25，2については，（ごくごく定性的に言えば）足し算をしてから掛け算をするように順序を工夫したから，と言えます。

そこで以下では，疑問点3 ~ 6について考えていきます。なお，3については，上述の数式ベースの説明では一応の回答は得られていると思いますが，図を用いて補強できればと思います。

前半と後半の和と差を利用して偶数次と奇数次に分ける

(17)，(18)式で述べたように，周波数間引き型のCooley-Tukeyアルゴリズムでは，時間領域波形を前半と後半に分け，(i) 前半と後半の和(f_e)，(ii) 前半と後半の差に回転因子を乗算したもの(f_o)，という2つの信号を新たに作ります。これらf_eとf_oは元の時間領域波形とは独立にしておいても良いのですが，図2，3に示すように，f_eとf_oをそれぞれ新たな前半部，後半部として繋げ，元のデータfに上書きします。

例えば図4のように，1周期を“1”とする任意の時間領域信号f(t)，または1周期をN点としてf(t)を時間軸上で離散化したf[i]を考えます。このとき，前半と後半の和として得た新たな前半部には（直流成分）を含む偶数次高調波成分のみが，前半と後半の差として得た後半部には（基本波を含む）奇数次高調波成分のみが集まることが分かります。え？本当にそうなるの…？図5に，例として2次高調波と3次高調波の場合を示します。［追記: 2021-12-02］図5の2次高調波の部分が長いこと誤っていました。ようやく修正しました。

図5を見ると，2次高調波では，前半と後半でまったく同じ波形になっています。したがって，前半と後半の和として得た新たな前半における振幅は，元の波形の2倍になります。一方，前半と後半の差をとると，当然ながらキャンセルしてしまいますので，新たな後半（回転因子を乗算する前）はずっと零に貼り付いている波形となります。2次だけでなく，直流を含む他の偶数次高調波も同様です。

新たな前半のデータ点数は当たり前ですが元の半分，つまりn = N / 2となっています。このnを基準として考えると，元の2次高調波成分が1周期を描いています。言わば，元の2次高調波が新たな前半の中の基本波に変換されています。

一方，3次高調波では，前半と後半で逆位相となっています。したがって，前半と後半の和として得た新たな前半ではキャンセルされて零となり，一方，前半と後半の差として得た新たな後半（回転因子を乗算する前）では，元の波形の2倍となります。3次だけでなく，基本波を含む他の奇数次高調波も同様です。

新たな後半のデータ点数nを基準として考えると，元の3次高調波成分が1周期半を描いています。言わば，元の3次高調波が新たな後半の中の1.5次高調波に変換されています。ここで，偶数次とは異なり，新たな後半の始まりと終わりで波形が繋がらない（次数間高調波になってしまう）ことを覚えておいて下さい。

そこで，図6のように，これまで「新たな後半」と呼んでいた部分に，−0.5次，つまり逆相0.5次の回転因子を乗算します。これは(17)式中の回転因子W_N¹[i]に相当します。次式のように表現し直すと分かり易いかもしれません。

$W_{N}^{1}[i] = e^{-j \frac{2 \pi i}{N}} = e^{-j \frac{2 \pi i}{2 n}} = e^{-j \frac{\pi i}{n}} = W_{n}^{0.5}[i] \tag{19}$

そう，半周期だけ，時計方向に回転させる作用となります。

すると，1.5次高調波成分が何と（と言うか，回転因子を乗算しているので，複素数としては当たり前ですが），1次，つまり基本波成分に変換されます。3次高調波成分から半分になった1.5次高調波成分だけでなく，他の(2k + 1)次高調波成分が，時間軸を半分にして逆相0.5次の回転因子を乗算することにより，k次高調波成分に変換されます。

図7は，以上を踏まえて，2次高調波と3次高調波のみをそれぞれ含むディジタル信号に対してバタフライ演算を適用した場合のイメージです。

2次に代表される偶数次については，バタフライ演算によってデータ点数が半分となり，2k次高調波がk次高調波に変換されます。

そして，3次に代表される奇数次については，バタフライ演算によってやはりデータ点数が半分となり，逆相0.5次高調波の回転因子W_N¹[i] (= W_n^0.5[i])を演算することにより，(2k + 1)次高調波がk次高調波に変換されます。

図7では，新たな前半，後半として次のステージのバタフライ演算に進んでいきます。

各次高調波成分の計算

入力される時間領域のディジタル信号f[i]に含まれる各次高調波成分は，どのように求まるのでしょうか？ここでは，データ点数N = 16としてアルゴリズムを図解しながら考えていきます。

図8に周波数間引き型Cooley-Tukeyアルゴリズムを図解したものを示します。入力される時間領域のディジタル信号f[i]はf[0]からf[15]までの16点です。含まれる高調波成分は0次（直流），1次（基本波）から8次高調波までを仮定します*26。

最初のバタフライ演算で，前半には0，2，4，6，8次高調波成分，後半には1，3，5，7次高調波成分がそれぞれ抽出されます。前半については，半分のデータ点数（8点）を基準とすれば，0，1，2，3，4次高調波と見なすことができます（「変換される」と表現しても差し支えないと思います）。また，後半については，0.5，1.5，2.5，3.5次高調波となった後，回転因子W_N¹[i] (= W_n^0.5[i])と乗算することにより，半分のデータ点数（8点）を基準とした0，1，2，3次高調波に変換されます。

2段階目以降のバタフライ演算で，順次，元のデータの各次高調波成分が，その時点での偶数次と奇数次に別れていきます。データ点数が2点になった段階では，それら2つの和と差を計算します*27。これによって，最終的には1点ずつのデータに別れます。1点なので，データが存在すれば直流分（0次高調波成分）です。

最後に，得られた16点のデータにビット反転並べ替えを適用すれば，図8の最下段の結果が得られます。

さて，偶数次と奇数次の代表として，3次と4次の高調波を考えます。最下段の3次高調波成分に赤丸，4次高調波成分に青丸を付けてバタフライ演算を逆に辿って行くと，元の3次，4次高調波に至ることが分かります。他の次数についても同様に辿って行くことができます。

ビット反転されるメカニズム

図8にも関連しますが，バタフライ演算によって得られるフーリエ係数F[k]がkを2進数表記した上でビット反転した位置に格納されてしまうメカニズムについて考えてみます。図9は，ビット反転される様子を描いた概略図です。

k次高調波成分があるとします。バタフライ演算では，kが偶数であれば前半に含まれ，kが奇数であれば後半に含まれます。また，前半と後半のそれぞれでkは2で除算され，偶数次（前半）の場合はそのまま，奇数次（後半）の場合は小数点以下の「.5」の部分が回転因子によって削り取られます。言わば，奇数次についてはkが2で除算された後に，小数点以下が切り捨てられることになります。

さて，kを0から7までの数として，3桁の2進数で表現してみましょう。x，y，zをそれぞれ0か1であるとして，

$k = (xyz)_{b} \tag{20}$

とします。最初のバタフライ演算では，k = (xyz)_bが偶数であるか否かは，zが0であるか否かと等価です。図9のようにz = 0であれば前半として左へ，z = 1であれば後半として右へ行きます。kは2で割られ，仮にkが奇数だったとしても小数点以下は切り捨てられるので，これは右シフトに相当し， (xyz)_b次から (0xy)_b次となります。すると，今度はyが0か否かによって第2段階のバタフライ演算で前半となるか後半となるかが決まります。同様に，第3段階ではxが0か否かによって，前半となるか後半となるかが決まります。このようにして，0 = 000_b次から7 = 111_b次高調波成分が，LSBから順次MSBに向かって0か1か評価され，並べ替えられていきます。結果として，ビット反転された順序となります。

これは，FFTのデータ点数が2のべき乗であることが，2進数とも相性が良いことを利用した考え方です。データ点数が3のべき乗である場合にも，バタフライ演算のようなアルゴリズムを作ることができますが，その場合，FFT結果は「3進数でのビット反転」として並ぶことになるので，コンピュータでは取り扱い難いものとなるでしょう。

このように，2で除算して小数点以下を切り捨てた後の偶奇によってどんどん分岐していき，最終的に1ヶ所に1つの高調波成分が格納される形となります。最後に2つが1つに分かれる前，2次と6次とか，1次（基本波）と5次がペアとして残るので，何となく不思議な気もしますね。

FFT結果の高周波側に鏡像が含まれるメカニズム～正相と逆相～

本節では，これまでとは若干異なる議論となります。さて，FFTを応用する場面では，入力されるディジタル信号f[i]が実部のみを持っていることがほとんどだと思います。筆者も，MATLABやNumPy，Microsoft Excelでさまざまな信号（主に，電力系統やパワエレ機器の電圧・電流）をFFTしてきました。その中で経験的に，高周波側に鏡像が現れることを知りましたが，これまでその理由についてあまり深く考えたことがありませんでした（いずれにしても，k > N / 2の周波数成分についてはサンプリング定理から正しく測れないため）。例えば，図10を考えます。

図10: 入力されるディジタル信号f[i]が実部のみを持っていた場合のFFT結果

図10は冒頭に挙げた図1において，f[i]の虚部を零としたもの，すなわち(1)式の実部をとった

$\begin{align} f(\theta) = 0.8 \cos \theta &+ 0.6 \sin 7 \theta \\ &+ 0.4 \sin 12 \theta + 0.3 \cos (17 \theta + \pi / 6) \end{align} \tag{21}$

です。図10を見ると明らかなように，高周波側に鏡像が生じています。これは何故でしょうか？

結論から言うと，これまでに述べてきたFFTのアルゴリズムは「正相成分を正しくスペクトルに分解する手法」でした。入力されるディジタル信号f[i]が実部のみを持っている──とは，実は正相成分と逆相成分がともに存在し，互いに複素共役となっている状態と考えることができます。そう，高周波側の鏡像は，逆相成分を表しています。

図11に，正相成分と逆相成分が複素共役になっている場合の概略図を示します。正相成分と逆相成分の和をとると，虚部が打ち消し合い，実部だけが残ります（振幅は2倍になります）。CCW (counter clockwise)は反時計回り，CW (clockwise)は時計回りです。図10を再びよく見ると，高周波側の鏡像は，実部は低周波側と同じ，虚部は逆符号になっていることが分かります。つまり，複素共役になっています。

図12は図1と同じですが，入力されるディジタル信号f[i]が正相成分のみを持っている場合のFFT結果です。高周波の鏡像は消えており，各成分の振幅が2倍であることが分かります。

一方，図13はf[i]が逆相成分のみであった場合のFFT結果です。f[i]は次式で与えられます。

$f(\theta) = 0.8 e^{-j \theta} + 0.6 e^{-j(7 \theta - \pi / 2)} + 0.4 e^{-j(12 \theta - \pi / 2)} + 0.3 e^{-j(17 \theta + \pi / 6)} \tag{22}$

高周波の「鏡像」と表現していたものは，鏡像というよりは「逆相成分の実体」だったと言えるでしょう。

本記事でこれまでに説明してきたFFTのアルゴリズムでは，逆相成分が高周波側へと移動していきます。そして，逆相の基本波成分が最後の値F[N]となります。具体的にアニメーションで確認してみます。図14に逆相基本波成分を周波数間引き型Cooley-TukeyアルゴリズムにてFFTする場合の計算途中経過をアニメーション化したものを示します。

バタフライ演算によって，逆相基本波成分は奇数次高調波の1つとして後半に入りますが，−1次が−0.5次となり，さらに逆相0.5次（−0.5次）の回転因子を乗算されることにより，再び逆相基本波成分（−1次）となってしまいます。次の段階のバタフライ演算でも同じことが起こり，最後まで逆相基本波成分のまま，最後の要素F[63]となります。その他の成分もk = N −1 を逆相基本波成分としてkが小さくなる方向に並びます。余談ですが，PythonではリストやNumPyのarrayの添え字を負にすることによって，リストやarrayの終わりから逆方向に参照できます。これは，逆相成分を求めるのに都合が良さそうです。

ここで，これまで議論してきた実部と虚部を，α-β座標系でのα軸電圧・電流，β軸電圧・電流とそれぞれ考えれば，三相の電圧・電流の各次高調波における正相・逆相成分をFFTから抽出することができそうです。FFTをこのように捉えた考え方はあまり見たことがないので，文献を探してみようと思っています。

3のべき乗の場合

電気屋にとって3という数字はお馴染み

最後に，データ点数が2のべき乗ではなく，3のべき乗の場合のFFTについて簡単に触れたいと思います。本記事の中でも何度か言及していますが，データ点数Nが2のべき乗でなくともCooley-Tukeyアルゴリズムのようなことが可能です。

まず，N = 2^mの場合の(17)，(18)式を(23)，(24)式として再掲します。

$\begin{align} f_{e}[i] &= f[i] + f[n + i] \\ f_{o}[i] &= \left ( f[i] - f[n + i] \right ) W_{N}^{1}[i] \end{align} \tag{23}$

$\begin{align} F[2 k] &= \displaystyle \sum_{i = 0}^{n - 1} f_{e}[i] W_{n}^{k}[i] \\[5pt] F[2 k+ 1 ] &= \displaystyle \sum_{i = 0}^{n - 1} f_{o}[i] W_{n}^{k}[i] \end{align} \tag{24}$

ここから「類推」して，N = 3^mの場合の式を書くことができます。ここではnを

$n = \displaystyle \frac{N}{3} \tag{25}$

とします。入力されるディジタル信号f[i]を前半，後半の2つではなく，前部，中部，後部の3つに分けます*28。

$\begin{align} f_{0}[i] &= f[i] + f[n + i] + f[2n + i] \\[5pt] f_{1}[i] &= \left ( f[i] + a^{2} f[n + i] + a f[2n + i] \right ) W_{N}^{1}[i] \\[5pt] f_{2}[i] &= \left ( f[i] + a f[n + i] + a^{2} f[2n + i] \right ) W_{N}^{2}[i] \end{align} \tag{26}$

$\begin{align} F[3 k] &= \displaystyle \sum_{i = 0}^{n - 1} f_{0}[i] W_{n}^{k}[i] \\ F[3 k+ 1 ] &= \displaystyle \sum_{i = 0}^{n - 1} f_{1}[i] W_{n}^{k}[i] \\ F[3 k+ 2 ] &= \displaystyle \sum_{i = 0}^{n - 1} f_{2}[i] W_{n}^{k}[i] \end{align} \tag{27}$

(26)式と(27)式はまさにバタフライ演算の亜種であり，N = 3^mの場合の周波数間引き型Cooley-Tukeyアルゴリズムと言えます。さて，改めて(26)式を見てみると，何やら謎のaという係数が乗算されています。そして，(26)式自体，電気屋としてはとてもとても見覚えがあるものです。そう，お察しの通り，

$a = \displaystyle e^{j \frac{2 \pi}{3}} = -\frac{1}{2} + j \frac{\sqrt{3}}{2} \tag{28}$

です。これはまさに，対称座標法そのものの式となっています。つまり，入力されるディジタル信号f[i]の前部，中部，後部を三相と見なした場合の零相，正相，逆相成分を抽出することで，再帰的な計算によるFFTを可能としています。

計算例

数式だけではつまらないので，実際にやってみましょう。データ点数N = 81 (= 3⁴)として，(1)式と同様に正相成分*29のFFTを試みます。

図15に結果を示します。1次（基本波），7次，12次，17次にスペクトルが現れており，正しくFFTできていることが分かります。

まとめ

以上，本記事では，FFTとして周波数間引き型Cooley-Tukeyアルゴリズムについて，筆者として悪戦苦闘した過程と結果を（若干アットランダムに）計算例を示しながら説明しました。まず，一般的な数式を用いたFFTの導出からバタフライ演算のシグナルフローグラフとブロック線図について述べ，その後，各種の疑問にできるだけ回答しました。特に，電気屋として持っている思考の道具である，座標変換などの概念との対比を試みました。

ところで，残っている疑問点もまだまだあります。例えば，表1の通り，DFTの計算量がO(N²)，FFTの計算量がO(N log N)であるとされていますが，これは，足し算よりも掛け算の計算量が十分に大きいという前提に立っていると思います。しかし，もしハードウェア浮動小数点乗算器を十分に持っている計算機や，マルチコアによる並列演算ができる場合，また，データの移動や読み書きに制約があるような系の場合，掛け算のみに焦点を当てるのは不適当で，場合によってはO(N²)やO(N log N)という枠組みが当てはまらなくなることもあるのではないかと言う気がしています*30。

本来はこの記事の中で，Python (JuputerLab)やM5StackへのDFTとFFTの実装や，計算時間の比較についても書くつもりでおりましたが，この時点でかなりの分量となってしまったため，マイコン類への実装については次記事に譲りたいと思います。

筆者としては，これまで，FFTをブラックボックスとして取り扱い続けてしまったことを恥じております。本来であれば，今から10 ~ 15年前にはこの記事のような理解を得ているべきでしたが，怠け続けていた年月があまりにも長いため，このような体たらくとなっております。しかし…！そう，今からでも学びの一歩一歩を再び踏みしめる所存です。

*1:後述しますが，実は3のべき乗でもそれ以上の自然数のべき乗でも，アルゴリズムは存在します。ただし若干効率は下がると思います。

*2:図1では，FFTの結果を見やすくするため，それぞれのバタフライ演算の後でデータを2で除算しています。

*3:最近，1つ目の要素を青，2つ目の要素を橙にするのが流行っていませんか？ MatplotlibもMicrosoft Excelもそうですよね。

*4:このデータは正相成分のみをもっています。

*5:ちょっと乱暴な書き方ですが…ご容赦下さい…。

*6:フーリエ級数展開なのかフーリエ変換なのか，はたまた連続なのか離散なのか，ぼやかして言うときに「フーリエ解析」という言い方は便利だと思っています。

*7:Microsoft Excelでも「分析ツール」でFFTが可能です。

*8:本ブログの筆者の数学力は驚くほど低いことが知られています。

*9:気象庁｜強震観測について｜計測震度の算出方法

*10:2021年9月号｜ Interface – CQ出版

*11:もちろん，前述のnumpy.fft.fftは使わず，自分でFFTの計算を書きます。

*12:Raspberry PiでBLDCモータを120°通電制御する - The Negligible Lab

*13:kが負になった場合は，逆相成分を抽出することになります。

*14:標準的な教科書などではω₀やΩ₀などと添え字を“0”にしている場合が多い気がしますが，「基本波」であること鑑み，本記事では添え字を“1”としました。

*15:なぜか，離散化するとf[i]のように大括弧（角括弧）で書くか，それともf_iのように下付きにする慣習がありますね。

*16:Nで割らないのが一般的な書き方のようです。したがって，(6)式のF[k]は(4)式のF[k]よりもN倍大きくなります。

*17:FFTに関する教科書や解説記事では「[i]」を付けていない場合がほとんどだと思いますが，iの関数だと分からなくなってしまわないように，本記事では「[i]」を付けることにします。

*18:この時点で既にd-q変換との関連性が疑われてきますね。

*19:図2，3では元のf[i]を上書きしていくようなアルゴリズムを描いています。計算機としてメモリの使用量を節減するためにこのようなアルゴリズムとなっていると筆者は理解しています。その副作用として，計算結果であるフーリエ係数F[k]が格納される位置がビット反転した値となってしまいます。

*20:所詮は有限のデータ点数である！気合を入れればLTspiceでも（制御ブロックとして）実装できるのではないか⁉

*21:もちろんコピーを作ってから図2，3のアルゴリズムを実行すれば元のデータは保存されます

*22:当たり前ですが，図3において回転因子Wの入ったブロックの数です

*23:最近，クイックソートの計算量がO(N log N)であることは常識か否か云々，といった議論がTwitterで勃発しており，不用意にランダウの記号を使うのはシロートといてちょっと怖いです😅

*24:これについても疑問があり，後述します。

*25:前部，中部，後部の3つに分けることによって，何とデータ点数が3のべき乗の場合のFFTも作ることができます。

*26:9次以降はサンプリング定理からエイリアスを生じてしまうため。

*27:後半に対する回転因子W_N¹[i]がi = 0の場合に1となるため。

*28:序盤，中盤，終盤あるいは初旬，中旬，下旬と言っても良い？

*29:ここで言う正相成分は，(26)式のf₁とは概念として異なりますね…。より適切な表現があれば良いのですが…💦

*30:繰り返しますが計算量についてはドシロートです😅

The Negligible Lab

Astray in the forest of electrical and electronic circuits. Adrift in the gap between time and frequency domains. 独立独歩を是とする。

FFTを訪ねて［前編］電気屋としての理解を探る

はじめに

FFTと私

インターフェース誌の記事をきっかけとして

電気屋としての思考の道具

M5StackとRaspberry Piへの実装

目次

参考資料

一般的な考え方

連続時間でのフーリエ級数展開とDFT

DFTからFFTへ

偶数次と奇数次に分ける

再帰的な計算

シグナルフローグラフとブロック線図，ビット反転並べ替え

計算量について

電気屋としての疑問と考え方

さまざまな疑問点

前半と後半の和と差を利用して偶数次と奇数次に分ける

各次高調波成分の計算

ビット反転されるメカニズム

FFT結果の高周波側に鏡像が含まれるメカニズム～正相と逆相～

3のべき乗の場合

電気屋にとって3という数字はお馴染み

計算例

まとめ

はじめに

FFTと私

インターフェース誌の記事をきっかけとして

電気屋としての思考の道具

M5StackとRaspberry Piへの実装

目次

参考資料

一般的な考え方

連続時間でのフーリエ級数展開とDFT

DFTからFFTへ

偶数次と奇数次に分ける

再帰的な計算

シグナルフローグラフとブロック線図，ビット反転並べ替え

計算量について

電気屋としての疑問と考え方

さまざまな疑問点

前半と後半の和と差を利用して偶数次と奇数次に分ける

各次高調波成分の計算

ビット反転されるメカニズム

FFT結果の高周波側に鏡像が含まれるメカニズム ～正相と逆相～

3のべき乗の場合

電気屋にとって3という数字はお馴染み

計算例

まとめ

FFT結果の高周波側に鏡像が含まれるメカニズム～正相と逆相～