CTセンサのスパイクに疑似的なHampelフィルタで対応してみた

入居しているシェアハウスの家電製品の動作をネットワーク上で確認できるように、AC100Vの電流をCTセンサで計測して、ＡＤＳ１０１５使用　ＰＧＡ機能搭載１２ｂｉｔＡＤコンバーターADS1015でAD変換して、ESP32でデータを送信するというシステムを作ってみたところ、信号にスパイクが度々発生したので、ソフトウェア的にフィルタを作成して対処しました。備忘録として残します。

システム構成

AC100Vを検出するCTセンサ
CTセンサから生じる電流を電圧に変換する100Ω抵抗
CTセンサの電圧をAD変換するADコンバータ
いつ買ったのか覚えていないESP32開発ボード

これらを接続してAC100Vの電流を計測します。

ADコンバータの接続方法ですが、Adafruitのページが参考になります。

なぜ参考になるのかと言えば、秋月電子通商のADS1015を使ったADコンバータは、回路がAdafruitのADS1015搭載 12BitADC 4CH 可変ゲインアンプ付きとほぼ同じだからです。

ADコンバータは、GNDを0Vとして1端子の電圧を計測する絶対電圧計測と、2端子間の電圧を計測する差動電圧計測の2パターンがありますが、絶対電圧計測はGND経由でノイズが入って面倒だったりするので、今回は差動電圧計測 を行います。

コード

同様にESP32のコードやライブラリもArduinoのものを流用します。

https://learn.adafruit.com/adafruit-4-channel-adc-breakouts/arduino-code
まずはArduinoIDEのシリアルプロッタでAC100Vの電流を計測してみます。

#include <Wire.h>
#include <Adafruit_ADS1015.h>
 
Adafruit_ADS1015 ads1015(0x48);

void setup(void)
{
  Serial.begin(115200);
  ads1015.setGain(GAIN_ONE);
  ads1015.begin();
}
 
void loop(void)
{
  int16_t results;
 
  results = ads1015.readADC_Differential_0_1();

  Serial.println(results); 

}

100V1Wの電球を接続して計測してみると以下の波形が得られます

数値としては-1～1でADコンバーターの性能としては限界ギリギリになりますが、100V1Wつまり0.01Aの有無を検出できていることがわかります。今回は電力計測が目的ではないので、これで十分でしょう。
電力を計測したい場合は基準電源を用いてダミーロードで実測値を出し、補正を加えてやる必要があると思いますが、ONOFFが分かればいい今回はそこまで求めません。

原因不明のスパイク

しかし計測を続けていくうちに問題が発生しました。

謎のスパイクが計測されています。数値としては下の画像のようになります。

このスパイクの出所がまったくわからず、対策できませんでした。
ESP32内臓のADCでスパイクが出るのなら、これは回路側の問題だと検討がつきますが、今回は秋月のADC基板に社外のI2cのライブラリを組み込んで使っているので、問題の究明が困難です。I2cではデジタル通信で送信しているので、その辺りでのバグも考えられます。

やはり同じADS1015だからといってAdafruitのライブラリを秋月のADCで使うのは無理があったのでしょうか、

しかし、このスパイク、かなり特徴的なのでソフトウェア的に除去できそうです。

スパイクの考察

まずこのスパイクがどのようなものかを考察してみます。

周期的ではない
大きさはADCで検出できる値の1000倍以上
データ数は僅か
スパイクの大きさは一定ではない

要するに突然巨大でランダムな外れ値が出るというのがこのスパイクのようです。

このスパイクを除去するのに、ノイズを平滑化するために多用される移動平均フィルタ等の平均値を用いたフィルタは使えません。1000倍ものスパイクが入っているので、1000回サンプリングして平均をとって周囲に散らしたところで値が1を超えたものになり、ADCの信号が潰れてしまいます。

このスパイクは明らかにエラーと考えて、エラーをいかに判別するか、そして妥当なデータを入力するかを考えたほうが良さそうです。

そこでこのスパイクを、外れ値除去に特化したHampelフィルタで除去してみたいと思います。

Hampelフィルタ

Hampelフィルタとは標準偏差を利用したフィルタで

標準偏差がσの母集団では、正常な値は中央値から+3σ,-3σの範囲に収まる

ということを利用したフィルタです。画像は平均0、標準偏差1の正規分布のグラフです。この仕組みを利用してフィルタを作ります。

具体的には、一定回数のサンプリングを行い、中央値、標準偏差値σを算出し、中央値から3σ異常離れた値は異常値として除去し中央値に置き換えます。単純な一定値以上、又は以下のカットと違って、統計的に外れ値を除去できるフィルタですので、今後大電流の機器を使用しても同じようにスパイクの除去が期待できます。

今回重要なのは、AC100Vの電流を計測しているという点です。
AC100Vの電圧と電流はsin波で表されますので、実用上は中央値も平均値も0にしておいて問題がありません。

つまりsin波なので±の振れ幅が同じなので中央値は0、そして平均は0

\mu =\dfrac{1}{n}\sum ^{n}_{i=1}x_{i}=0

標準偏差の2乗である分散は次にような式、つまり2乗和の平均になります

\sigma ^{2}=\dfrac{1}{n}\sum ^{n}_{i=1}\left( x_{i}-\mu \right) ^{2}=\dfrac{1}{n}\sum ^{n}_{i=1}x_{i}^{2}

このあたりを踏まえてコードを書くと次のようになります。

#include <Wire.h>
#include <Adafruit_ADS1015.h>

#define N 1000　//サンプリング数
 
Adafruit_ADS1015 ads1015(0x48);

long data_buff[N];//フィルタ用データバッファ

void setup(void)
{
  Serial.begin(115200);
  ads1015.setGain(GAIN_ONE);
  ads1015.begin();
}
 
void loop(void)
{
  int16_t results;
  
  long sum_x2 = 0; //sum(x^2)

  for(int i = 0; i< N;i++){
    results = ads1015.readADC_Differential_0_1();
    data_buff[i] = results;
    //分散を計算 平均も中央値も0だからこれで可能
    sum_x2 += (long)results * (long)results;
  }
  //ハードウェア的に平均は0なのでこういった実装が可能
  double sigma = sqrt(sum_x2 / N);//標準偏差
  //σが小さくなることが度々あったので、最低値保障　これより小さいと1を検出できないことがある
  //スパイクの排除にはこれで十分
  if(sigma < 1)sigma = 1.0;
  
 for(int i = 0; i< N;i++){
    if(data_buff[i] > 3.0 * sigma||data_buff[i] <-3.0 * sigma){
      data_buff[i] = 0;
      }
 }

 for(int i = 0; i< N;i++){
  Serial.println(data_buff[i]);
 }
}

このコードは、

1000回サンプリングを行いデータバッファに保存する
分散を計算する
分散から標準偏差を求める
標準偏差が小さすぎる場合、標準偏差を1にする
データバッファの値を、中央値0から+3σ以上、-3σ以下の場合は中央値である0にする
データバッファの値を出力する

という仕組みになっています。サンプリングを完了してから出力するという処理の都合上、リアルタイムに出力できるわけではありませんが、

今回特殊な処理があるとすれば

標準偏差が小さすぎる場合、標準偏差を1にする

の部分です。

標準偏差σの3倍以上の場合は中央値にする処理をしているので、標準偏差σが極小の場合、例えばσ=0.1の場合3σでも0.3

ADコンバーターはint出力なので絶対値が1以下の値を出力しません から、3σが1よりも小さくなってしまうと、中央値である0に置き換えられ、データの1の値が消えることによりデータが失われてしまいます。

そこでσが1より小さくなった場合は1に置き換えるという処理をすることで、データの喪失を防いでいます。

これで電流を計測すると次にようになります。

ほぼ波形は変わっておらず、スパイクは除去されました。

まとめ

想定外のスパイクでしたが、hampelフィルタを通すことで、リアルタイム性が失われる代わりにスパイクを除去することができました。

ADコンバーターの仕様による特殊処理も含めていい学習になったと思います。