FFTを訪ねて［後編］PythonとC++で作ってみる
はじめに
11月に入っていよいよ秋が深まって参りましたが，いかがお過ごしでしょうか？ さて，前記事にて，高速フーリエ変換(Fast Fourier Transform, FFT)に対する筆者なりの理解の仕方をまとめました。周波数間引き(decimation-in-frequency)型のCooley-Tukeyアルゴリズムや，その中で要となるバタフライ演算について述べた他，α-β変換やd-q変換など，電気屋としての思考の道具との関連性について考えてみました。
その続きとなる本記事では，FFTを実現する関数*1をPythonとC++で実装した様子や，それぞれのアルゴリズムを
- Raspberry Pi (3 Model A+ & Zero WH)
- Windowsパソコン (Anaconda3)
- Androidスマートフォン (PyDroid 3 on Xperia XZ1)
- M5Stack (Arduino IDE)
- NumWorksの関数電卓 (MicroPython)
にて動かしてみた結果について書いていこうと思います。また，FFTの威力を実感するため，普通の離散フーリエ変換(Discrete Fourier Transform, DFT)も同様に実装して比較することとしました。
なお，PythonとC++を選んだ経緯ですが，まず，(i) JupyterLabでアルゴリズムを対話的に確認したかったためPythonで作った，(ii) 単体でも動くようにJupyter Notebook (.ipynb)ファイルではなく単体の.pyファイルとした，(iii) M5StackやESP32 DevKit Cなどのため，Arduino IDEでも使えるようにC++に移植した（Mbedでも使えるかも）。また，複素数ライブラリについて若干の改造が必要ながら，C++はRaspberry Piでも使える，と言ったところです。
本記事でも前記事と同様に，下記のインターフェース誌の記事を参考にさせて頂きました。とは言え，PythonやC++に合わせて最初から作り直しています。
特に，上記記事のp. 53に記載のリスト1，2を参考にさせて頂いており，本記事では略して「イ誌リスト1」，「イ誌リスト2」として言及することにします。
動作風景
まずは実際に動いているところをご覧下さい。M5Stackのアナログ入力端子(GPIO35)に三角波と矩形波の信号を入力し，analogRead関数でサンプリングしています。得られた信号の時間領域での波形と，それをFFTしたスペクトルを液晶ディスプレイに表示するプログラムを作ってみました。動画1をご覧ください*2。
それぞれ奇数次高調波と思われるスペクトルを確認できますね。
図1に示すように，Nucleo L476RGを信号源として使いました。本記事ではNucleoにはFFTを実装しませんでしたが，Mbedに移植してみても良いかと思います。
目次
Pythonでの実装
Pythonではもちろん，NumPyのnumpy.fft.fft関数で簡単かつ高速にFFTできます*3。しかしそれではもちろん面白くないので，自分で周波数間引き型Cooley-Tukeyアルゴリズムを実装することにします。
Pythonでの数値計算や信号処理というと，やはりNumPyを使うことになると思いますが，本記事では，(i) NumPyの配列演算を活用したバージョン（NumPy版）と，(ii) NumPyを使わず，cmathモジュールだけでFFTを実現するバージョン（cmath版）をそれぞれ作ってみることにしました。特に後者はNumPyを使えないMicroPythonのような環境に役立つのではないかと思っています。
なお，FFTについては表1に示す3タイプを作ることにしました（cmath版はType IIIのみ）。
| タイプ | アルゴリズム | 備考 |
|---|---|---|
| FFT Type I | 再帰呼び出しを用い，フーリエ係数の並べ替えにはテンポラリな配列*4を用いる方式 | イ誌リスト1に相当 |
| FFT Type II | 再帰呼び出しを用い，フーリエ係数の並べ替えビット反転並べ替えを用いる方式 | イ誌リスト2に相当 |
| FFT Type III | 再帰呼び出しではなく多重forループを用い，フーリエ係数の並べ替えビット反転並べ替えを用いる方式 |
Pythonでの関数呼び出しにおける引数は，ミュータブルな型に対してはいわゆる「参照渡し」です。NumPyのndarrayやPythonのリストはミュータブルであるため，関数に渡した場合，関数の中でそのndarrayやリストを変更すると，呼び出し元にもその変更が影響します。本記事ではこれを利用して，与えられたndarrayやリストを書き換えていくような関数として実装しました。Cooley-Tukeyアルゴリズムでは元の信号を書き換えながら演算するため，この方法との相性が良いと考えます。
なお，本記事ではアルゴリズムを学ぶために必要な最小限の機能しか実装していません。例えば，与えられた入力信号のデータ点数Nが本当に2のべき乗であるかチェックしたり，Pythonの場合は引数の型がそもそもNumPyのndarrayやリストではなかった場合の例外処理など，実用の際に必要となる機能はありませんのでご注意ください。
NumPy版
時間領域の入力信号f[i]がNumPyのndarrayに格納されているとします。まず，FFTとの比較用としてDFTを作ります。フーリエ係数F[k]をF[0]からF[N]までひとつひとつ計算していくプログラムです。次に，イ誌リスト1，2と同様に，再帰呼び出しを用いてバタフライ演算を繰り返し，元のf[i]を書き換えながらフーリエ係数F[k]を得るプログラムを書いていきます。
Pythonでは複素数を（NumPyをインポートしなくても）ナチュラルに扱うことができます。また，NumPyでは配列と配列の和や積を計算する際に，1つ1つの要素に対するforループを作る必要がありません。その2点で，イ誌リスト1，2（実部と虚部を別の変数としています）とは異なる作りになっています。
なお，以下で示すそれぞれの関数を呼ぶ前に，定番のおまじないであるimport numpy as npが実行されていると仮定します。
DFT
まずはDFTです。前記事のおさらいになりますが，DFTは次式のように入力信号f[i]からフーリエ係数F[k]を得る操作です。F[k]をスペクトルとも呼ぶことにします。
ただし，WNk[i]は回転因子と呼ばれ，
と表されます*5。これは単に(2)式右辺のeの肩の上の文字が小さくて見えづらいので，WNk[i]と略記しているものです。
これをPythonのプログラムに愚直に置き換えます。直流成分であるF[0]から初めて，F[1]，F[2]，…，F[N]までfor文で真面目にひとつひとつ積和演算を繰り返してしていきます。ただし，NumPyを活用しているのでforループは1重で済んでしまいますね。内積を計算（積和演算）する@演算子の威力です。
# DFT, with NumPy def DFT(f): N = len(f) i = np.arange(N) F = np.zeros_like(f) # To contain calculation results # DFT core calculation for k in range(N): # Prepare complex sinusoid, rotation operator W = np.exp(-2j * np.pi * k * i / N) # multiply-accumulate computation F[k] = W @ f # Return result return F
DFTでは入力信号f[i]を破壊することなく計算結果F[k]が得られるので，return FとしてスペクトルFを返すことができます。ただし，以降で述べるCooley-TukeyアルゴリズムによるFFTと関数の使い方を揃えるため，return Fの代わりに
# Overwrite results to input
f[:] = F
として，敢えて入力信号f[i]を上書きする実装も考えられ，本記事ではそちらを使うことにしました。なお，f = Fではなくf[:] = Fとした理由については付録にて後述します。
FFT Type I
いよいよFFTをPython + NumPyで実装します。まずはイ誌リスト1と同様のアルゴリズムであるType Iです。
NumPyの配列演算の力を活用してバタフライ演算とスペクトルの並べ替えを実装すれば，何とforループがなくなってしまいます。スライスを活用すれば，f[:n]が前半を，f[n:]が後半を表すという，シンプルな書き方ができます。
# FFT Type I, with NumPy # Caution! Original input data is overwritten and destroyed. def FFT_1: N = len(f) if N > 1: n = N // 2 # Do not use '/', which gives float. i = np.arange(n) # Prepare complex sinusoid, rotation operator W = np.exp(-2j * np.pi * i / N) # Butterfly computation f_tmp = f[:n] - f[n:] f[:n] += f[n:] f[n:] = W * f_tmp # Recursively call this function itself FFT_1(f[:n]) # First half FFT_1(f[n:]) # Second half # Simple permutation F = np.empty_like(f) F[0::2] = f[:n] F[1::2] = f[n:] # Overwrite results to input f[:] = F
また，スペクトルの並べ替え部分では，イ誌リスト1と同様にテンポラリな配列Fを用意しています*6。F[0::2]が偶数次成分F[0]，F[2]，F[4]，…，F[N]を，F[1::2]が奇数次成分F[1]，F[3]，F[5]，…，F[N − 1]を表しますので，イ誌リスト1のC言語プログラムと比較してシンプルに書けます。
FFT Type II
こちらはイ誌リスト2と同様に，バタフライ演算をすべて終えてからビット反転並べ替えを行う方法です。バタフライ演算を実行する関数_FFT_2_butterflyと，メイン関数から呼び出される本体である関数FFT_2の2つに分かれています*7。
# FFT Type II, with NumPy def _FFT_2_butterfly(f): N = len(f) if N > 1: n = N // 2 # Do not use '/', which gives float. i = np.arange(n) # Prepare complex sinusoid W = np.exp(-2j * np.pi * i / N) # Butterfly computation f_tmp = f[:n] - f[n:] f[:n] += f[n:] f[n:] = W * f_tmp # Recursively call this function itself _FFT_2_butterfly(f[:n]) # First half _FFT_2_butterfly(f[n:]) # Second half def FFT_2(f): # Butterfly computations _FFT_2_butterfly(f) # Bit-reversal permutation i = 0 for j in range(1, N - 1): k = N >> 1 while(True): i ^= k if k > i: k >>= 1 continue else: break if j < i: f_tmp = f[j] f[j] = f[i] f[i] = f_tmp
ビット反転については，どうしてこの書き方で反転できるのか，筆者はまだよく理解していません…。何か分かった時点で付録に追記します。
FFT Type III
イ誌リスト1，2とは異なり，多重のforループを用いることにより関数の再帰呼び出しを行わない方法でも実装してみました。スタックの節減になると考えられる他，後述のcmath版では，このタイプのみを実装しています*8。
周波数間引き型Cooley-Tukeyアルゴリズムにおけるバタフライ演算のブロック線図を前記事から再掲します。
図2において，ピンク，水色，緑でハッチングしたそれぞれのバタフライ演算を，forループで愚直に廻していけば，関数の再帰呼び出しをしなくてもFFTが可能です*9。ステージ1（ピンク）ではデータ全体を前半と後半に分けてバタフライ演算を実行します。ステージ2ではデータ全体を2つに分け，それぞれの中でバタフライ演算を実行し，さらにステージ3ではデータ全体を4つに分け…というように，多重のforループを回していきます。
# FFT Type III, with NumPy def FFT_3(f): N_orig = len(f) m = int(np.log2(N_orig)) # Butterfly computations for s in range(0, m): # Stage counter div = 2**s # Number of divided data segments N = N_orig // div # Number of data points in each segment n = N // 2 # Do not use '/', which gives float # Prepare complex sinusoid, rotation operator i = np.arange(n) W = np.exp(-2j * np.pi * i / N) # Butterfly computations of each segment for j in range(div): # Counter for divided data segments # Butterfly computation f_tmp = f[N * j : N * j + n] \ - f[N * j + n : N * (j + 1)] f[N * j : N * j + n] \ += f[N * j + n : N * (j + 1)] f[N * j + n : N * (j + 1)] = W * f_tmp # Bit-reversal permutation i = 0 for j in range(1, N_orig - 1): k = N_orig >> 1 while(True): i ^= k if k > i: k >>= 1 continue else: break if j < i: f_tmp = f[j] f[j] = f[i] f[i] = f_tmp
多重forループのイメージとしては表2となります。なお，Norig (N_orig)は入力データf[i]のデータ点数です。
| 階層 | ループ変数 | ループ回数 | ループ変数が表す内容 | 処理内容 |
|---|---|---|---|---|
| 1 | s | log2 Norig | バタフライ演算のステージ | 図1の各ステージでデータ全体を何分割するか計算 |
| 2 | j | 2s | 各ステージで何番目のバタフライ演算か | |
| 3 | i | Norig / 2s + 1 | バタフライ演算を実行中の要素番号 | 図1の各ステージで分割されたそれぞれのデータに対してバタフライ演算を実行 |
ただし，NumPyの配列演算を活用している場合，バタフライ演算を具体的に実行している3階層目にはforループが不要となります。後述のcmath版ではforループとして実装します。
バタフライ演算部で入力信号fに与えられているスライスが複雑に見えますが，これは図1における例えば緑色で囲まれた部分1つの前半と後半を抽出しているだけです。
cmath版
様々な理由でNumPyを使えない環境もあるかと思います。NumWorksの関数電卓のMicroPythonもその一例ですね。筆者として試したことはありませんが，M5Stackなどで使えるMicroPythonやCircuitPythonも同様であると推測します。そこで，cmathモジュールのみでFFTするような実装も考えてみました。
入力信号f[i]はNumPyのndarrayではなくPythonのリストとして格納しておきます。
DFT
まずはDFTです。NumPyとは異なりcmathには配列演算の機能がないので，forループを1つ増やして，ひとつひとつ計算します。
# DFT, with cmath def DFT_cmath(f): N = len(f) F = [0] * N # Define list to contain calculation results, 0-initialized # DFT core calculation for k in range(N): for i in range(N): # Prepare complex sinusoid, rotation operator W = cmath.exp(-2j * cmath.pi * k * i / N) # multiply-accumulate computation F[k] += W * f[i] # Return result return F
NumPy版と同様に，return Fの代わりに
# Overwrite results to input
f[:] = F
とすることもできますね。
FFT Type III
付録で述べますが，NumPyのndarrayとPythonのリストでは，スライスに対する操作の挙動が異なります。cmath版では入力信号f[i]をリストとして与えますが，この際にイ誌リスト1，2のような関数の再帰呼び出しの方法では意図通りに計算できなかったため*10，再帰呼び出しを用いないType IIIのみを実装しました。NumPyの配列演算を活用できないため，forループが1つ増え，表2の通りに3重になっています。
ただし，ビット反転並べ替えの箇所は，一言一句，NumPy版と同じになっています。
# FFT Type III with cmath def FFT_3_cmath(f): N_orig = len(f) m = int(cmath.log(N_orig, 2).real) # FFT core calculation for s in range(0, m): # Stage counter div = 2**s # Number of divided data segments N = N_orig // div # Number of data points in each segment n = N // 2 # Do not use '/', which gives float. # Butterfly computations of each segment for j in range(div): # Counter for divided data segments for i in range(n): # Counter for each element in data segment # Prepare complex sinusoid, rotation operator W = cmath.exp(-2j * np.pi * i / N) # Butterfly computation f_tmp = f[N * j + i] - f[N * j + n + i] f[N * j + i] += f[N * j + n + i] f[N * j + n + i] = W * f_tmp # Bit-reversal permutation i = 0 for j in range(1, N_orig - 1): k = N_orig >> 1 while(True): i ^= k if k > i: k >>= 1 continue else: break if j < i: f_tmp = f[j] f[j] = f[i] f[i] = f_tmp
C++での実装
筆者はもともと，Arduino IDEにてM5Stack向けのFFTの実装を作ることを考えていたので，C++のプログラムを考えていました。ただし，Linux環境であるRaspberry Piでも当然C++プログラムをコンパイル，実行できます。そこで，より一般的と思われるRaspberry Pi向けの実装をまずご説明します。
C++では時間領域の入力信号を浮動小数点複素数型（例: complex<double>型）の配列に格納し，DFTやFFTする関数にはその先頭要素へのポインタと要素数を渡します。Pythonとは異なり，CやC++では，先頭要素へのポインタとして配列を関数に渡すのですが，関数の側では配列の要素数の情報が失われてしまうため，PythonにおけるN = len(f)のような便利なことができません。したがって，要素数を一緒に渡すのが定石のようです*11。
Raspberry Pi版
DFT
まずはやはりDFTですね。C++では#include <complex>を置けば，complex<double>型として浮動小数点複素数を取り扱うことができるようになりますので，イ誌リスト1，2のように実部と虚部を分けて記述する必要はなくなります。なお，#include <cmath>も必要です。また，using namespace std;を書いている前提となっていますのでご注意ください。
NumPyのような配列演算はないので，forループで要素をひとつひとつ取り出して計算するPythonのcmath版に近い形になります。
// DFT, in C++ void DFT(complex<double> *f, int N) { // Declare complex array to store output complex<double> F[N]; // Straight-forward calculation of Fourier coefficients for (int k = 0; k < N; k++) { for (int i = 0; i < N; i++) { // Prepare complex sinusoid, rotation operator complex<double> W = polar(1.0, -2 * M_PI * k * i / N); // Perform multiply-accumulate operation F[k] += W * f[i]; } } // Overwrite input with output for (int i = 0; i < N; i++) { f[i] = F[i]; } }
なお，4行目でFをcomplex<double>型の配列として生成していますが，complexクラスのdouble型に特殊化されたコンストラクタにて自動的に0 + j0に初期化されるようです*12。
FFT Type I
FFT Type Iでは，複素数としてcomplex<double>型を利用している点が異なりますが，イ誌リスト1に近い実装となります。
// FFT Type I, in C++ void FFT_1(complex<double> *f, int N) { if (N > 1) { int n = N / 2; for (int i = 0; i < n; i++) { // Prepare complex sinusoid, rotation operator complex<double> W = polar(1.0, -2 * M_PI * i / N); // Butterfly computation complex<double> f_tmp = f[i] - f[n + i]; f[i] += f[n + i]; f[n + i] = W * f_tmp; } // Recursively call this function itself FFT_1(&f[0], n); // First half FFT_1(&f[n], n); // Second half // Simple permutaton complex<double> F[N]; for (int i = 0; i < n; i++) { F[2 * i] = f[i]; F[2 * i + 1] = f[n + i]; } for (int i = 0; i < N; i++) { f[i] = F[i]; } } }
FFT Type II
FFT Type IIも同様です。バタフライ演算を行う関数_FFT_2_butterflyは他のファイルに含まれる関数から呼ばれることがないため，static属性を付しています。
// FFT Type II, in C++ static void _FFT_2_butterfly(complex<double> *f, int N) { if (N > 1) { int n = N / 2; for (int i = 0; i < n; i++) { // Prepare complex sinusoid, rotation operator complex<double> W = polar(1.0, -2 * M_PI * i / N); // Butterfly computation complex<double> f_tmp = f[i] - f[n + i]; f[i] += f[n + i]; f[n + i] = W * f_tmp; } // Recursively call this function itself _FFT_2_butterfly(&f[0], n); // First half _FFT_2_butterfly(&f[n], n); // Second half } } void FFT_2(complex<double> *f, int N) { // Perform butterfly computations _FFT_2_butterfly(f, N); // Bit-reversal permutation int i = 0; for (int j = 1; j < N - 1; j++) { for (int k = N >> 1; k > (i ^= k); k >>= 1); if (j < i) { complex<double> f_tmp = f[j]; f[j] = f[i]; f[i] = f_tmp; } } }
FFT Type III
FFT Type IIIも同様です。実装の根本的な部分ではPythonのcmath版のFFT Type IIIとほとんど共通しています。4行目で要素数の2を底とする対数を計算するところが若干違いますね。また，累乗を計算する際に，Pythonでは**演算子で計算できましたが，C++ではpow関数が必要です。しかし違いはそれぐらいのものです。
// FFT Type III, in C++ void FFT_3(complex<double> *f, int N_orig) { // Calculate number of data divisions int m = int(log2(N_orig)); // FFT core calculation for (int s = 0; s < m; s++) { int div = int(pow(2, s)); // Number of divided data segments int N = N_orig / div; // Number of data points in each segment int n = N / 2; // Butterfly computations of each segment for (int j = 0; j < div; j++) { for (int i = 0; i < n; i++) { // Prepare complex sinusoid, rotation operator complex<double> W = polar(1.0, -2 * M_PI * i / N); // Butterfly computation complex<double> f_tmp = f[N * j + i] - f[N * j + n + i]; f[N * j + i] += f[N * j + n + i]; f[N * j + n + i] = W * f_tmp; } } } // Bit-reversal permutation int i = 0; for (int j = 1; j < N_orig - 1; j++) { for (int k = N_orig >> 1; k > (i ^= k); k >>= 1); if (j < i) { complex<double> f_tmp = f[j]; f[j] = f[i]; f[i] = f_tmp; } } }
Arduino IDE (M5Stack, ESP32)版
Raspberry Pi版をもとにして，M5StackやESP32 DevKit Cで使えるよう，Arduino IDE向けに移植します。ほとんど共通していますが，浮動小数点複素数型で利用するライブラリが異なります。
Raspberry Pi版では#include <complex>とした上で，complex<double>型を用いています。一方，Arduino IDE版では，Rob Tillaartさんの下記のライブラリを用います。
このライブラリはArduino IDEのライブラリマネージャからインストールできます。しかし，GitHubのREADME.mdにあるように，ArduinoオリジナルのAVR用とは異なるESP32用のコンパイラでは問題が生じます。と言うのも，ESP32用を含む非AVR用のコンパイラでは「Complex」という語が予約されており，コンパイルが通らないためです。
Apparently the name "Complex" is already in use (reserved) by some non-AVR compilers so it won't include the Complex.h file. Problem seen on Due and Teensy3.5
そのため，GitHubのREADME.mdに記載の「Solution」*13に則り，ライブラリのフォルダ名，その中のComplex.h，Complex.cppの「Complex」を「CComplex」にリネームし，かつ，CComplex.cpp内で#include <Complex.h>となっている箇所を#include <CComplex.h>に修正します。
その後，表3のように置き換えれば，Raspberry Pi版をArduino IDE版として移植できます。
| Raspberry Pi版 | Arduino IDE版 | 備考 |
|---|---|---|
| #include <complex> | #include <CComplex.h> | |
| complex<double> | Complex | 実部・虚部はそれぞれdoubleではなくfloatとなる。 |
| W = polar(...); | Complex W; W.polar(...); | 極形式で初期化するコンストラクタがないため，宣言と値の設定の2行に分ける。 |
なお，文字コードの問題で，Raspberry Piで作成したソースコードをArduino IDEで利用しようとすると，コンパイルエラーが出ることがあります。問題が生じたら，外部のエディタなどで適宜文字コードを変更します。
この他，Arduino IDEではインデントのタブ幅がデフォルトで2スペースなので，お好みでタブ幅を調節して下さい*14。
動作検証
それではいよいよ，実際に動かしてみることにします。
テストプログラム
以上で示したPythonとC++のプログラムはDFTやFFTを実行する関数の部分のみです。そこで，入力信号f[i]を作ってそれぞれの関数に与えたり，それぞれの関数の実行時間を計測したり，結果を人間に分かり易く表示したり，といったテスト機能を作りました。長くなってしまうので，本記事にはそのまま貼り付けませんが，代わりにGitHubに置きましたので，ご興味がある方は試してみて下さい。
時間領域の入力信号
それぞれのプラットフォームで実装した各アルゴリズムの実行結果や実行時間を比較するため，時間領域の入力信号f[i]を統一しました*15。
図3にデータ点数N = 128の場合の波形を示します。この波形は前記事から例として使っているもので，
ただし，
です。1次（基本波），7次，12次，17次高調波成分をもっており，それぞれの振幅は0.8，0.6，0.4，0.3です。また，正相成分のみとなっています。
実行時間の測り方
DFTに対するFFTの威力を確かめるため，
- Python:
import timeとした上でtime.perf_counter関数にて - C++:
#include <chrono>とした上でchrono::system_clock::now関数にて*16
実行時間を計測します。また，実装に誤りがないか否か確認するため，numpy.fft.fftによる計算結果と比較します。もちろん，numpy.fft.fftの実行時間も測ります。
Raspberry Pi，Windowsパソコン，AndroidスマートフォンはOSがさまざまなタスクを走らせている中でDFTやFFTを実行するので，システムの負荷によって，毎回同じ計算時間とはなりません。そこで，同じ計算を50回繰り返した平均値を計算するようにしました。一方，M5StackやNumWorksの関数電卓はベアメタル（OSなし*17）であるため，毎回（ほぼ）同じ計算時間となります。
それぞれのプラットフォームでの実行の様子
ここでは各プラットフォームでテストプログラムを実行しているイメージを示します。
Raspberry Pi
まずはRaspberry Piです。Raspberry Pi 3 Model A+と同Zero WHで同じプログラムを動かして検証してみたいと思います。
図4にRaspberry Pi 3 Model A+でPythonプログラムを実行している様子を示します。いわゆるヘッドレス運用であるため，ターミナルエミュレータRLoginからSSHでログインしています。DFT，FFTの各Typeとも誤差の範囲*18で同じフーリエ係数（スペクトル）を算出することを確認しました。
同様に，C++での実行の様子を図5に示します。
Windowsパソコン
筆者はWindowsパソコンにAnaconda3をインストールし，JupyterLabとSpyderをPython開発環境として使っています*19。Raspberry Piで動作するPythonプログラムはそのままWindows上のPythonでも（OSに依存する機能を使っていない限り）そのまま動作します。図6に，Spyderを使って実行している様子を示します。
筆者はいわゆる「Sandyおじさん」であり，いまだにIntel Core i7-2600Kを常用しています…💧 最新のPCであればもう少し快適なのかもしれませんね。なお，筆者の怠惰により，Windows上でC++をコンパイル・実行できる環境が整っていないため，今回はPythonのみとなります。
Androidスマートフォン(Pydroid 3)
AndroidにおけるPython開発環境として，Pydroid 3というアプリがあります*20。もちろんGoogle Playストアから簡単にインストールできます。スマホでPythonが動くので，ちょっとしたことを試すのにとても面白いですし，pipでいろいろな追加パッケージをインストールすれば本格的な使い方もできるようです*21。NumPyもpipからインストールできます。試したことはありませんが，GUIアプリも作れるそうです。
Raspberry Piで動いているPythonプログラムをそのまま走らせることができます*22。
図7に，筆者のXperia XZ1にてPydroid 3を使い，Pythonプログラムを編集および実行している様子を示します。ソースコードエディタの画面と実行結果を表示するターミナルの画面を1枚の画像として繋げました。
スマートフォンでこんなことができてしまうとは，いい時代になったものです。
M5Stack
Arduino IDEを用いてC++で実装したDFT，FFTの各アルゴリズムをテストします。M5Stackには液晶画面が付いているので，せっかくですのでRaspberry PiやPydroidのターミナル（ハッカーが触る黒い画面操作❗）のようなテキストをシリアル通信でパソコンに送るのではなく，グラフを描いてみることにしました。液晶画面のドット数は320 × 240なので，データ点数はデフォルトでN = 256とします。
 (a) DFT |  (b) FFT Type I |
 (c) FFT Type II |  (d) FFT Type III |
図6に実行している様子を示します。図8(a) ~ (d)は4つとも同じに見えますが，それぞれDFT，FFT Type I，II，IIIを実行した様子です。計算が正しければ，どれも同じスペクトルを表示するはずです。画面の下端に実行時間を表示しています。
NumWorksの関数電卓
フランスNumWorks社の関数電卓はリーズナブルな価格で入手できるパワフルなグラフ関数電卓です。本ブログでも過去にナイキスト線図を描く方法について書いています。
NumWorksの関数電卓にはMicroPythonにてプログラミングできる機能もあります。
今回は，このMicroPythonを用いて，DFT，FFTのアルゴリズムを実装しました。なお，NumPyは残念ながらありませんので，cmath版のDFTとFFT Type IIIのみを実装します。
NumWorksの関数電卓本体のボタンをポチポチ押してプログラミングすることももちろんできますが，あまりにも効率が悪いので，図9のように，NumWorksのウェブサイトにログインして，オンラインのエディタとシミュレータで動作を確認してから，USBケーブルで電卓本体に転送します。
NumWorksの関数電卓のMicroPythonの実行環境にはメモリが少ないらしく，データ点数N = 128までしか実行できませんでした。図10のように，実行時間を表示しましたが，データ点数N = 64までであれば，Matplotlibでスペクトルをプロットすることもできます。それにしても，電卓でFFTできるとはなかなか面白いです。
実行時間の比較
それぞれのプラットフォームにおけるDFT，FFTの各アルゴリズムの実行時間を比較してみました。結果は表3をご覧下さい。Raspberry PiのC++に関しては，データ点数N = 2,048の際に，コンパイラg++の最適化オプションを変えた場合ついても比較してみました。
う～ん，数字がチカチカして見づらいですが，4コアであるRaspberry Pi 3 Model A+と1コアである同Zero WHの計算パワーの違いが見て取れますね。また，データ点数Nの増加に伴って，cmath版のDFTでは爆発的に計算時間が伸びていくことが分かりますね。AndroidスマートフォンであるXperia XZ1の計算パワーは流石にWindows PCには及ばないものの，Raspberry Pi 3 Model A+よりは上回っています。しかしこれも，ラズパイ4になった場合にはどうなるか分かりませんね。NumWorksの関数電卓は，残念ながら桁違いに遅くなっています*23。
これだけでは分かり難いので，図11にグラフを描きました。列がすべて埋まっているRaspberry Pi 3 Model A+を代表として，データ点数Nと各実装の計算時間をプロットしています。数値の幅が広いので両対数グラフにしていますが，横軸の基数は2，縦軸の基数は10です。なお，C++のコンパイラg++には特に最適化オプションを与えておりません。
Python-cmath版のDFTが圧倒的に長い計算時間を誇り，トップを独走しております。また，逆に高速な方に注目すると，numpy.fft.fft関数の圧倒的なパワーを実感できます。それら2つの間に各アルゴリズムの実装が分布していますが，概ね3つのグループに分けられます。
グループ1は，Python-NumPy版とC++版のDFTです。何とわずかですが，Python-NumPy版の方がC++版よりも高速に計算しているようです。これは筆者のC++プログラムの実装が拙い（マルチスレッドやマルチプロセッシングを活用していない，など）ため，NumPyで積和演算した方が高速だったからではないかと推測します。
次に高速なグループ2には，Python-NumPy版とPython-cmath版のFFTの各Typeが分布しています。ここではcmath版がかなり健闘していることが分かりますね。いやむしろ，同じFFT Type IIIであれば，NumPyよりもcmathが優っています。この原因はよく分かりませんが，仮説としては，ndarrayとリストでは，スライスの扱い方が異なり，場合によってはリストの方が高速に取り扱えるからではないかと推測します。
さらに高速なグループ3には，C++でのFFTの各Typeが分布しています（Type IとIIIの折れ線グラフがほとんど重なっています）。
PythonとC++のFFTの各Typeを比較すると，面白いことに，それぞれの中で再帰呼び出しを使わないType IIIのランクが異なります。Pythonでは，Type IIIがType I, IIよりも顕著に高速になっています。特にデータ点数Nが小さい場合に違いが目立ちますね。一方，C++では各Typeで大きな差はないようです。おそらくですが，Pythonでは関数の呼び出し時に，呼び出し元の情報をスタックフレームオブジェクトとして生成する（と筆者は理解しています）のですが，これに要する計算コストがC++の関数呼び出しの際のそれよりも大きいためではないかと推測します。
なかなかおもしろい結果が得られたのではないかと思っております。が，しかし…Pythonでフーリエ解析が必要になったとき，MicroPythonなどの特殊な環境でない限り，自分でFFTを実装することは馬鹿げていると実感しました。素直にnumpy.fft.fft関数を呼ぶべきですね…😅
まとめ
以上，本記事ではPythonとC++でDFTやFFTのアルゴリズムを実装し，様々なプラットフォームで実行してみました。どうも今回はプログラミングのブログのようになってしまいましたね…。しかし，筆者自身はプログラミングに関しては*24ド素人であるため，プログラミングにお詳しい方々の目からすれば，噴飯ものかもしれません😅
このところ，ぱわみのある記事をあまり書いていないので，そろそろLTspiceを用いたパワー系のお話でもできればと思っていますが，恐らく次回はRaspberry Piで地震計を作るお話になろうかと思います。
付録: スライスについて
NumPyの配列(ndarray)とPythonのリスト
本記事では，NumPyの配列演算を活用したFFTとcmathのみを用いたFFTを作りました。cmath版はDFTの他はFFT Type IIIのみとなっていますが，ここではその理由を説明しながら，NumPyのndarray（配列）とPythonのリストにおけるスライスの挙動の違いについて述べたいと思います。
cmath版では入力信号f[i]をリストに格納するのですが，再帰呼び出しを用いるFFT Type I, IIを作ってみると，なぜか誤った結果ばかりが返ってきました。理由を探っていくと，スライスの挙動が異なるためであると判りました。
次の2つのコードを見てみましょう。コードA-1はNumPyのndarrayの挙動を表したものです。ndarrayとしてaを定義し，そのスライスをbに代入しています。(iii)の箇所でbの先頭要素b[0]に新たに6を代入して書き換えています。(v)の箇所を実行すると，bはもともとa[2:4]であったことが引き継がれており，a[2]が6に書き換わっていることが分かります。これは，NumPyのスライスは元のndarrayの「ビュー」であり元のndarrayとメモリを共有しているためです（ただし，aとbのオブジェクトIDは異なっています）。
# Code A-1, NumPy ndarray # (i) Define ndarray and its slice a = np.array([1, 2, 3, 4, 5]) b = a[2:4] # (ii) Show content of b print(b) # [3 4] # (iii) Change first element of b b[0] = 6 # (iv) Show content of b print(b) # [6 4] # (v) Show content of a print(a) # [1 2 6 4 5]
一方，コードA-2のように，同じことをPythonのリストで試してみると，(v)の箇所の結果の通り，元のaは変わっていません。リストをスライスすると，新たなオブジェクトが生成され，それは元のリストとは切り離されてしまうためです。
# Code A-2, Python list # (i) Define list and its slice a = [1, 2, 3, 4, 5] b = a[2:4] # (ii) Show content of b print(b) # [3, 4] # (iii) Change first element of b b[0] = 6 # (iv) Show content of b print(b) # [6, 4] # (v) Show content of a print(a) # [1, 2, 3, 4, 5]
上記を踏まえて，コードA-3のようにリストとndarrayのスライスを関数の引数として与えるとどうなるか試してみます。Pythonの関数では，ミュータブルなオブジェクトを引数として与えると，いわゆる「参照渡し」となります。したがって，ndarrayもリストも関数の中で要素を書き換えると，呼び出し元でもその要素が変わります（CやC++でポインタが与えられたような挙動ですね）。しかし，スライスが与えられると，リストでは厄介なことに，書き換えられません。
# Code A-3, Slices as function arguments # (i) Define function to change first element of argument def func(f): f[0] = 100 # (ii) Define ndarray and list a = np.array([1, 2, 3, 4, 5]) l = [1, 2, 3, 4, 5] # (iii) Apply function to ndarray's slice and see influence func(a[2:4]) print(a) # [1 2 100 4 5] # (iv) Apply function to list's slice and see influence func(l[2:4]) print(l) # [1, 2, 3, 4, 5]
そうなのです。(iv)の箇所のように関数にリストのスライスを与えても，スライスを作った時点で元のリストとは異なる別オブジェクトが生成されてしまうので，関数の中でその要素を書き換えても，元のリストには影響しません。
イ誌リスト1，2に基づくFFT Type I, IIでは，バタフライ演算の結果得られる入力信号を前半と後半に分け，つまり，スライスしてからそれぞれを引数として自身を再帰的に呼び出します。しかし，リストのスライスでは，関数の中で要素を書き換えても元のリストに影響しないため，このやり方では誤った計算結果となってしまいます。関数に戻り値を持たせるなど，工夫次第で回避できそうな気もしますが，筆者の怠惰ゆえに再帰呼び出しを用いないFFT Type IIIのみを実装することにしました。
全要素への代入
DFTとFFT Type Iでは，
# Overwrite results to input
f[:] = F
という箇所が登場します。最後の行ではスライスを用いてfの全要素を指定しています。ここで，スライスを用いず，f = Fとしてはいけない理由ですが，このようにしてしまうと，関数の中で変数fが新たにオブジェクトFを指すようになってしまいます。呼び出し元のfが指していたオブジェクトは関数の中ではどこからも参照されなくなってしまい，書き換えられません。関数から出れば，呼び出し元のfはもともとのオブジェクトを指しています。したがって，呼び出し元のfを書き換えたいという意図は叶いません。
一方，f[:]= Fとして全要素を表すスライスを左辺に置けば，fの各々の要素にFの各々の要素が代入され，呼び出し元のfも書き換わります。以下のコードA-4をご覧下さい。
# Code A-4, Assingment to slice for all elements # (i) Assign new ndarray to argument def func1(f): t = np.array([8, 9, 10]) f = t # (ii) Assign new ndarray to of all-element slice of argument def func2(f): t = np.array([8, 9, 10]) f[:] = t # (iii) Define list and ndarray a = [1, 2, 3] b = np.array(a) # (iv) Apply func1 to both func1(a) print(a) # [1, 2, 3] func1(b) print(b) # [1 2 3] # (v) Apply fun2 to both func2(a) print(a) # [8, 9, 10] func2(b) print(b) # [8 9 10]
これを実行すると，関数func1では呼び出し元のリストやndarrayは書き換わらず，func2では書き換わります。
このように，Pythonのスライスという機能は大変便利ですが，その実を理解するのは少し厄介です。リストやndarrayの挙動を詳しく理解しておく必要がありそうですね。
*1:プログラミングにおけるfunctionを「関数」と呼ぶのって，やっぱり変ですよね…。引数も戻り値も数値で，計算以外の処理（画面の表示などの入出力）をしないならば納得ですが…。慣れとは怖いもので，違和感を感じなくなっていました。
*2:動画1，図1ではNucleoで作った信号の周波数や振幅をきちんと制御していませんので，機会を見て作り直したいと考えています。
*3:しかも，データ点数は2のべき乗でなくても構いません。
*4:本記事ではndarrayを「配列」とも呼びます。
*5:複素正弦波(complex sinusoid)とも呼ぶことにします
*6:イ誌リスト1ではFではなく実部がYr，虚部がYiという名前でした。
*7:ビット反転並べ替えの箇所においては，筆者のスキルではNumPy（もしくはPythonそのもの）をうまく活用できませんでした。
*8:NumPyのndarrayとPythonのリストではスライスの挙動が異なるためです。
*9:Wikipediaの記事をはじめ，筆者が見つけた解説記事の大半では再帰呼び出しではなく，この多重ループでのプログラム例を説明していますが，これは筆者にとって理解の妨げになりました…。
*10:もう一工夫すれば計算できたのかもしれません。
*11:関数に渡す前は，どのように要素数の情報を保持しているのだろうか…？ 例えばint a[] = {1, 2, 3, 4}とした場合，sizeof(a) / sizeof(a[0])を計算すると，きちんと4が得られます。どうやってsizeof(a)を得ているのだろうか…？
*12:complex::コンストラクタ - cpprefjp C++日本語リファレンス
*13:と言うよりworkaroundですね。
*14:メニューバーから「ツール」➔「自動整形」と操作するか，CTRL + Tを押せば，2スペースに自動調整され（てしまい）ます。
*15:ただし，M5Stackでは実部と虚部がfloat型(32 bit)，他ではdouble型（Pythonの「float型」はC/C++のdouble型(64 bit)に相当）となっています。
*16:ただしArduino IDEではmicros関数を用います
*17:…と言ってもFreeRTOSはあります。
*18:筆者は数値計算やIEEE 754に明るくないので，アルゴリズムと誤差の評価が不適切かもしれません。
*19:南先生の「Pythonによる制御工学入門」をきっかけとして使い始めました。
*20:ロシア製アプリのようです。
*21:とりあえずMatplotlibは動きました。何とpigpioも動き，ネットワーク上のRaspberry PiのGPIOを操作できます。
*22:ただし，コマンドライン引数の与え方が分からないため，プログラム本文中にてデータ点数を指定しています。
*23:テストプログラムに問題があるのかもしれません。
*24:いや，他の分野についても
