初版

M5Atom Matrix"だけ"で作る簡易震度測定システム

製作品

本記事のシステム概要 ==== 本記事で説明する簡易震度測定システムは、M5AtomMatrixに、ArduinoIDEを使ってプログラムを書き込むことで実現できます(ソースコードは本記事の最後)。 本記事で扱うシステムが実現する機能は下記の通りです。 - 起動してから計測した [計測震度](https://www.data.jma.go.jp/eqev/data/kyoshin/kaisetsu/calc_sindo.html)の最大値を表示する機能 　→最大震度の整数・小数部分をそれぞれ緑・青のLEDの点灯している数で表現します 　→何らかの異常があるとき、中央の列のLEDが赤に、異常がないときは白に点灯します - 異常な震度を計測した時(センサにぶつかるなど)、最大震度をリセットできる機能 　→ボタンを10秒程度押すと、最大震度が0にリセットされ、ボタンを押してから20～30秒後から再度最大震度を表示し始めます 本記事中では、以下の点に関する工夫について説明しました。 - [arduinoFFT](https://github.com/kosme/arduinoFFT)を活用した逆フーリエ変換の処理の実装 - FreeRTOSのマルチタスクおよびキューを活用した、加速度の一定周期での計測および震度の計算処理の実装 開発の背景 ==== [以下の記事の通り](https://elchika.com/article/1ccbfaa7-a19a-436a-b7d0-3d6adbbc9409/)、筆者は以前から自宅の震度を計測するシステムを作ろうと考えていました。 [M5AtomとRaspberryPiで自宅の震度測定システムを作る](https://elchika.com/article/1ccbfaa7-a19a-436a-b7d0-3d6adbbc9409/) https://elchika.com/article/1ccbfaa7-a19a-436a-b7d0-3d6adbbc9409/ 上記の記事は、大きな地震が起きた後、「とりあえず版でもいいのでとにかく早く震度が測定できるシステムを作って公開する」ということで作成しましたが、今回は本来作りたかった形式に近い形で震度計測システムを作ることができたので本記事を執筆することとしました。 震度計測の仕組み ==== [過去の記事](https://elchika.com/article/1ccbfaa7-a19a-436a-b7d0-3d6adbbc9409/)でも記載していますが、[気象庁のHP](https://www.data.jma.go.jp/eqev/data/kyoshin/kaisetsu/calc_sindo.html)によると、震度は下記の手順で求めらます。 1. 水平2方向+上下方向、合わせて3方向の加速度を測定する →[気象庁のHPに上がっている加速度のデータ](https://www.data.jma.go.jp/eqev/data/kyoshin/jishin/index.html)は、加速度は100Hzで記録されています。 1. 加速度データをフーリエ変換し、フィルタをかけ、逆フーリエ変換をかける 1. フィルター処理済みの３成分の波形をベクトル的に合成する 1. 合成した加速度をソートし、上から0.3秒目のデータの加速度をaとしたとき、以下の数式で求められる値Iの小数第3位を四捨五入し、小数第2位を切り捨てる $$I=2log_{10}a+0.94$$ 例えば、I=5.195だったとき、小数第3位を四捨五入すると5.20となり、計測震度は5.2となる。一方、I=5.194ならば小数第3位を四捨五入すると5.19となり、計測震度は5.1となることに注意が必要です。 システムの検討 ==== 今回作成する震度計測システムは、以下の点を考慮して作成することとしました。 - できるだけ構成要素を少なくする 　- 同一のシステムを構築しやすくする 　　→本記事の読者が簡単にまねできる 　　→地震被害が発生しやすい場所ごとに設置し、局所的な震度を計測できる 　　　(水槽や本棚など) 　- 設置しやすくする 　　→電源の取得方法や構成要素の固定方法等の検討項目を減らすことができる - 24時間365日つけっぱなしでも罪悪感の沸かない装置を使用する 作成する震度計測システムの要求機能は下記の通りとすることとしました。 - 震度を測定する機能：(必須) 　-震度の妥当性を検証する機能(加速度の生値を後から確認する機能)：(今後の課題①) - 計算した最大震度を表示する機能：(必須) 　-外出先など、離れた場所で過去/最大の震度を確認できる機能：(今後の課題②) - 異常な震度を計測した時(センサにぶつかる等)、最大震度をリセットできる機能：(必須) 以上の条件を検討した結果、今回はM5AtomMatrixを使用することとしました。理由は以下の3点です。(最近発売されたATOM S3の方が、液晶があってリッチな情報表示ができるのでいいかも。AtomMatrixと値段変わらないし。) - これ単体で加速度の測定・演算処理・MatrixLED(25個)での表示が可能であり、要求機能の必須機能を実現できる見込みがあること - 自宅にWiFiによるインターネット接続環境があれば、M5AtomMatrixのWiFi接続機能によりデータの送信が可能であること(今後の課題①,※1) - [ATOMIC TFカードキット](https://www.switch-science.com/catalog/6475/)を使用すれば加速度データの生値を記録できること(今後の課題②,※2) ※1：筆者が引っ越したばかりで、現状WiFiによるインターネット接続環境が常設されていないため、本記事執筆時点での実装はあきらめた。 ※2：筆者の手元にも[ATOMIC TFカードキット](https://www.switch-science.com/catalog/6475/)はあるが、本記事のシステムの必須機能ではないため、記事の分かりやすさを優先して本記事のシステムでは扱わないこととした。なお、M5Stackで正確な周期で加速度を測定・記録するためには、マルチタスク処理を実装する必要があるなど、少し難易度が高いが、[参考となるページ](https://westgate-lab.hatenablog.com/entry/2020/03/28/113621)は調べればそれなりにヒットする。 難しいポイントと解決方法 ==== 上記の機能を実現する上で、実装上難しいポイントは「課題①：フーリエ変換・逆変換を行う」「課題②：100Hz(一定周期)で加速度を測定する」あたりです。それぞれポイントを説明します。 課題①：フーリエ変換・逆変換を行う ==== フーリエ変換は、自分で実装するととても大変&実行速度がかかってしまうため、[下記のライブラリ(arduinoFFT)](https://github.com/kosme/arduinoFFT)を使用することとしました。 [https://github.com/kosme/arduinoFFT](https://github.com/kosme/arduinoFFT) arduinoFFTの機能について調べると、高速フーリエ変換等の基本的な機能が備わっています。しかしながら、一般的にあまり用途がないためか、arduinoFFTを用いても、逆フーリエ変換の機能は実装されていない(？※)ようです。そこで本システムの開発にあたり、逆フーリエ変換の機能を実装する必要がありました。 ※：実装されているような気もするが、あまり情報がない・・・ 色々と調べてみた結果、Arduinoで逆フーリエ変換をしている例は全く見つからず、また逆フーリエ変換を実装した事例も見つかりませんでしたが、結局[WikiPediaの記事](https://ja.wikipedia.org/wiki/高速フーリエ変換)を参考に実装してみることにしました。 [WikiPediaの「高速フーリエ変換」のページ](https://ja.wikipedia.org/wiki/高速フーリエ変換)を参考にした逆フーリエ変換の実装の方針は下記の通りです。なお、加速度データなどの時系列データを$f(t)$、$f(t)$をフーリエ変換したデータを$F( \omega )$とし、$F( \omega )$の実部を$Re( \omega )$、虚部を$Im( \omega )$、虚数単位を$i$としたとき、下式を満たすものとします。 ### 前提条件 $F( \omega )=Re( \omega )+i \times Im( \omega )$ ### 実装の方針 1. $F( \omega )$の共役な複素数$\overline{F( \omega )}$を求める。なお、共役な複素数とは、ある複素数の虚部の符号を反転したものであり、下式を満たす。 $\overline{F( \omega )}=Re( \omega )-i \times Im( \omega )$ 1. $\overline{F( \omega )}$にフーリエ変換処理を施した結果の実部を取り出す。 なお、上記はフーリエ変換するデータ$f(t)$が実数の場合のみ成立します。また、2の手順後のデータは、FFTの処理の実装法によっては、データ数で除算するなどの処理を施して大きさを調整する必要がある場合があります。 ### 実装結果 以上を踏まえて実装した、逆フーリエ変換を行う関数は下記の通りです。ただし、本ソースコードのコメント中にもあるように、逆フーリエ変換を施した後のデータの順番についてはおかしくなっていることが確認できていますが、本記事の震度測定システムではデータの順番は重要ではないため無視しています(XYZ軸方向の加速度データの時間的対応さえ取れていれば問題ない)。 ```arduino:逆フーリエ変換をするための処理 #include "arduinoFFT.h"//要インストール #define SAMPLING_FREQUENCY 100//サンプリング周波数[Hz] #define FFT_SAMPLES 1024 double vReal[FFT_SAMPLES],vImag[FFT_SAMPLES];//FFT用の変数(実部と虚部) arduinoFFT FFT = arduinoFFT(vReal, vImag, FFT_SAMPLES, SAMPLING_FREQUENCY);//FFT用のオブジェクト //データを代入してフーリエ変換する関数 void FFT_IFFT(double *data) { //データの代入 for (int i = 0; i < FFT_SAMPLES; i++) { vReal[i] = data[i]; vImag[i] = 0; } FFT.Compute(FFT_FORWARD);//FFTをかける } //逆フーリエ変換をする関数 void FFT_IFFT(void) { for (int i = 0; i < FFT_SAMPLES; i++) { vImag[i] = -vImag[i];//共役な複素数にする=複素数成分の符号を反転する vReal[i] = vReal[i]; } //TODO：なんでかわからないけど、FFT_REVERSEにしないと(FFT_FORWARDだと)おかしな値になる。 //FFT_REVERSEにすると、大きさは正しそうだ(から問題ない)が、なぜかIFFT後のデータの順番が逆になる。 FFT.Compute(FFT_REVERSE);//FFTをかける } ``` 課題②：100Hz(一定周期)で加速度を測定する ==== 100Hzの一定周期でデータを取得するだけであれば、下記のようにすればOKです。 1. 現在時刻(マイクロ秒単位)を取得 1. データ取得処理の実施 1. 現在時刻を取得し、データ取得処理にかかった時間(マイクロ秒単位)を計算 1. 一定周期になるよう、マイクロ秒単位でスリープ ```arduino:処理を一定周期で行う例 #define SAMPLING_FREQUENCY 100//サンプリング周波数[Hz] void loop() { long t=micros();//経過時間計測用の変数 M5.IMU.getAccelData(&accX, &accY, &accZ); delayMicroseconds(1000/SAMPLING_FREQUENCY*1000 - (micros() - t));//この処理で10.000ms待つ } ``` しかしながら、本システムのように、一定周期でデータを取得しつつ、データがある程度たまったら計算処理を行う必要がある場合、データ取得サイクルに係る時間内に計算が終わらない可能性を考慮する必要があります。今回のシステムでいうと、100Hzで加速度を取得するためには、10ms以内にデータ取得等の処理を終わらせる必要がありますが、フーリエ変換等を用いる震度計算処理は、データ数によっては10ms以内に終わらない可能性が考えられます。 よって、本システムのように、一定周期で実行される必要のある処理と、時間がかかる処理を行う場合、マルチタスク処理で実装することが望ましいと思われます。 本システムにおいては、[下記のページ](https://lang-ship.com/blog/work/m5stickc-esp32-multi-task/)を参考に、マルチタスクおよびキューを用いて実装することとしました。キューを用いることで、加速度データが一定量(本記事のソースコード中では1024個＝10.24秒分のデータ)溜まったタイミングで震度計算処理を実行することができます。 [https://lang-ship.com/blog/work/m5stickc-esp32-multi-task/](https://lang-ship.com/blog/work/m5stickc-esp32-multi-task/) なお、計算震度はデータが一定数溜まるごとに計算しますが、一度に処理するデータ数は下記の要因を考慮して決定する必要があります。 - データを増やし過ぎてはならない理由 - データ数はO(N)、フーリエ変換処理に係る時間がO(NlogN)で増加する →データ数を増やしていくと、いずれデータがたまる時間よりもデータを処理する時間の方が長くなり、処理しきれなくなる - 加速度データおよび処理結果を保持する変数に応じたスタックサイズを確保する必要がある →データ数を増やし過ぎるとメモリが足りなくなる - データを減らし過ぎてはならない理由 - 震度は取得したベクトル合成データの大きい方から0.3秒目のデータを使用するため、計算に使用するデータ取得期間が短すぎると実際の震度よりも小さい値となる 本記事のシステムでは、以下の検討の元、データが1024個溜まる(10.24秒経過する)度に震度を計算することとしました。 - 高速フーリエ変換(FFT)を適用するデータの個数は、2のべき乗である必要がある - 同一の加速度データを用いて、データ数$2^4$～$2^{14}$個程度で実験的に計算したところ、$2^{10}$個以上であればそれなりに正確な震度(震度5強で震度の誤差0.1程度)が測定できた 作成する手順 === 以上を踏まえ、本システムの作成方法は下記の2ステップのみです。なお、下記の手順はArduinoIDEを用いてM5AtomMatrixを開発できる環境を把握している人向けです。ご不明な点があればコメントへ質問の書き込みをお願いします。 ①プログラムを作成する前に、ArduinoIDEの「ツール」→「ライブラリを管理」から「arduinoFFT」をインストールする ②ArduinoIDEの「ファイル」→「新規ファイル」で新しいファイルを作成し、下記のソースコードをコピペし、書き込む ソースコード ==== ```arduino:shindokei.ino //加速度値から震度を計算するシステム。 //★震度の計算方法の詳細は下記のページを参照 //https://www.data.jma.go.jp/eqev/data/kyoshin/kaisetsu/calc_sindo.html //★マルチタスク処理やキューの送受信に関して参考にしたページ //https://lang-ship.com/blog/work/m5stickc-esp32-multi-task/ #include "M5Atom.h" #include "arduinoFFT.h"//要インストール using std::sqrt; using std::pow; using std::exp; #define SAMPLING_FREQUENCY 100//サンプリング周波数[Hz] //★震度計算に使用するデータのサイズを指定。2のべき乗(1024以下)に設定すること。 //★「FFT_SAMPLES/SAMPLING_FREQUENCY」秒間のデータを処理した後、 //加速度の大きい方から「0.3*SAMPLING_FREQUENCY」番目の値を元に震度を計算するため、 //「FFT_SAMPLES/SAMPLING_FREQUENCY」秒が「0.3秒」より十分大きくないと正確な震度が計算できない。 //★FFT_SAMPLESの値に応じてスタックサイズを大きく設定する必要がある。 //★FFTに係る時間のオーダーはO(NlogN)、データ取得に係る時間はO(N)なので、 //O(NlogN)>O(N)より、FFT_SAMPLESの値を増やし過ぎるとマイコンの計算能力を超える。 //☆以上を考慮の上で適切な値を設定すること。 #define FFT_SAMPLES 1024 //データ処理用の変数 //ローカル変数として定義したいが、スタックオーバーフローしたので外に出した(どうすべきかよくわからん) double vReal[FFT_SAMPLES],vImag[FFT_SAMPLES];//FFT用の変数(実部と虚部) arduinoFFT FFT = arduinoFFT(vReal, vImag, FFT_SAMPLES, SAMPLING_FREQUENCY);//FFT用のオブジェクト double Filter_const[FFT_SAMPLES];//FFTした後の値にかけるフィルタ値用の変数(震度計算に特有のフィルタ値) double X_tmp[FFT_SAMPLES],Y_tmp[FFT_SAMPLES],Z_tmp[FFT_SAMPLES];//フィルタ処理後の値を一時的に格納する変数 double result[FFT_SAMPLES];//ベクトル合成後の加速度を格納する変数 float accX = 0, accY = 0, accZ = 0;//加速度を格納する変数 //フーリエ変換後に適用するフィルタ値の計算を行う関数 void calc_Filter_const(void){ for(int i=0;i<FFT_SAMPLES;i++){ int i_tmp=i; if (i>FFT_SAMPLES/2){//ナイキスト周波数以降の値は折り返し i_tmp=FFT_SAMPLES-i; } if(i==0){//直流成分は除去する(後述のFCがZeroDivisionにならないよう、例外処理する) Filter_const[i]=0.0; } else{//気象庁の定義通りのフィルタを計算する double freq=i_tmp*SAMPLING_FREQUENCY/(double)FFT_SAMPLES; double FH=1/sqrt(1+0.694*pow(freq/10,2.0)+0.241*pow(freq/10,4.0)+0.0557*pow(freq/10,6.0)+0.009664*pow(freq/10,8.0)+0.00134*pow(freq/10,10.0)+0.000155*pow(freq/10,12.0)); double FL=sqrt(1-exp(-pow(freq/0.5,3.0))); double FC=1/sqrt(freq); Filter_const[i]=FH*FL*FC; } } } //加速度をフーリエ変換して、フィルタをかけて、逆変換する関数 void FFT_IFFT(double *data,double *res) { //加速度値の代入 for (int i = 0; i < FFT_SAMPLES; i++) { vReal[i] = data[i]; vImag[i] = 0; } FFT.Compute(FFT_FORWARD);//FFTをかける //IFFTの処理：FFTの結果を共役な複素数にして、再度FFTをかけるといいらしい for (int i = 0; i < FFT_SAMPLES; i++) { vImag[i] = -vImag[i]*Filter_const[i];//共役な複素数にする=複素数成分の符号を反転する vReal[i] = vReal[i]*Filter_const[i]; } //TODO：要検討：なんでかわからないけど、FFT_REVERSEにしないと(FFT_FORWARDだと)おかしな値になる //TODO：要検討：FFT_REVERSEにすると、値は正しそうだから問題ないが、なぜかIFFT後のデータの順番が逆になる FFT.Compute(FFT_REVERSE);//FFTをかける for (int i = 0; i < FFT_SAMPLES; i++) {//処理結果を格納する res[i]=vReal[i]; } } //qsort用の関数(降順) int compare_sort(const void *a, const void *b){ if (*(double*)a > *(double*)b){return -1;} else if (*(double*)a < *(double*)b){return 1;} else{return 0; } } //フィルタ処理した加速度値から震度を計算する関数 double calc_shindo(double *data_X,double *data_Y,double *data_Z){ //加速度データにフィルタ処理をかける FFT_IFFT(data_X,X_tmp); FFT_IFFT(data_Y,Y_tmp); FFT_IFFT(data_Z,Z_tmp); //ベクトル合成する for (int i = 0; i < FFT_SAMPLES; i++) { result[i]=sqrt(pow(X_tmp[i],2.0)+pow(Y_tmp[i],2.0)+pow(Z_tmp[i],2.0)); } //TODO：もうちょっと分かりやすいソートの手段がないか要検討 qsort(result, FFT_SAMPLES, sizeof(double), compare_sort);//降順に並び変える double shindo=std::log10(result[int(SAMPLING_FREQUENCY*0.3)])*2+0.94;//震度を計算する //小数第3位で四捨五入して、小数第2位以下を切り捨てる処理 //★shindo=3.00501…のときの例を示す //★基本的に小数第1位を四捨五入した方が楽なので、100倍して四捨五入して100で割ることで第3位を四捨五入する //☆100倍する //shindo*100=300.501 //☆小数第1位で四捨五入：+0.5してからdouble→intの変換 //単にdouble→intの変換をすると小数点以下切り捨てとなる。 //そこで+0.5してから型変換することで、四捨五入を実現する //(小数点以下が0.5以上ならば小数部分が1.0以上となり切り上がり、 //0.5未満なら0.99…となるため切り捨てられる) //int(shindo*100+0.5)=int(301.001)=301 //☆小数第2位を切り捨てる：int型で10で割る //ここまでで計算した値を10で割って、1の位の値を切り捨てる(最終的な小数第2位以下の値の切り捨てに相当) //(int)301/10=30 //☆最終的な震度の値と一致するように戻す：double型に変換して10で割る //double(30)/10=3.0 shindo=double(int(shindo*100+0.5)/10)/10.0; return shindo; } //キュー関連のあれこれ struct QueueData {//受け渡しデータ用の構造体 double X[FFT_SAMPLES]; double Y[FFT_SAMPLES]; double Z[FFT_SAMPLES]; }; QueueData queue_data; QueueHandle_t xQueue; TaskHandle_t calc_shindo_task_Handle; int queue_counter=0;//受信した加速度データの数を数えるカウンタ int test_counter=0;//テスト用カウンタ //震度の計算等を行うタスク void calc_task(void *pvParameters) { double max_shindo=0.0; bool reset_flag=true;//リセットボタンで再起動したときの衝撃で過大な震度がでるため、起動直後の1エポック分の震度は計算しない(最大値に反映しない) while (1) { QueueData queue_data; if (xQueueReceive(xQueue, &queue_data, 0) == pdTRUE) {// キューを受信 double shindo=calc_shindo(queue_data.X,queue_data.Y,queue_data.Z);//震度を計算する //最大加速度リセットボタンを押した際の処理 //★リセットフラグを用いることで、リセットボタンが押されていた間の震度は最大震度の比較対象外とする //(ボタン押しの衝撃で大きな震度となるため) if(reset_flag){//ボタンが押されていた間の次のエポックでフラグをFalseにする reset_flag=false; } else{//ボタンが押されていた間の次のエポックの震度は、最大震度の比較処理を行わない if (shindo>max_shindo){ max_shindo=shindo; } } //ボタンが押されていたら最大震度を0にして、リセットフラグをTrueにする //本当はこのタスク上でボタンの押し下げ判定をしない方がいい(別タスクにすべき) M5.update(); if(M5.Btn.wasPressed()){ reset_flag=true; max_shindo=0.0; } //最大震度の整数部分と小数第1位の数字を取得する int I,F;//整数部と小数部の数字を格納する変数 int tmp;//計算用の変数 //max_shindo=1.2だが、内部的には1.199…だったりすと、 //max_shindo/0.1=11.99…となってしまい、int型に変換すると11になってしまう(本来は12になってほしい)。 //そこで、max_shindo/0.1に、ほんの少しだけ値を追加することでこのような事態を回避する tmp=int(max_shindo/0.1+0.5);//震度の小数第2位以下が切り捨てられていなかった場合に備え、+0.5することとした I=tmp/10; F=tmp-I*10; //最大震度をLEDで表示する：整数部分の数を緑、小数部分を青で表示 M5.dis.clear(); for(int i=0;i<I;i++){ M5.dis.drawpix(i, 0x00ff00);//緑：震度の整数部分 } for(int i=0;i<F;i++){ M5.dis.drawpix(24-i, 0x0000ff);//青：震度の小数部分 } disp_status(0xffffff);//ステータス：白(正常) for(int i=7;i<11;i++){ M5.dis.drawpix(i, 0xffffff);//白：震度8以上になることはないので、8～10のLEDは白く表示 } Serial.printf("max:%.1lf,current:%.1lf\n",max_shindo,shindo); } delay(1); } } //真ん中の列のLEDでステータスを表示する関数 //色を指定すると真ん中の列のLEDが指定色で光る void disp_status(int color){ for (int i = 10; i < 15; i++){ M5.dis.drawpix(i, color); } } void setup() { M5.begin(true, false, true); if (M5.IMU.Init() != 0){ Serial.println("ImuInitError"); disp_status(0xff0000);//ステータス：赤(異常) } else{ disp_status(0x00ff00);//ステータス：緑(初期化中) } //フーリエ変換後に適用するフィルタ値の計算 calc_Filter_const(); //キュー作成 xQueue = xQueueCreate(2, sizeof(QueueData)); //Core0の優先度1で震度計算タスク起動 xTaskCreateUniversal( calc_task, // タスク関数 "calc_shindo", // タスク名 32768, // スタックサイズ NULL, // 引数 1, // 優先度(大きい方が高い) &calc_shindo_task_Handle, // タスクハンドル PRO_CPU_NUM // 実行するCPU(PRO_CPU_NUM or APP_CPU_NUM) ); delay(10); } //メインの処理(加速度データの取得) void loop() { long t=micros();//経過時間計測用の変数 if(queue_counter==FFT_SAMPLES){//データがFFT_SAMPLES個そろったらキューを送信 queue_counter=0; int ret = xQueueSend(xQueue, &queue_data, 0);//キューを送信 if (!ret) { //キュー送信失敗 Serial.println("QueueSendError"); M5.dis.clear(); disp_status(0xff0000);//ステータス：赤(異常) } } M5.IMU.getAccelData(&accX, &accY, &accZ); //MPU6886が出力する値の単位は[G]のため、1G=980galに直す queue_data.X[queue_counter]=accX*980; queue_data.Y[queue_counter]=accY*980; queue_data.Z[queue_counter]=accZ*980; queue_counter++;//キュー送信用のカウンタの加算 delayMicroseconds(1000/SAMPLING_FREQUENCY*1000 - (micros() - t));//この処理でSAMPLING_FREQUENCY[Hz]になるよう(10.000ms)待つ } ``` まとめ ==== 以上のようにして、M5AtomMatrixのみを用いて、簡易的に震度を計測するシステムを作成しました。なお、本記事で作成したシステムは、[以前作成した、Raspberry PiとM5AtomLiteを組み合わせたシステム](https://elchika.com/article/1ccbfaa7-a19a-436a-b7d0-3d6adbbc9409/)を用いて、震度を計算する処理の妥当性を検証しています。具体的には、以下の手順で検証しました。 ・気象庁のHPで公開されている加速度データを、震度を計算する関数に入力する ・Raspberry Pi(以前作成したもの)・M5AtomMatrix(今回作成したもの)で計算された震度と、元としたデータの震度(気象庁が計算したもの・真値)と、を比較する 以上の手順でいくつかのデータで確認したところ、いずれも計算震度が一致したため、それなりに妥当な計算処理をしていると結論付けました。 なお、記事中でも記載していますが、本記事は**最小限の構成で震度計を作成する**ことを目的としています。よって、例えば以下のような改良をするのもよいかと思います。 ・最近発売したATOM S3を使用して、過去の震度のグラフを表示する ・[M5AtomLiteとIMU Unitを用いて加速度を取得](https://elchika.com/article/c8b95b9a-3b09-4a3e-9e48-212a258fd917/)し、計算した震度をWiFi等を用いて別のデバイスに送信・表示する ・ATOMIC TFカードキットを使用して、生の加速度データや震度データを記録する