晩秋から年末へと時が移ろいつつありますが，いかがお過ごしでしょうか？ 今朝は筆者の住む街でも氷点下を記録しました。寒くなって参りましたね。今回も残念ながらぱわみのある記事ではありませんが，お付き合い下さい。
さて，前々々回*1，前々回，前回を通じて，Raspberry Piを用いたマイ地震計を作りました。これは，図1のようにブレッドボード上にモーションセンサMPU6050とOLEDディスプレイを置いたプロトタイプであったため，今回はこれを基板，すなわちRaspberry Piの“HAT”にします*2。
早速ですが，図2に完成イメージを示します。
モーションセンサMPU6050とOLEDディスプレイの各モジュールは後で取り外せるように，1列ピンソケットを介して実装しました。また，前回までのプロトタイプにはなかったタクトスイッチを2個追加しました。いずれ，このタクトスイッチで画面を切り替えて波形を表示するなど，工夫していきたいと思います。
続きを読むRaspberry Piによるマイ地震計を作るため，前々回にて計測震度計算プログラムを，前回にてタイマ割込みやマルチプロセッシングなどリアルタイム処理を作って参りました。今回は，Raspberry Piに加速度センサとOLEDディスプレイを接続し，プログラム全体を組み上げていきます。
前々回，前回挙げた開発課題を再掲します。
今回は，上記5，6，7にあたる処理をPythonにて作って行きます。また，加速度センサ（モーションセンサ）の分解能やノイズについて，初歩的な実験を行って確認しました。
続きを読む前回と今回で，Raspberry Piを用いたマイ地震計について書いております。前回は加速度の時系列データから計測震度を計算するPythonプログラムについて述べました。
前回挙げた開発課題を再掲します。
今回は，上記2，3，4にあたる処理をPythonにて作って行きます。これは，マイ地震計を作るために必要なリアルタイム処理です。
続きを読む世界の地震の10%以上が日本で発生している──。巷間よく言われていることですが，先日も東京23区と埼玉県で10年ぶりに震度5強を観測したばかりであり，地震について考えない日はないと言っても過言ではないでしょう。
地震が発生した際に，気象庁のホームページやNHKなどで各地の震度が発表されます。しかし，これを待たずして自作の地震計で震度が分かったら，さらに，震度が分かるだけでなく加速度の波形も記録できたら面白そうです。昨今は加速度センサが簡単に手に入るので*1，地震計というのは電子工作やプログラミングの練習にもちょうどいい題材になるのではないでしょうか？
今回，次回，次々回を合わせて，筆者がマイ地震計を作ろうと思い至った経緯や，開発にあたっての課題，Raspberry Piでの実現方法について書いていきます。今回はその第一弾として，加速度の時系列データから計測震度を求める計算方法について書いていこうと思います。
加速度センサから定周期でデータを読み出しながら震度計算するといったリアルタイム処理については次回に，加速度センサとそのI2C接続，プログラムの全体構成については次々回に譲ります。
いきなりですが，実際に起きた地震を捉えた観測例を示します。図1に，試作中のマイ地震計プロトタイプを示します。マルチプロセッシングを活用するので，4コアのCPUを持つRaspberry Pi 3 Model A+を使うことにし，ブレッドボード上にInvenSense社のモーションセンサMPU6050とOLEDディスプレイを置きました。いずれもI2CバスにてRaspberry Piに接続されますので，Raspberry PiのGPIO2 (SDA)とGPIO3 (SCL)に並列に接続しました（もちろん電源も必要です）。ジャンパワイヤの配線はいたってシンプルです。［追記: 2021-11-19］同じI2CバスにモーションセンサMPU6050とOLEDディスプレイを繋げると，加速度の定期的な読み出しのサイクルが乱される可能性があることが判りました。［追記: 2022-03-15］同じI2Cバスに接続してもOLEDディスプレイのためのフレームの隙間にMPU6050のためのフレームがほぼ10 ms周期で挿入されており，ほぼ問題ないという結論に至りました。詳しくは次々回の付録をご覧下さい。
試作中のPythonプログラムは10 ms毎に加速度をサンプリングしつつ，加速度データが300点集まった時点で（つまり，3秒毎に）短時間計測震度を計算します。
ある夜，寝る前にプログラムを走らせておいたところ，早朝に発生した地震による加速度が記録されており，また，加速度から震度も（単に「それらしい値」ですが）計算されておりました。図2に結果を示します。この地震でマイ地震計プロトタイプが計算した計測震度の最大値は2.8であり，震度階級では震度3に相当します。
これは，2021年11月8日03時08分頃に発生した茨城県南部を震源とするマグニチュード4.3の地震です*2。気象庁のデータによれば，筆者の住む街では震度2を記録しています。一方，マイ地震計プロトタイプでは震度3を計測しました。恐らくですが，2階に置いてあったため，家屋での共振によって大きめの加速度を観測したのではないかと推測します。
なお，筆者のマイ地震計プロトタイプは，震度1以下の地震では反応しないようにしました。と言うのも，次々回で詳述しますが，モーションセンサMPU6050のノイズが大きいため，全く揺れていない状態と震度1を区別できないためです。より高精度の加速度センサが手に入れば，震度1以下も測れるのではないかと思います。
*1:後から分かりましたが，地震計として本格的に使うためには，1 galを下回るような小さな加速度を精度よく測れる必要があり，趣味の電子工作で手に入るような加速度センサでは，本来は対応できないと考えます。
*2:https://www.data.jma.go.jp/multi/quake/quake_detail.html?eventID=20211108031241&lang=jp
11月に入っていよいよ秋が深まって参りましたが，いかがお過ごしでしょうか？ さて，前記事にて，高速フーリエ変換(Fast Fourier Transform, FFT)に対する筆者なりの理解の仕方をまとめました。周波数間引き(decimation-in-frequency)型のCooley-Tukeyアルゴリズムや，その中で要となるバタフライ演算について述べた他，α-β変換やd-q変換など，電気屋としての思考の道具との関連性について考えてみました。
その続きとなる本記事では，FFTを実現する関数*1をPythonとC++で実装した様子や，それぞれのアルゴリズムを
にて動かしてみた結果について書いていこうと思います。また，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に移植してみても良いかと思います。
*1:プログラミングにおけるfunctionを「関数」と呼ぶのって，やっぱり変ですよね…。引数も戻り値も数値で，計算以外の処理（画面の表示などの入出力）をしないならば納得ですが…。慣れとは怖いもので，違和感を感じなくなっていました。
*2:動画1，図1ではNucleoで作った信号の周波数や振幅をきちんと制御していませんので，機会を見て作り直したいと考えています。
高速フーリエ変換(Fast Fourier Transform, FFT)は，通常の離散フーリエ変換(Discrete Fourier Transform, DFT)に比べて非常に高速なアルゴリズムとして知られており，ランダウの記号を用いると，DFTの計算量がO(N2)であるところ，FFTではO(N log N)になるようです。このため，理学・工学の幅広い分野で応用されております。言わば，ディジタル信号処理の中の1つの大きな“柱”となっています。ただし，入力される時間領域のディジタル信号のデータ点数が2のべき乗（例えば64とか1,024とか4,096など）となっている必要があります*1。
図1に，FFTのうち，本記事で取り扱う周波数間引き型Cooley-Tukeyアルゴリズムの計算途中経過を可視化したGIFアニメーションを示します。Python (JupyterLab)とMatplotlibを用いて作成しました*2。青が実部，橙が虚部を表しています*3。
ここでは，時間領域波形f(θ)として，
を与えています*4。jは虚数単位です(j2 = −1)。また，θは図1の横軸全体で0から2πまで変化すると考えて下さい*5。FFTの結果，(1)式の通り，1次（基本波），7次，12次，17次にスペクトルが立っていることが分かります。データ点数はN = 64 (= 26)です。
ところで，筆者はかなり以前から電力・電気系の（自称）技術者（俗にいう「電気屋」）として電力系統の高調波(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の本質を探訪してみようと思い立ちました。
*1:後述しますが，実は3のべき乗でもそれ以上の自然数のべき乗でも，アルゴリズムは存在します。ただし若干効率は下がると思います。
*2:図1では，FFTの結果を見やすくするため，それぞれのバタフライ演算の後でデータを2で除算しています。
*3:最近，1つ目の要素を青，2つ目の要素を橙にするのが流行っていませんか？ MatplotlibもMicrosoft Excelもそうですよね。
*4:このデータは正相成分のみをもっています。
*5:ちょっと乱暴な書き方ですが…ご容赦下さい…。
*6:フーリエ級数展開なのかフーリエ変換なのか，はたまた連続なのか離散なのか，ぼやかして言うときに「フーリエ解析」という言い方は便利だと思っています。
*7:Microsoft Excelでも「分析ツール」でFFTが可能です。
*8:本ブログの筆者の数学力は驚くほど低いことが知られています。
電力系統の周波数は，言わずもがな，東日本では50 Hz，西日本では60 Hzです。しかし，いつもピッタリ50 Hzや60 Hzとなっているわけではなく，それぞれ50 Hz，60 Hzを中心としながらも，常に揺らいでいます。発電が消費よりも大きければ周波数は上昇し，逆に消費が発電よりも大きければ周波数は低下します。現状，水車や蒸気タービンで回転している同期発電機が大勢を占めていることから，負荷が重ければ減速し，軽ければ加速する──とも言い換えられますね*1。
もちろん，50 Hz，60 Hzから離れ過ぎると問題が生じるので，地域によって異なりますが，±0.2 ~ ±0.3 Hzに収まるよう（±0.1 Hzの滞在率が概ね99%以上*2），水車に流す水量や，ボイラに投入する燃料を調整したり，周波数を観測して特定の発電所を制御したり，時刻，曜日，季節に応じて運転・停止する発電所を選択したりといった措置が取られています。
諸外国には，電力系統の周波数をリアルタイムにインターネット上に公開しているTSO (transmission system operator)もあります。例えば，スウェーデンのSvenska Kraftnätなどがそうです*3。
しかし，日本ではこれは実施されていません。電力広域的運営推進機関(OCCTO)ホームページでも，地域間連系線（周波数変換所や直流送電を含む）の融通量を見ることはできますが，周波数は見ることができません。実際，どういう風に周波数が変動しているのか？
*1:太陽光や風力など，系統連系インバータを介した発電設備が増えてくるとこの挙動が変わってくるとして，その対策に関する研究が最近では盛んです。
*2:https://www.occto.or.jp/iinkai/chouseiryoku/jukyuchousei/2020/files/jukyu_shijyo_19_02_02.pdf