The Negligible Lab

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

ESP32マイコンで商用電源品質監視装置を作る

C++ Googleスプレッドシート ESP32 Arduino IDE ディジタル信号処理 電子工作

はじめに

本ブログも実質的に3年目に入りました。このところ半年以上も新しい記事を投稿できておりませんでしたが，一応，生存はしております。筆者が多忙を理由に怠惰ゆえに本ブログを放置していた間に，世界史を画するような大事件がいくつも起こりましたね…。そのような社会情勢から，今後，電力需給が逼迫する事態が頻発する可能性が高まってくるでしょう。電源品質監視装置に興味があるという方も増えてくるのではないかと思っています。

過去の試み ～浮遊静電容量を用いた周波数監視装置～

さて，本ブログでは以前に，浮遊静電容量を用いて電力系統の周波数を測る装置について書いております。

negligible.hatenablog.com

これは，M5Stackのアナログ入力端子に電位として浮いている電線や金属塊を繋げ，電力系統周波数と同じ周波数のいわゆる「ハムノイズ」を拾うことで商用電源の周波数を測ろうという試みでした*1。100 Vコンセントなどに接続することなく，完全に浮いた状態で周波数を測れるという面白さがありましたが，上記の記事にあるように，如何せん波形については正しく取得できないことから，この装置は周波数しか計測できません。やっぱり電圧波形そのものをきちんと取得して，実効値や高調波を観てみたい，つまり，「電力系統の息吹を感じたい」との思いから，変圧器を用いた新しい監視装置を作ろうと考えました。

完成イメージ

図1に完成した商用電源品質監視装置の写真を示します。

図1: 完成した商用電源品質監視装置

AC 100 Vという超高圧（？）を取り扱うことから，人が充電部に触ることができないよう，システム全体をアクリルケースに収めました。マイコンはM5Stackにも搭載されているEspressif Systems社のESP32で，秋月電子通商のモジュール(AE-ESP32-WROOM-32E-MINI)を使っています。ケース前面右上に240 × 240ドットの液晶ディスプレイを取り付け，周波数，実効値，波形，スペクトルを表示するようにしました。ディスプレイへの描画に関しては，らびやんさんのLovyanGFXを活用させて頂いております。

動画1に完成した装置の動作画面を示します。

youtu.be
動画1: 完成した商用電源品質監視装置の動作画面

スペクトルに関しては，本ブログでも以前に書きましたFFTのプログラムをそのまま活用しています*2。マイコンで自作FFT関数がリアルタイムでうねうね動いている──こりゃ超気持ちいい！

この装置は，計測した周波数，実効値，スペクトル（3, 5, 7, 9, 11次高調波含有率とTHD）を10秒毎にGoogle スプレッドシートにPOSTします。図2にGoogle スプレッドシートに作った監視画面，いわゆるダッシュボードを示します。これは今まさに本記事の執筆時点での様子です（早く寝ましょう💧）。

図2: Google スプレッドシートのダッシュボード

これまでに，2022年6月の電力需給逼迫に伴う周波数低め運用や，落雷に起因すると思われる電圧変動，また，宅内での家電機器の入切に伴う変動などを観測しております。

本記事の内容

本記事では，変圧器をいわゆるPT (potential transformer)として用いた電圧波形の検知方法に注目したESP32周りの回路構成について述べます。その後，周波数，電圧実効値，スペクトルの計算方法について書いていきます。さらに，Arduino IDEにて作成したC++でのプログラム，特に定周期割り込みなどについて説明します*3。また，例によって十字配線ユニバーサル基板を用いた実装について簡単に触れます。最後に，電力需給逼迫時や落雷発生時など，過去の特異な観測例をいくつかご紹介します。

目次

ハードウェア構成

本章では，AC 100 Vコンセントの電圧をESP32マイコンに取り込むための分圧ネットワークや，十字配線ユニバーサル基板を使った基板設計などについて述べます。今回のプロジェクトでは，電子回路・ハードウェアとしては極めて単純な作りになっています。厚めのアクリル板を使った透明ケースに収めてちょっとばかり見栄えにこだわったところがこれまで（自作スマートリモコンなどは基板だけでした）のとの違いになりますね。

ESP32マイコン周り

図3にシステム全体のハードウェア構成を示します。AC 100 Vの入力から電源ユニット（後述）にてDC 5 Vを得て秋月電子通商のESP32モジュール(AE-ESP32-WROOM-32E-MINI)に給電します。ESP32モジュールは内部に3端子レギュレータを持っており，ESP32および周辺のデバイスを動作させるための3.3 Vを作っています。また，240 × 240ドットのカラーLCDをSPIバスにてESP32マイコンに接続しています。

図3: システム全体のハードウェア構成

分圧ネットワーク

一方，AC 100 Vを分岐して100 V/6 Vの変圧器に接続します。なお，6 Vでなくても良いのですが，たまたま手元にあったものを利用しています*4。これはもともとは電源トランスとして作られたものですが，今回はこれをいわゆるPT (potential transformer)，つまり計器用変成器として扱います*5。

さて，ESP32マイコンのアナログ入力ピン*6に対しては当然ながら正の電圧(0 V ~ +3.3 V)しか入力できませんが，AC 6 Vからどのようにしてこの範囲のアナログ電圧を得るのでしょうか？ そう！ 変圧器の2次側（6 V側）はいわゆるフローティング電源として扱えるため，その一端を3.3 Vのちょうど半分，1.65 Vの電位に固定してしまえば良いのです。そのために，R1 = R2 = 330 Ωの抵抗器2本で3.3 Vの電源電圧を分圧し，変圧器2次側の一端に接続します。また，6 Vは大き過ぎるため，R3 = 1.8 kΩとR4 = 120 Ωの抵抗器で分圧します。また，気休め程度ですが，R1, R2, R3と並列に0.1 μFのコンデンサを接続してノイズ低減を図ります*7。

その結果，3.3 V / 2 = 1.65 Vを中心とする振幅0.53 Vの交流信号が得られます。電圧の範囲としては1.12 V ~ 2.18 V，ESP32の持っている12 bit A-D変換の結果としては1,391 ~ 2,707となります。ESP32マイコンのA-D変換器については，0 Vや3.3 Vに近い領域での線形性の悪評が高いですが，これを避けるよう，真ん中付近に小さめの振幅の信号として見えるように設計しました。図4に分圧ネットワークのLTspiceシミュレーション結果を示します。単純な回路ですが，一応，シミュレーションにて意図通り動作していることを確認しました。ここではESP32マイコンのアナログ入力端子の入力抵抗を1 MΩと想定しています（実際はもっと大きいと考えられます）。

図4: 分圧ネットワークのLTspice解析結果

なお，図3にあるような電源トランスにおいて「6 V」と謳っている場合，これは定格負荷を取った場合の値というのがよくある話です。今回は無負荷に近いため，これより高い電圧が見えている可能性があります。実際，測ってみると7.9 V（高い！）となっていました。さらに言えば，R1 ~ R4には誤差5%の炭素皮膜抵抗器を使っています（誤差1%のものが手元になかったため）。したがって，上記の計算とは合わない可能性が非常に高くなっています。まぁ，0 Vや+3.3 Vをはみ出ることはなく，真ん中(1.65 V)付近に信号があればよいという（テキトーな）設計です。最終的にはHIOKIのハンディテスタの表示している実効値と一致するように校正しました。

基板設計とハードウェア実装

まずはブレッドボード上でプロトタイピング，つまりは後述のソフトウェアが完全に動作するまでは基板化はせず，デバッグを繰り返していました。図5にプロトタイピングの様子を示します（ブレッドボードにAC 100 Vを印加する実績を解除しました）。

図5: ブレッドボードによるプロトタイピングの様子

この状態でソフトウェアをデバッグし，完成した（と思った）ところで，これまでと同様に十字配線ユニバーサル基板を用いて基板に実装しました。今回は最初からアクリルケースに収めるつもりでしたので，LCDについては基板に直接実装せず，アクリルケースにねじ止めし，ESPマイコンとはハーネス（と言うほどでもありませんが）を作って接続しました。同様に，100 V/6 Vの変圧器（図5の左端に写っているもの）も基板上には取り付けられないので，アクリルケースにねじ止めし，メイン基板まで電線とコネクタで接続します*8。

図6: 十字配線ユニバーサル基板の設計

図6に十字配線ユニバーサル基板の設計を示します。分圧ネットワークはESP32モジュールの下に入れてスペース確保を狙いました。また，AC 100 Vを十字配線ユニバーサル基板で取り扱うため，AC 100 V部の絶縁にはカットパターンを1本余計に入れ，かつ，低圧部との間にもカットパターンを引きました。恐らく製品としてはこれでは不十分だと思いますが，とりあえず3ヶ月ほどは連続安定稼働を続けています（いつか閃絡するのだろうか…）。

図7: AC 100 V-DC 5 V, 700 mA電源モジュール

www.amazon.co.jp

AC 100 VからDC 5 Vを得る電源モジュールには，Amazonで手に入れた小さなモジュールを使っています。図7に写真を示します。3個セットで800円を切るというコストパフォーマンスです。これも3ヶ月ほど連続安定稼働を続けています。

アクリルケースにはこれまたAmazonで入手したものを用いています。

www.amazon.co.jp

ドリルで慎重に孔をあけ，メイン基板とLCDを取り付けます。また，AC 100 Vのケーブルが通る孔もあけます。

サンプリングと信号処理

本章では，上記の分圧ネットワークによってESP32マイコンのGPIO35に入力された電圧信号をサンプリングし，周波数，実効値，スペクトルを計算する信号処理について概説します。

全体構成

図8にメインとなる信号処理部のブロック図を示します。信号処理部は主に312.5 μs毎（50 Hzの1周期あたり64回*9）にアナログ信号をサンプリングし記録する部分（タスクsampleVoltage()関数にて実行しており，以降，「サンプリング部」と称します）と，約50 ms毎に周波数，実効値，FFTを計算し，波形と共にLCDに表示する部分（loop()関数にて実行しており，以降，「表示処理部」と称します）に分けられます。また，図示していませんが，Google スプレッドシートにPOSTするタスクが別途動作しています。

図8: 信号処理部のブロック図

サンプリングした信号に対して，カットオフ周波数1 HzのディジタルHPF（2次バターワース特性）にてDC成分*10を取り除き，また，カットオフ周波数500 HzのディジタルLPF（同じく2次バターワース特性）にて高周波のノイズを取り除きます。HPFの後の信号とLPFの後の信号をグローバル変数として確保した配列（8周期分のリングバッファ）に順次格納します。また，LPFを通った後の信号から周波数と実効値を計算します*11。なお，LPFを通った後の信号の零クロスから周期を計算し，その逆数として周波数を得ます。

零クロスの検知方法と周波数の計算

周波数は周期の逆数として測ります。周期の計測のためには零クロスを検知する必要があります。

図9に零クロスの検知方法を示します。本プロジェクトでは立ち上がり零クロスのみを検知します。312.5 μs毎に訪れるサンプリングのたびに，サンプリングした電圧の前回の値が負で今回の値が正であれば零クロスがあったと見做します。これはノイズに弱いため，前述のようにカットオフ周波数500 HzのLPFを通った後の信号を使っています（フィルタを通っても周波数は変わらないため）。また，零クロスのあり得る範囲として窓(20 ± 2.5 ms)を設定しています（したがって，44.4 ~ 57.1 Hz（周期で言えば22.5 ms ~ 17.5 ms）が測定範囲となります）。零クロスを検知したら，

  • 前回の零クロスを検知してから今回の零クロスを検知するまでのサンプリング回数N
  • 零クロス検知の1サンプル前の信号の値vn
  • 零クロス検知時の信号の値vp

の3つをzeroCrossingというクラスのオブジェクト*12の配列（零クロス100回分のリングバッファ）に格納します。

図9: 零クロス検知方法

定周期サンプリングしているため，真の零クロスを検知していませんが，零クロス近傍での信号が直線と見なせるのであれば，図6中の拡大図のように，真の零クロスから零クロスを検知したサンプリングまでの時間Tzcを次式で推定できます（vn < 0であるため分母は正です。）。

T_{zc} = \displaystyle \frac{v_{p}}{v_{p} - v_{n}} T_{s} \tag{1}

M回分（M周期分，前述の通りM = 100です）の零クロスに要した時間TMを計算するにあたって，まず100周期分のそれぞれのNとサンプリング周期Tsの積の和を得ます。さらに，リングバッファ上の最初（0回目）の零クロスに関するTzcを加算し，最後（(M − 1)回目）の零クロスに関するTzcを減算します。数式で表せば次式のようになります。

T_{M} = \displaystyle \sum_{i = 0}^{M - 1} T_{s} N[i] + T_{zc, 0} - T_{zc, M - 1} \tag{2}

TMが得られれば，次式から周波数fを容易に計算できます。

f = \displaystyle \frac{M}{T_{M}} \tag{3}

実効値の計算

実効値Vrmsに関しては直接的にroot mean squareを計算しています。図5のLPFの後の信号をvLPFとすれば，

V_{\mathrm{rms}} = \displaystyle \sqrt{\frac{1}{N} \sum_{i = 0}^{N - 1} \left (v_{\mathrm{LPF}}[i] \right)^{2}} \tag{4}

です。ここでNは1周期分のN = 64としても良いのですが，ノイズ対策として格納しているすべてのデータ，つまり8周期分(N = 64 × 8 = 512)としました。

ところで，周波数が変動した場合，上記の「512サンプリング = 8周期分」という前提が崩れてしまうことになります。そう，8周期分のつもりが，7.9周期とか8.1周期になってしまうことがあり得ることになります。しかしこれは気にしないことにしました。周波数はほとんどの瞬間に49.9 ~ 50.1 Hzの間に収まるはずですので，312.5 μs × 512 = 160 msと8周期分の差は概ね320 μs，すなわち1サンプル分ほどに収まります。したがって，実効値の計算にはさほど影響しないと考えます。

スペクトルとTHDの計算

図5のHPFの後の信号vHPFを8周期分のリングバッファに格納します。ノイズ対策として8周期分（512ポイント）を平均した1周期分（64ポイント），かつ，虚部に零を詰め込んでcomplex<float>型の配列を作り，過去に実装したFFT関数にポインタとして渡しています。

negligible.hatenablog.com

64ポイントしかデータがありませんので，高調波の次数としては0次（直流），1次（基本波）から31次までしか得られません。なお，HPFを通しているため直流成分は零となります。2次以上の各次については1次（基本波）に対する比をパーセントで表して改めて別の配列に格納します。また，THD（total harmonic distortion，総合ひずみ率）を次式から計算します（ただし，31次高調波まで）。

\mathrm{THD} = \displaystyle \frac{\sqrt{\displaystyle \sum_{k = 2}^{31} V_{k}^{2}}}{V_{1}} \tag{5}

なお，周波数が変動した場合，「64サンプリング = 1周期分」という前提が崩れてしまうことになりますが，上記の実効値の場合と同様の理由から気にしないことにしました。周波数はほとんどの瞬間に49.9 ~ 50.1 Hzの間に収まるはずですので，312.5 μs × 64 = 20 msと1周期分の差は概ね40 μs = 0.128サンプル分に収まります。したがって，FFTの計算にはさほど影響しないと考えます。

ソフトウェアの実装

本章では，上記のサンプリングや信号処理をはじめとするソフトウェア実装に関して，その一部をソースコード例を示しながら概説します。すべてを網羅した説明は紙面（？）の都合上できませんので，GitHub上のソースコードをご参照ください。

github.com

ArduinoフレームワークとFreeRTOS

ArduinoフレームワークでESP32マイコンをプログラミングする場合，無意識にFreeRTOSというリアルタイムOSを利用することになります。Arduinoの象徴たるsetup()関数とloop()関数に加えて，FreeRTOSによって「タスク」を作る・使うことができます。これは，要するにマルチタスク，マルチプロセッシングです。ESP32マイコンはCPUコアを2つ持っているため，2つの処理を同時に走らせることができます（ペリフェラルなど共有部分の取り扱いがどうなっているのか，筆者も完全には理解していません…💧）。

図10: ソフトウェア実装構成

図10にソフトウェアの実装構成を示します。Arduinoではsetup()関数を抜けた後，loop()関数が繰り返し呼ばれる状態になります。まず，setup()関数でLCDの初期化，WiFiへの接続，リアルタイムクロックの設定（ntpを用いたntp.nict.jpとの同期）を実行します。

また，setup()関数の中で，電圧をサンプリングするタスクsampleVoltage()と，Google スプレッドシートにデータをPOSTするタスクpostGSSを作り，loop()関数が走るCore 1ではなくCore 0で走らせることによって処理を分散します。

loop()関数内ではサンプリングした電圧や零クロス情報から周波数や実効値を（calcFreq()関数やcalcRMS()関数を呼んで）計算し，結果をLCDに表示する他，リアルタイムクロックの時刻の「秒」を監視して，10秒毎にタスクであるpostGSS()関数に通知を送ります（Google スプレッドシートにデータをPOSTするため）。

いろいろ書いていますが，本プロジェクトではLang-Shipさんの記事を大いに参考にさせて頂いており，この場を借りて感謝申し上げます。

以下，いくつかのトピックについて，ソフト実装の方法を簡単に説明します。

lang-ship.com

定周期サンプリング

まずはAC 100 V系統の電圧を定周期でサンプリングしなければなりません。FreeRTOSで定周期サンプリング，すなわち定周期割り込みを作る場合，

  1. 定周期で動作させる処理を記述した関数を作り，while(true)による無限ループ内の先頭に通知を待つxTaskNotifyWait()関数を置く
  2. setup()関数内で上記の関数をタスクとする
  3. 同じくsetup()関数内でタイマ割込みを設定し割り込み関数を設定する
  4. 上記の割り込み関数を作り，xTaskNotifyFromISR()関数で上記1，2で作成したのタスク（関数）に通知を送る

このような流れになります。コード例を抜粋すれば下記となるでしょうか。

// 割り込み関数
void IRAM_ATTR onTimer() {
  BaseType_t taskWoken;
  // タスクに通知を送る
  xTaskNotifyFromISR(taskHandle_sample, 0, eIncrement, &taskWoken);
}

（中略）

// 定周期サンプリング関数
void sampleVoltage(void *pvParameters) {
  uint32_t ulNotifiedValue;
  while (true) {
    // 通知を待つ
    xTaskNotifyWait(0, 0, &ulNotifiedValue, portMAX_DELAY);
    （電圧のサンプリング，HPF，LPF，零クロス検知などの処理）
  }
}

（中略）

// setup関数
void setup() {
  （中略）

  // タスクの作成
  xTaskCreateUniversal(
    sampleVoltage,        // タスクとして設定する関数名
    "sampleVoltage",      // タスク名（任意）
    8192,                 // スタックメモリ
    NULL,                 // パラメータ
    24,                   // 優先度
    &taskHandle_sample,   // このタスクのタスクハンドラ
    PRO_CPU_NUM           // タスクを実行するCPUコア
  );

  // タイマ割込みの設定
  hw_timer_t * timer = NULL;
  timer = timerBegin(0, 40, true);
  timerAttachInterrupt(timer, &onTimer, true);
  timerAlarmWrite(timer, 625, true);
  timerAlarmEnable(timer);

  （中略）
}

26行目以降の「タスクの作成」の箇所で，優先度を最高の「24」に設定しています。また，このタスクを動作させるCPUコアとしてPRO_CPU_NUM (CPU0)を指定しました。これは，loop()関数がAPP_CPU_NUM (CPU1)で実行されるので，それとは別のコアに設定したいためです（ただし，WiFiを含む無線通信に関する機能がCPU0で動くため，注意が必要です）。優先度に関しては最大にする必要があるのか分かりませんでしたが，定周期性が崩れてしまうと周波数，実効値，スペクトルがいずれも正しく計算できないことから，ひとまず最大にしています（そしてそのまま…💧）。

37行目以降の「タイマ割込みの設定」にて割り込み発生の周期とISR (interrupt service routine)としてどの関数を呼ぶのかを設定しています。39行目でtimerBegin(0, 40, true)とした箇所にて，Timer0を使用し，80 MHzを40分周して0.5 μs毎にカウントアップするように設定します。次の行で，timerAttachInterrpt(timer, &onTimer, true)とすることで，冒頭のonTimer()関数をISRとして設定します。さらに次の行でtimerAlarmWrite(timer, 625, true)にて625カウント毎，つまり，0.5 μs × 625 = 312.5 μs毎に割り込みが発生するようにします。これはちょうど20 msの1 / 64となっていますので，50 Hzの1周期で64回サンプリングすることになります。最後にtimerAlarmEnable(timer)でタイマを有効にします。

すると，1行目以降のonTimer()関数が312.5 μs毎に呼ばれます。ここで，キーワード「IRAM_ATTR」はこの関数をESP32の内部RAMに置くようにというコンパイラへの指示になります。さて，このonTimer()関数内で電圧のサンプリングやHPF，LPFなどの信号処理を実行しても良さそうなのですが，割り込み関数(ISR)の実行時間が長くなるとクラッシュしてしまう（筆者も理由は完全に理解していません。何かのウォッチドッグタイマかな…？）ため，サンプリングと信号処理を他のタスク，この場合はsampleVoltage()関数に実装し，onTimer()関数からはこのタスクに通知を送るだけで，できるだけ早く抜けるようにします。一方，タスクの方は作られると常時動き続けるのでwhile(true)などで無限ループを作っておく必要があります。その先頭にxTaskNotifyWait()関数を置き，ISRであるonTimer()関数からの通知を待ちます。通知が来たところで電圧のサンプリングやHPF，LPF，零クロス検知などの信号処理を実行します。

ここでもLang-Shipさんのブログ記事を（大いに）参考にさせて頂いています。 lang-ship.com

ディジタルフィルタ

GPIO35に入力された電圧をサンプリングするのですが，オフセットやノイズを取り除くため，HPFとLPFを施します。かなり以前ですが，ディジタルLPFを実装した様子をtweetしました。

この時に作成したLPFは1次でしたが，今回は2次のバターワース特性(ζ = 0.707)のフィルタを使うことにしました。フィルタはC++らしくクラスとして実装します。例えばHPFは以下のような感じです。もっとC++らしく関数テンプレートでも使えばfloat型とdouble型の両方に適用できたのかもしれませんね。

// 2次HPF
class SecondOrderHPF {
  public:
    SecondOrderHPF(float Ts, float zeta, float omega_n);
    ~SecondOrderHPF(void);
    float apply(float u);

  private:
    float m, n;
    float a1, a2, b0, b1, b2;
    float y_old[2], u_old[2];
};

// コンストラクタ
SecondOrderHPF::SecondOrderHPF(float Ts, float zeta, float omega_n) {
  // 双一次変換に戻づく係数の設定
  m = 4.0 * zeta / (omega_n * Ts);
  n = 4.0 / pow(omega_n * Ts, 2.0);

  a1 = 2.0 * (1.0 - n) / (1.0 + m + n);
  a2 = (1.0 - m + n) / (1.0 + m + n);
  b0 = n / (1.0 + m + n);
  b1 = -2.0 * n / (1.0 + m + n);
  b2 = b0;

  // 初期値
  y_old[0] = y_old[1] = 0;
  u_old[0] = u_old[1] = 0;
}

// デストラクタ
SecondOrderHPF::~SecondOrderHPF(void) {}

// ディジタルフィルタの適用
float SecondOrderHPF::apply(float u) {
  float y = b0 * u + b1 * u_old[0] + b2 * u_old[1]
    - a1 * y_old[0] - a2 * y_old[1];
  y_old[1] = y_old[0];
  y_old[0] = y;
  u_old[1] = u_old[0];
  u_old[0] = u;
  return y;
}

ディジタルフィルタの係数については，本ブログの過去記事をそのまま適用しています。

negligible.hatenablog.com

これ，SymPyに頼って力業で求めたものですが，本当に計算しておいて良かったです…。いろんなところで役に立っています。

周波数の計算

アルゴリズムについては前述の通りです。ソフトウェアの実装としては以下のようになっています。最初の零クロスを線形補完で求め，2番目以降は零クロス間のサンプリング数を数え，最後の零クロスを同じく線形補完で求めます。途中の零クロスを線形補完で求めても意味がないので，サンプリング数を求めるだけになっています。周波数の計算に用いる零クロスの数N_ZEROCROSSは100になっていますので，求められる周波数は過去2秒間を窓とする移動平均になりますね。

なお，このcalcFreq()関数はloop()関数から呼ばれます。したがって，calcFreq()が周波数を計算しようと走っている間にタイマ割込みが発生し，sampleVoltage()関数によるタスクが動き出して100個の零クロス情報を格納しているオブジェクトの配列zcが変わってしまうことがあります。これをうまく解決しようと思ったのですが…とりあえず大きな誤差は生まなそうなので諦めることにしました。にはは…💧

// 零クロス情報から周波数を計算
float calcFreq() {
  float T = 0;

  // すべての零クロス情報から零クロス時間の合計を求める
  T += zc[0].V_pos / (zc[0].V_pos - zc[0].V_neg);
  for (int j = 1; j < N_ZEROCROSS; j++) {
    T += zc[j].N;
  }
  T -= zc[N_ZEROCROSS - 1].V_neg / (zc[N_ZEROCROSS - 1].V_pos 
    - zc[N_ZEROCROSS - 1].V_neg);
  T *= T_SAMPLE_US * 1e-6;

  // 周波数を計算してreturnする
  float f = 1.0 / T * (N_ZEROCROSS - 1);
  return f;
}

実効値の計算

同様に，アルゴリズムについては前述の通りです。ソフトウェアの実装としては以下のようになっています。

// 実効値の計算
float calcRMS() {
  float V_rms = 0;
  for (int j = 0; j < N_SAMPLE; j++) {
    V_rms += pow(v_lpf[j], 2);
  }
  V_rms /= N_SAMPLE;
  V_rms = sqrt(V_rms);
  return V_rms;
}

うむ。straight-forward万歳です。ここではサンプリングしてHPF，LPFを施した電圧情報が格納されているfloat型の配列v_lpfをグローバル変数にしていますが，要素数と一緒にポインタ渡しすれば，もっと汎用的な関数として使えると思います。

LCDの描画に関するソフトウェア実装

せっかくカラーLCDを使うので，眺めていて楽しい画面を作りたいですよね。ただ数値を表示するだけではやっぱり物足りませんので…。そこで，らびやんさんのLovyanGFXを活用して，周波数と実効値のアナログメータ，波形，スペクトルを表示することにしました。

画面デザインと完成イメージ

本プロジェクトを始めるにあたり，最初に考えたのが画面デザインです。やっぱり見た目から入らないとモチベーションが湧きませんよね。図11にその際に構想した画面イメージを示します。

図11: 画面イメージ

今回は表示したい情報が多いので扇型のアナログメータはスペース的に難しくなります。そこで，直線型（と言うのかな）のアナログメータを2つ用いて周波数と実効値を視覚化することにします。さらに，1周期分の波形とスペクトルを表示します。これがリアルタイムにうねうね動いたらなかなか楽しいゾ──と想像しながらかいていました。なおこの画像はLCDの解像度に合わせて240ドット × 240ドットとしてMicrosoftペイントで描きました。

この完成イメージが冒頭の動画となります。

アナログメータ

今回の直線型アナログメータもC++らしくクラス(hMeter)として作成しました。パラメータを変えることで，周波数用と実効値用の両方に使っています。なお，アナログメータの左隣の数値はloop()関数内でベタに書いています。sprintf()関数の作りなどに柔軟性が欲しかったため，このクラスには入れませんでした。

class hMeter {
  public:
    hMeter(int _x0, int _y0, float _u_min, float _u_max);
    ~hMeter();
    void drawFrame();
    void drawLabels(); 
    void update(float u);

  private:
    const int width = 120;
    const int height = 24;
    const int colBG = TFT_BLACK;
    const int colFrame = TFT_DARKGREY;
    const int colLabel = TFT_WHITE;
    const int colHand = TFT_RED;
    int x0, y0, x_old;
    float u_min, u_max;
};

extern LGFX lcd;

// コンストラクタ
hMeter::hMeter(int _x0, int _y0, float _u_min, float _u_max) {
  // Initializing instant variables
  x0 = _x0;
  y0 = _y0;
  x_old = 0;
  u_min = _u_min;
  u_max = _u_max;
}

// デストラクタ
hMeter::~hMeter() {}

// フレームを描く
void hMeter::drawFrame() {
  
  // 水平線
  lcd.drawFastHLine(x0, y0, width, colFrame);

  // 大目盛り
  lcd.drawFastVLine(x0, y0 - height / 2, height + 1, colFrame);
  lcd.drawFastVLine(x0 + width / 2, y0 - height / 2, height + 1, colFrame);
  lcd.drawFastVLine(x0 + width, y0 - height / 2, height + 1, colFrame);

  // 中目盛り
  lcd.drawFastVLine(x0 + width / 4, y0 - height / 4, height / 2 + 1, colFrame);
  lcd.drawFastVLine(x0 + 3 * width / 4, y0 - height / 4, height / 2 + 1, colFrame);

  // 小目盛り
  for (int i = 0; i < 20; i++) {
    lcd.drawFastVLine(x0 + i * width / 20, y0 - height / 8,
    height / 4 + 1, colFrame);
  }  
}

// ラベルを描き込む
void hMeter::drawLabels() {  
  lcd.setFont(&fonts::Font2);
  lcd.setTextColor(colLabel);
  lcd.setTextDatum(top_center);
  lcd.drawNumber(static_cast<int>(u_min), x0, y0 + height / 2 + 2);
  lcd.drawNumber(static_cast<int>((u_max + u_min) / 2),
    x0 + width / 2, y0 + height / 2 + 2);
  lcd.drawNumber(static_cast<int>(u_max),
    x0 + width, y0 + height / 2 + 2);
}

// 表示をアップデートする
void hMeter::update(float u) {
  lcd.drawFastVLine(x0 + x_old, y0 - height / 2, height, colBG);
  drawFrame();
  int x = constrain(static_cast<int>((u - u_min) / (u_max - u_min) * width),
    0, width);
  lcd.drawFastVLine(x0 + x, y0 - height / 2, height, colHand);
  x_old = x;
}
波形

波形についてもクラス(tdPlot)として作りました。サンプリングが1周期あたり64点なのですが，波形表示部の横幅を128ドットとし，線形補完して（単にdrawLineして）描いています。右端のフレームの色が変わってしまっているように見えるので，まだバグが残っているのかもしれません。実用上問題ないため放置しています…。

class tdPlot {
  public:
    tdPlot(int _x0, int _y0, float _u_min, float _u_max);
    ~tdPlot();
    void drawFrame();
    void drawLabels();
    void update(float *u);

  private:
    const int N = 64;
    const int width = 128;
    const int height = 45;
    const int colBG = TFT_BLACK;
    const int colFrame = TFT_DARKGREY;
    const int colLabel = TFT_WHITE;
    const int colPlot = TFT_YELLOW;
    int x0, y0;
    float u_min, u_max;
    float u_old;
};

extern LGFX lcd;

// コンストラクタ
tdPlot::tdPlot(int _x0, int _y0, float _u_min, float _u_max) {
  // Initializing instant variables
  x0 = _x0;
  y0 = _y0;
  u_min = _u_min;
  u_max = _u_max;
}

// デストラクタ
tdPlot::~tdPlot() {}

// フレームを描く
void tdPlot::drawFrame() {
  lcd.drawRect(x0 - 1, y0 - 1, width + 2, height + 2, colFrame);
  lcd.drawFastHLine(x0, y0 + height / 2 + 1, width, colFrame);
}

// ラベルを描く
void tdPlot::drawLabels() {
  lcd.setFont(&fonts::Font2);
  lcd.setTextColor(colLabel);
  lcd.setTextDatum(middle_right);
  lcd.drawNumber(static_cast<int>(u_max), x0 - 2, y0);
  lcd.drawNumber(static_cast<int>((u_max + u_min) / 2),
    x0 - 2, y0 + height / 2);
  lcd.drawNumber(static_cast<int>(u_min), x0 - 2, y0 + height);
}

// 波形を更新する
void tdPlot::update(float *u) {
  // 古い波形を削除して水平線を書き直す
  lcd.fillRect(x0, y0, width, height / 2 + 1, colBG);
  lcd.fillRect(x0, y0 + height / 2 + 2, width, height / 2 - 1, colBG);
  lcd.drawFastHLine(x0, y0 + height / 2 + 1, width, colFrame);

  // 最初の点を描く
  int y = constrain(static_cast<int>((u[0] - u_min) / (u_max - u_min) * height),
    1, height - 1);
  lcd.drawPixel(x0, y0 + height - y, colPlot);
  int y_old = y;

  // 残りの点を描く
  for (int i = 1; i < N; i++) {
    int y = constrain(static_cast<int>((u[i] - u_min) / (u_max - u_min) * height),
      1, height - 1);
    lcd.drawLine(x0 + (i - 1) * 2, y0 + height - y_old, x0 + i * 2,
      y0 + height - y, colPlot);
    y_old = y;
  }

  // 最後の点を描く
  y = constrain(static_cast<int>((u[0] - u_min) / (u_max - u_min) * height),
    1, height - 1);
  lcd.drawLine(x0 + width - 2, y0 + height - y_old, x0 + width,
    y0 + height - y, colPlot);
}
スペクトル

スペクトルに関してもほぼ同様ですが，棒グラフを地道に描くことになります。

class spectrumPlot {
  public:
    spectrumPlot(int _x0, int _y0, float _u_max);
    ~spectrumPlot();
    void drawFrame();
    void drawLabels();
    void update(float *u);

  private:
    const int N = 64;
    const int width = 128;
    const int height = 45;
    const int colBG = TFT_BLACK;
    const int colFrame = TFT_DARKGREY;
    const int colLabel = TFT_WHITE;
    const int colPlot = TFT_SKYBLUE;
    int x0, y0;
    float u_min, u_max;
    float u_old;
};

extern LGFX lcd;

// コンストラクタ
spectrumPlot::spectrumPlot(int _x0, int _y0, float _u_max) {
  // Initializing instant variables
  x0 = _x0;
  y0 = _y0;
  u_max = _u_max;
}

// デストラクタ
spectrumPlot::~spectrumPlot() {}

// フレームを描く
void spectrumPlot::drawFrame() {
  lcd.drawRect(x0 - 1, y0 - 1, width + 2, height + 2, colFrame);
  lcd.drawFastHLine(x0, y0 + height / 2 + 1, width, colFrame);
}

// ラベルを描く
void spectrumPlot::drawLabels() {
  // 縦軸
  lcd.setFont(&fonts::Font2);
  lcd.setTextColor(colLabel);
  lcd.setTextDatum(middle_right);
  lcd.drawNumber(static_cast<int>(u_max), x0 - 2, y0);
  lcd.drawNumber(static_cast<int>(u_max / 2),
    x0 - 2, y0 + height / 2);
  lcd.drawNumber(0, x0 - 2, y0 + height);

  // 横軸
  lcd.setFont(&fonts::Font0);
  lcd.setTextColor(colLabel);
  lcd.setTextDatum(top_center);
  for (int i = 1; i < N / 2; i += 4) {
    lcd.drawNumber(i, x0 + i * 4 + 1, y0 + height + 2);
  }
}

// スペクトルを更新する
void spectrumPlot::update(float *u) {
  // 古いスペクトルを削除して水平線を書き直す
  lcd.fillRect(x0, y0, width, height / 2 + 1, colBG);
  lcd.fillRect(x0, y0 + height / 2 + 2, width, height / 2 - 1, colBG);
  lcd.drawFastHLine(x0, y0 + height / 2 + 1, width, colFrame);

  // スペクトルをプロットする
  for (int i = 1; i < N / 2; i++) {
    int y = constrain(static_cast<int>(u[i] / u_max * height), 0, height);
    lcd.drawFastVLine(x0 + i * 4, y0 + height - y, y, colPlot);
    lcd.drawFastVLine(x0 + i * 4 + 1, y0 + height - y, y, colPlot);
  }
}

Google スプレッドシートとApp Script

これまでの筆者のIoTプロジェクトと同様に，観測したデータをGoogle スプレッドシートにPOSTします。方法については本ブログの過去の記事で取り上げたものと同様です*13。

観測例

これまでにいくつか特異な観測例を得ていますのでご紹介します。

2022年6月の電力需給逼迫時の周波数低め運用

2022年は早くも6月27日に関東甲信で梅雨明けが発表され，直後に猛暑に襲われました。このため，冷房需要の急騰から電力需給が逼迫し，日中に周波数が公称値の50 Hzを離れ，概ね49.9 Hzを中心に運用されることがありました。図12に6月29日の観測例（Google スプレッドシートに作ったダッシュボード）を示します。気温が上昇してきた午前10時頃に運用の中心が49.9 Hzに変わった様子が見られます。

図12: 2022年6月29日の観測例（ダッシュボード）

図12 ~ 14に2022年6月27日，29日，30日の電力系統周波数の観測結果（それぞれ1日分）を示します。昨今の世界情勢と6月に猛暑を迎えたことから大変厳しい需給状況であったことが窺えますね。なお，橙色のプロットは5分間の移動平均です。

図13: 2022年6月27日の電力系統周波数の観測結果

図14: 2022年6月29日の電力系統周波数の観測結果

図15: 2022年6月30日の電力系統周波数の観測結果

落雷の影響かもしれない電圧変動

周波数の観測だけであれば，本ブログの過去の記事*14から進歩はありません。やはりと言うか，電圧（実効値）を観測することで新たな発見がありました。

図16: 落雷の影響かもしれない電圧変動の観測例

図16は，落雷の影響であるかもしれない電圧変動の観測例です。これは2022年6月3日の夕方ですが，このとき，ちょうど積乱雲が通過しており，暴風雨と雷が襲い掛かっておりました。15時38分頃にステップ状の電圧上昇，15時42分頃に同じくステップ状の電圧低下が観測されています。ちょうどこの時間帯に我が家に給電しているかもしれない配電用変電所に落雷があったことが「Yahoo! Japan 天気・災害」の「雨雲レーダー」から分かっており，落雷に伴う電圧変動と送り出し電圧の再調整があった可能性を示唆しています。とは言え，同じ時間帯に周波数変動も起きており，より大きな範囲での負荷変動があったのかもしれません。

ただし電圧はローカルなパラメータなので，電力系統全体の変動なのか，近傍に限られた現象なのかは峻別が困難です。次で述べるように，宅内での家電機器の入切も大きな影響を及ぼします。また，同じ100 V低圧系統に繋がっているご近所さんの影響も受けるでしょう。

電圧実効値の許容範囲

電気事業法施行規則第38条には電圧の許容範囲が記載されており，100 V系統の場合には「101ボルトの上下6ボルトを超えない値」，つまり101 ± 6 V (95 V ~ 107 V)が許容される変動の範囲となります。

観測を続けていると，確かにこの範囲に収めようとしている様子が見て取れますが（変圧器の負荷時タップ切替器を使って配電用変電所の送り出し電圧を調整しているように見える），時々，106 Vギリギリ，いや，107 Vを超えてしまうような場合を観測しています。図17は106 Vまで上がった場合の例です。このとき，自作装置の表示を疑ってHIOKIのテスタでも同様に測定しましたが，やはり同じでした。

図17: 電圧実効値が106 Vに達する

平日で多くの事業所で就業開始となる午前9時前後，お昼休みの始まる12時前後，休憩があると思われる15時前後，終業時刻となる17時前後にやはり電圧変動が観測され易い傾向にありますね。もちろん，系統全体の負荷が変わるので周波数も影響を受けて同じ時刻に変動しますが，これまでのところ電圧の方が変動を受けやすいように見えます。いずれにしてもとても面白い観測ができています。

宅内でのkW級の負荷投入

これは本プロジェクトがまだブレッドボード上のプロトタイピングの段階でしたが，4月のまだ寒い日に1,200 Wのオイルヒータの投入に伴う電圧実効値の変動をはっきり捉えることができました。動画2をご覧下さい。

youtu.be
動画2: 1,200 Wのオイルヒータの入切に伴う電圧変動を捉える

このヒータやヘアドライヤーなど大きな負荷を投入すると照明が暗くなったりするので分かってはいましたが，これはまさに「おおぅ！？ こんなに変わるの？」と率直な感想が漏れました。

高調波の観測例

スペクトルを画面に表示していますが，ずっと観測していると，3次高調波が変動しつつも常に概ね2%存在しており，5次高調波が昼と夜で大きく変動することが判りました。また，需給が逼迫するようなときには，なぜか5次高調波よりも7次高調波が大きく現れてくるという傾向があることも判りました。とても不思議です。

3次および9次，15次など3n次高調波は三相系統においては零相成分であるため，変圧器を超えて影響することはなく，100 V系統だけで決まっているのではないかと思われますが（ただし三相平衡を仮定した場合），5次，7次などはローカルな低圧系統とグローバルな高圧，特別高圧，超高圧系統のどの範囲から来たものなのか，判然としません。

自宅とは違う場所ということで，某月某日都内某所に自作のこの監視装置を持ち込んで，1時間ほど高調波を測ってみました。

図18: 某月某日都内某所での高調波の観測結果

7次，11次高調波が5次，9次高調波よりも大きいという，自宅とは全く異なる傾向が見られました。自宅で11次が見られないのは測定系の作りが間違っているからかも──と疑っていましたが，この結果から，存在するなら見えることが確認できました。

まとめ

以上，本記事ではESP32マイコンを用いた商用電源品質監視装置のハードウェア構成，サンプリングと信号処理，ソフトウェア構成について述べ，また，観測例についていくつか紹介しました。

この観測装置は筆者の机の上で常時稼働しており，周波数，電圧実効値，波形，スペクトルが常にうねうねと動き続けている様子がいつでも見られます。まさに，「電力系統の息吹を感じる」ことができます。筆者のような特異体質を持っている場合には癒し効果まで期待できそうです。

今後も観測を続け，夏・冬の需給逼迫時，地震発生に伴う突発的な周波数や電圧の変動などについて観測できた場合は，本記事への追記という形などでレポート致します。

*1:筆者よりも先行してのるさんが同様の装置を発表されています。M5Stackだけで商用電源周波数の変動計測 - | ProtoPedia

*2:FFTを訪ねて［後編］PythonとC++で作ってみる - The Negligible Lab

*3:ただしソフトの作りについては，本ブログ以上に詳しく記述されている方がいらっしゃり，例えば筆者はLang-shipさんの記事を大いに参考にさせて頂いております。ESP32の高精度タイマー割り込みを調べる | Lang-ship

*4:7，8年前にやはり電源品質監視装置を作ろうとして集めていた部品です。

*5:例えばこの変圧器の6 V側にダイオード整流器と3端子レギュレータを取り付けてESP32マイコンの電源とすることもできると思いますが，ダイオード整流器の高調波電流が変圧器の漏れインダクタンスに電圧降下を生じると，高調波スペクトルの計測に誤りが生じる恐れがあります。

*6:今回はGPIO35を使います。

*7:ブレッドボードでのプロトタイピング時に電源モジュールからと思われるノイズがそこそこ見えていたため。

*8:当初，これとは異なるより小さな変圧器を使う計画でしたが，テストしてみると無負荷でも2次側の波形が歪んでいたため，今の変圧器を使うことにしました。

*9:本当は128回，サンプリング周期156.25 μsにしたかったのですが，何故か正常に動作しませんでした。

*10:分圧器に使用した抵抗器の誤差などに起因する。

*11:実効値についてはLPFの前の信号を使うべきだったかもしれません

*12:メソッドを作っていないので単に構造体でも可。

*13:M5Stackで遊ぶ ～浮遊静電容量で系統周波数を測る～ - The Negligible Lab

*14:M5Stackで遊ぶ ～浮遊静電容量で系統周波数を測る～ - The Negligible Lab