RPi ZeroW で高精度AD-HATを使う ~地震計の製作~

# 背景
去年のアドベントカレンダーで紹介した地震計は1Hz前後の低周波ノイズを除去しきれなかった。

https://qiita.com/pirotan628/items/01bf7c0ab0b28334aaed

今年は同じテーマで別のAD変換ボードを使ってみることにした。蓑虫クリップ以外のハンダ付けが無かったので、自作と言えるかは微妙。

# ハードウェア
使用したのは千石で手に入る「[High-Precision AD HAT](https://www.waveshare.com/wiki/High-Precision_AD_HAT)」。量子化ビット数が32bitと高分解能で、買ってすぐに使えるのが決め手。

マイコンはAD変換ボードがぴったりサイズの Raspberry Pi Zero W を選択（当初 Raspberry Pi 4 でテストしたが、ノイズが大きく断念した。手待ちの電源に制約があり、本体と電源どちらの問題か切り分けができなかった）。

取得したデータをmseed形式で保存したいと思い、計測から1次処理まで1台で済ませるにも丁度良い組み合わせとなった。

GPIOにHATを接続するだけでハードウェアの準備は完了。ジオフォンをIN1,IN2番端子に繋いで観測する。今回の用途では差動入力で微小電圧の測定を行う。これらはソフトウェア側で設定できる。


<img width=300 src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/88076/fd55fff0-e8b8-0981-36cc-24d0a666cdc8.jpeg">

<img width=300 src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/88076/8027a9e4-a994-a966-7d4d-be11e3cbd29f.jpeg">


# ソフトウェア
OSは Raspberry Pi OS Lite (32bit; Kernel Ver. 6.1)。Raspberry Pi はオーバーレイファイルシステムで動かし、データはNASをマウントして保存することにした。

ADS1263の制御には発売元のWikiに従い[githubから入手できるサンプルプログラム](https://github.com/waveshare/High-Pricision_AD_HAT)の完成度が高く、ほとんどそのまま使える。

AD値の読み出しからmseed形式での保存までをPythonで用意した。サンプリングは100Hzで行い、タイムスタンプを付けるためにインターバルタイムハンドラを使用した。

使用した環境では "gpiomem-bcm2835" が存在せず、ADS1263.py でimportされる公式の config.py ではRaspberryPiZeroがJestonと認識されてしまう。そこで、124行目を下記のように変更した。

```config.py
# 変更前：
if os.path.exists('/sys/bus/platform/drivers/gpiomem-bcm2835'):

# 変更後：
if os.path.exists('/sys/bus/platform/drivers/gpiomem-bcm2835') or os.path.exists('/sys/bus/platform/drivers/rpi-gpiomem'):
```


<details><summary>ソースコード全体はこちら</summary><div>

```seismo2.py
#!/usr/bin/python
# -*- coding:utf-8 -*-

import time
import ADS1263
import RPi.GPIO as GPIO

import signal
import sys
import os

from datetime import datetime, timedelta, timezone

from obspy import Stream, Trace
from obspy.core import UTCDateTime
import numpy as np


REF = 5.08          # Modify according to actual voltage
                    # external AVDD and AVSS(Default), or internal 2.5V
ADC = ADS1263.ADS1263()

JST = timezone(timedelta(hours=+9), 'JST')

dataFileLoc = "/mnt/nas/obs_geophone/raw" #location of the mseed file

stnName = "HOM"

statsZ = {}
statsZ['network'] = "HOME" #fake network name
statsZ['station'] = stnName
statsZ['component'] = 'Z'

## Station location
statsZ['stla'] = 34
statsZ['stlo'] = 135

## Event location
statsZ['evla'] = 34
statsZ['evlo'] = 135
statsZ['evdp'] = 0

statsZ['sampling_rate'] = 100.0

dt_now = datetime.now(JST)
fnb = dt_now.strftime('%Y%m%d%H%M%S')
fn = fnb + ".csv"

rawdata = 0
ns = 0
tr = []
fnb="0000000000"

# mseed file length by min
#minutes = 1
minutes = 30
#minutes = 60
#minutes = 180

sdlen = minutes * 60 * statsZ['sampling_rate']

channelList = [0]  # The channel must be less than 10

def writedata():
    global ns
    global statsZ
    global fnb
    global tr

    fnb = dt_now.strftime('%Y%m%d%H%M%S')
    statsZ['starttime'] = UTCDateTime(dt_now)
    tra = np.array(tr)
    tr = []
    print(fnb,tra,ns)

    ns = 0    
    trZ = Trace(data=tra, header = statsZ)
    st = Stream(traces=[trZ])
    st.write(os.path.join(dataFileLoc, f"{fnb}-Z.mseed"), format="MSEED")

    return 0

# インターバルタイマーハンドラー
def intervalHandler(signum,frame):
    global ns
    global tr
    global rawdata
    global dt_now

    ns = ns + 1
    if ns == 1: dt_now = datetime.now(JST)
    val = rawdata
    tr.append(val)

    if ns == sdlen:
        dummy = writedata()

def getVol(raw):
    if(raw>>31 ==1):
        data = -1*(REF*2 - raw * REF / 0x80000000)
    else:
        data = (raw * REF / 0x7fffffff)
   
    return int(data*100000000)


def get30Int(raw):
    if(raw>>31 ==1):
        data = -1 * (0xffffffff - raw)
    else:
        data = raw
   
    return int(data/4)  # STEIM2 can treat only 30bit...


def readData():
    global rawdata
    
    # The faster the rate, the worse the stability
    # and the need to choose a suitable digital filter(REG_MODE1)

#    if (ADC.ADS1263_init_ADC1('ADS1263_400SPS') == -1):
    if (ADC.ADS1263_init_ADC1('ADS1263_100SPS') == -1):
        exit()

    ADC.ADS1263_SetMode(1) # 0 is singleChannel, 1 is diffChannel
   
    while(1):
        ADC_Value = ADC.ADS1263_GetAll(channelList)    # get ADC1 value

#        rawdata = getVol(ADC_Value[0])
        rawdata = get30Int(ADC_Value[0])

    ADC.ADS1263_Exit()

    return rawdata

#=====================================================================

# インターバル起床するハンドラを設定(SIGALRMハンドラー)
signal.signal(signal.SIGALRM,intervalHandler)

# 秒の変わり目をポーリング
while True:
    dt_now = datetime.now()
    if dt_now.microsecond < 1000: # 秒の変わり目から1ms未満ならOK
        print('microsec = %d' % dt_now.microsecond)
        fnb = dt_now.strftime('%Y%m%d%H%M%S')
        break

# インターバルタイマー設定
# 最初の開始タイミング1秒後、以後0.01秒間隔
signal.setitimer(signal.ITIMER_REAL, 1, 0.01)


try:
    while True:
        rawdata = readData()

except IOError as e:
    print(e)
   
except KeyboardInterrupt:
    print("ctrl + c:")
    print("Program end")
    ADC.ADS1263_Exit()
    dummy = writedata()
    sys.exit()
```
</div></details>

# 結果
得られた波形、PSDとスペクトログラムはこちら。

地動ノイズが少ない時間帯のデータ：
<img width=500 src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/88076/2f80d48f-faee-58d8-c962-2a3d50cc9ad1.png">

地動ノイズが多い時間帯のデータ：
<img width=500 src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/88076/757f28ca-6cc1-903a-2405-42dcc92c78e0.png">

40Hz付近のピークが非常に目立つ。fs=100Hzなので、fs/2で折り返した60Hzのエイリアスと思われる。それ以外の低周波ノイズはSX8725Cを用いたときに比べて目立たないが、10秒程度の間隔でステップノイズが混入している。

24時間分のドラム波形では、相変わらず交通由来の地動ノイズが良く見える。
<img width=600 src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/88076/f12224e9-a27a-490b-0e46-5363d2074ad6.png">

## 電源の比較
ここでスイッチングノイズによる影響を検討するため、電源を変えてみた。上のデータは写真右のELECOMを使用していた。比較として iPad mini 5 に付属の電源を使用した。

<img width=400 src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/88076/e3123bbe-e1bf-2d6d-b03f-d1f13e620dbc.jpeg">

### ELECOM
交通量が多い時間帯の背景ノイズは約 +/-1e5 ～ +/-1.5e5[counts]
<img width=500 src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/88076/bb1358f1-8023-e3d5-19e8-e1b973702540.png">

### Apple
交通量が多い時間帯の背景ノイズは約 +/-2e5[counts]
<img width=500 src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/88076/208c95fe-88c4-52d5-6cff-7c5972c93397.png">

パワースペクトルでみると10Hz以上は大差ないのですが、10Hz以下の比較ではAppleの充電器が1桁程度高いノイズレベルを示しています。高精度AD変換をRaspberryPiで行う場合、ラズパイオーディオと同様に電源ノイズを気にしないといけない事が判明した。

# まとめ
- High-Precision AD HAT を RPi ZeroW で制御して、ジオフォンからの微小な振動を記録する事ができた。
- ラズパイでアナログ信号を扱うには、電源ノイズを克服する事が重要である。可聴域およびさらに低周波側のノイズ対策のため、オーディオ関連のノウハウを勉強したい。

去年作成したデバイスより良い物を目指して新しいデバイスを作成してみたものの、結果はまたもやノイズとの戦いになった。前回と比べると地動ノイズに対して機器側のセルフノイズが小さいので、観測装置としては改善されたと思う。

今回の場合、電源さえ低ノイズにすればもう少し良い結果が期待できるという事で、バッテリを用いた試験もやってみたい（モバイルバッテリもノイズが出ていそうだけど、むしろ高周波ノイズは気にしないので、もしかすると有りかもしれない）。

2023/12/22追記:
実際にモバイルバッテリーでやってみた。

https://qiita.com/pirotan628/items/dc1815c7bac5872b54fb