ESP32向けmicropythonで、FFTをする自作モジュールを作る

概要

ESP32用のMicroPythonで、FFTをする自作モジュールを作る。

背景と目的

前回、ESP32用のMicroPythonで、C実装の自作モジュールを作る環境が整った。今回は、この環境を使ってFFTをする自作モジュールを作る。なお、PythonでそのままFFTの計算をコーディングする手もあるが、これはPCですら非常に遅いので、マイコンでは全く使い物にならないと思われる。Cで実装したものを呼ぶという今回の形が有効だろう。

詳細

0.実施条件

前回と同じ。

1.FFTプログラムの実装

1.1 課題整理

プログラム実装にあたり、以下の部分がよくわからなかった。FFTをやる以上、扱うデータは配列形式になるはずなので、Pythonのarrayでもlistでもtupleでもいいのだが、C側とやり取りするにはどうしたらいいか、というところを知っておく必要がある。

• 引数として配列を渡す方法
• 戻り値として配列を返す方法

1.2 引数として配列を渡す方法

このあたりを中心にWeb上をいろいろ調べたのだが、体系的に情報がそろっているところが見つけられなかったので、どうにか自力でやってみた。必要な部分だけ抜き出す。どこかに、ドキュメントってあるのだろうか。。。 計算用のdouble配列として、実部y_real[N]、虚部y_imag[N]を用意。FFTなので、初期値として、y_realに入力信号そのもの、虚部に0を入れている。 mp_obj_get_array_fixed_n関数で、Pythonのリストx_objの各要素N個をmp_obj_t型にバラして並べ、その後ろのforでmp_obj_get_floatによってdouble型の数値を取り出す必要があるようだ。

STATIC mp_obj_t example_mcpy_fft(mp_obj_t x_obj, mp_obj_t N_obj) {

    int N = mp_obj_get_int(N_obj);

    double y_real[N], y_imag[N];
    mp_obj_t *x_p;
    mp_obj_get_array_fixed_n(x_obj, N, &x_p);
    for (i = 0; i < N; i++) {
        y_real[i] = mp_obj_get_float(x_p[i]);
        y_imag[i] = 0.0;
    }

1.3 戻り値として配列を返す方法

同様に、このあたりを手掛かりに、やってみた。基本的には逆でいいのだが、実部、虚部を返さないといけないので、ちょっと手順が多い。 最終的な戻り値は、[実部リスト, 虚部リスト]の2次元配列的なものにする。 まず、double型配列y_real、y_imagを、mp_obj_tの配列y_real_p、y_imag_pにmp_obj_new_float関数を使って入れていく。さらに、mp_obj_new_list関数を使ってy_real_p、y_imag_pを2要素の配列y_pとしてまとめる。最後に、y_pをmp_obj_new_listでPythonリストに変換して返す。

    mp_obj_t y_real_p[N], y_imag_p[N];

    for (i = 0; i < N; i++) {
        y_real_p[i] = mp_obj_new_float(y_real[i]);
        y_imag_p[i] = mp_obj_new_float(y_imag[i]);
    }
    mp_obj_t y_p[2] = {
        mp_obj_new_list(N, y_real_p),
        mp_obj_new_list(N, y_imag_p)
    };

    return mp_obj_new_list(2, y_p);

1.4 FFT

いわゆるバタフライ演算をただ書いただけなので、コードは省略。

2. 動作確認

作成したCモジュールを組み込んだファームウェアをビルドし、ESP32に書き込む。そして、以下のコードを実行。 なお、FFTにもかかわらずpython側から呼び出している関数名がexample.add_intsというナンセンスな名前なのは、名前を変えようとしてもCのコンパイルエラーが出て解決できなかったから。前回のサンプルを改造して作っているのだが、なぜか名前を変えられない。。。

import example
import math
import utime

# 入力信号作成
N = 8
x = [0.0] * N
for i in range(N):
    x[i] = math.cos(2 * math.pi * i / N)

t0 = utime.ticks_us()
y = example.add_ints(x, N) # FFTを実行
t1 = utime.ticks_us()

# 結果表示
for i in range(N):
    i , y[0][i], y[1][i]
"elapsed_time:", (t1-t0) / 1000, "msec"

結果は以下。長さ8で、1周期分のコサイン波を入力としているので、実部の基本波成分に値が立ち、その他と虚部はほぼ0。一応正しくできていそう。所要時間は、約29msec。

(0, 3.258414e-07, 0.0)
(1, 4.0, 3.085174e-07)
(2, 3.17865e-08, 1.192093e-07)
(3, -5.006516e-08, 1.972448e-07)
(4, -3.894144e-07, 0.0)
(5, -5.006516e-08, -1.972448e-07)
(6, 3.178651e-08, -1.192093e-07)
(7, 4.0, -3.085174e-07)
>>> "elapsed_time:", (t1-t0) / 1000, "msec"
('elapsed_time:', 29.873, 'msec')

また、Nを大きくしたときの挙動を確認。N=512で、約67msec。

:
(506, -7.816201e-08, -7.783505e-07)
(507, -1.293481e-06, 4.859487e-07)
(508, -3.251388e-07, -1.517361e-06)
(509, -2.192639e-06, 1.296839e-06)
(510, -4.54177e-06, -2.580137e-06)
(511, 256.0, -2.165856e-05)
>>>
>>> "elapsed_time:", (t1-t0) / 1000, "msec"
('elapsed_time:', 67.453, 'msec')

FFTの計算オーダーは

と言われるので、N=8で24、N=512で4608と192倍差があるはずだが2倍ちょっとにしかなっていない。つまり、FFTの計算そのものよりも、上に挙げたデータ受け渡しのオーバーヘッドの方が全然大きいのだろう。

さらに、Nを1024にしてみたところ、メモリ不足でESP32が再起動してしまった。メモリ使用量を工夫しないと、あまり大きなNを扱うのは難しそう。でも、ESP32で(MicroPythonではないが)画像処理をしている例もあるので、この程度で足りなくなるとも思えないのだが。後で時間があれば、改善してみたい。

まとめと今後の課題

ESP32用のMicroPythonで、FFTをするモジュールを作ることができた。