背景：部屋の空気質(CO2、温度、湿度）をクラウドに上げて集計しているが、気圧が計測できていなかった。高精度の気圧センサ―らしいQMP6988を使って気圧を計測、集計する
課題：QMP6988はレジスタの値を読んだ後、補正が必要なのだが補正式が結構複雑。Arduino/MicroPyhon版のドライバであればGitHubに公開さえているがRust版は無いようだ（もし公開されていてもいざ使おうとすると引数の型変換等エラー出まくって苦労しそう）
取り組み：プログラミングの勉強も兼ねてRust版QMP6988のドライバを作ってみる
詳細：
補正式はQMP6988の仕様書にかかれているのだけど*1、結構分かりくい。最初読んだときは全く分からなかった。今も完全には理解していないけど、多分こういうことだろうと解釈して自分がPGしやすくなるよう適当にレイアウトを変えた。
補正式を見ながらゴリゴリ計算したら値が得られるのだろう。。気になるのは、GitHubに掲載されているドライバは補正で使っている係数の値が小数点になっていなかった。テーブルに書かれているE-NNって、10のマイナスNN乗と思うのだが、、誤差を抑えるため演算の最後にまとめて10のNN乗で割っているのか？
Rustで試行錯誤せず、まずはPythonで演算式を作った方がよさそうだ。(GitHubにもPython版ドライバがあるのだけど、struct moduleにより実装されており、ちょっと意味分からず）
Pythonでサンプル実装してみた。
#!/usr/bin/python3
# -----------------------------
# fc: 22 TEMP_TXD0
# fb: 2b TEMP_TXD1
# fa: 69 TEMP_TXD2
# -----------------------------
# f9: be PRESS_TXD0
# f8: 7e PRESS_TXD1
# f7: 9f PRESS_TXD2
# -----------------------------
#
#
# b8: 81 COE_b00_a0_ex
# ------------------------
# b7: 73 COE_a2_0
# b6: ec COE_a2_1
# b5: c4 COE_a1_0
# b4: e4 COE_a1_1
# b3: c6 COE_a0_0
# b2: f7 COE_a0_1
# b1: d4 COE_bp3_0
# b0: f6 COE_bp3_1
# af: 09 COE_b21_0
# ae: ee COE_b21_1
# ad: ec COE_b12_0
# ac: 1f COE_b12_1
# ab: ac COE_bp2_0
# aa: 04 COE_bp2_1
# a9: 9f COE_b11_0
# a8: 0d COE_b11_1
# a7: c6 COE_bp1_0
# a6: 12 COE_bp1_1
# a5: cd COE_bt2_0
# a4: 05 COE_bt2_1
# a3: 16 COE_bt1_0
# a2: 13 COE_bt1_1
# a1: c1 COE_b00_0
# a0: 5b COE_b00_1
#
#
# A S
# a1 -6.30E-03 4.30E-04
# a2 -1.90E-11 1.20E-10
# bt1 1.00E-01 9.10E-02
# bt2 1.20E-08 1.20E-06
# bp1 3.30E-02 1.90E-02
# b11 2.10E-07 1.40E-07
# bp2 -6.30E-10 3.50E-10
# b12 2.90E-13 7.60E-13
# b21 2.10E-15 1.20E-14
# bp3 1.30E-16 7.90E-17
#
conversion_factor = {
'a1': (-6.30E-03, 4.30E-04),
'a2': (-1.90E-11, 1.20E-10),
'bt1': ( 1.00E-01, 9.10E-02),
'bt2': ( 1.20E-08, 1.20E-06),
'bp1': ( 3.30E-02, 1.90E-02),
'b11': ( 2.10E-07, 1.40E-07),
'bp2': ( -6.30E-10, 3.50E-10),
'b12': ( 2.90E-13, 7.60E-13),
'b21': ( 2.10E-15, 1.20E-14),
'bp3': ( 1.30E-16, 7.90E-17),
}
cf = conversion_factor
TEMP_TXD0 = 0x22
TEMP_TXD1 = 0x2b
TEMP_TXD2 = 0x69
PRESS_TXD0 = 0xbe
PRESS_TXD1 = 0x7e
PRESS_TXD2 = 0x9f
COE_b00_a0_ex = 0x81
C0E_a2_0 = 0x73
C0E_a2_1 = 0xec
C0E_a1_0 = 0xc4
C0E_a1_1 = 0xe4
C0E_a0_0 = 0xc6
C0E_a0_1 = 0xf7
C0E_bp3_0 = 0xd4
C0E_bp3_1 = 0xf6
C0E_b21_0 = 0x09
C0E_b21_1 = 0xee
C0E_b12_0 = 0xec
C0E_b12_1 = 0x1f
C0E_bp2_0 = 0xac
C0E_bp2_1 = 0x04
C0E_b11_0 = 0x9f
C0E_b11_1 = 0x0d
C0E_bp1_0 = 0xc6
C0E_bp1_1 = 0x12
C0E_bt2_0 = 0xcd
C0E_bt2_1 = 0x05
C0E_bt1_0 = 0x16
C0E_bt1_1 = 0x13
C0E_b00_0 = 0xc1
C0E_b00_1 = 0x5b
C0E_b00_ex = COE_b00_a0_ex >> 4
C0E_a0_ex = COE_b00_a0_ex & 0x0F
a1 = cf['a1'][0] + cf[ 'a1'][1]*((C0E_a1_1<<8) + C0E_a1_0) / 32767
a2 = cf['a2'][0] + cf[ 'a2'][1]*((C0E_a2_1<<8) + C0E_a2_0) / 32767
bt1 = cf['bt1'][0] + cf['bt1'][1]*((C0E_bt1_1<<8) + C0E_bt1_0) / 32767
bt2 = cf['bt2'][0] + cf['bt2'][1]*((C0E_bt2_1<<8) + C0E_bt2_0) / 32767
bp1 = cf['bp1'][0] + cf['bp1'][1]*((C0E_bp1_1<<8) + C0E_bp1_0) / 32767
b11 = cf['b11'][0] + cf['b11'][1]*((C0E_b11_1<<8) + C0E_b11_0) / 32767
bp2 = cf['bp2'][0] + cf['bp2'][1]*((C0E_bp2_1<<8) + C0E_bp2_0) / 32767
b12 = cf['b12'][0] + cf['b12'][1]*((C0E_b12_1<<8) + C0E_b12_0) / 32767
b21 = cf['b21'][0] + cf['b21'][1]*((C0E_b21_1<<8) + C0E_b21_0) / 32767
bp3 = cf['bp3'][0] + cf['bp3'][1]*((C0E_bp3_1<<8) + C0E_bp3_0) / 32767
a0 = ((C0E_a0_1 << 12) + (C0E_a0_0 << 4) + C0E_a0_ex) /16
b00 = ((C0E_b00_1 << 12) + (C0E_b00_0 << 4) + C0E_b00_ex) /16
Dt = (TEMP_TXD2 << 16) + (TEMP_TXD1 << 8) + TEMP_TXD0 - 2**23
Dp = (PRESS_TXD2 << 16) + ( PRESS_TXD1 << 8) + PRESS_TXD0 - 2**23
print(a1, a2, bt1, bt2, bp1, b11, bp2, b12, b21, bp3, a0, b00)
print(Dt, Dp)
Tr = a0 + a1 * Dt + a2 * Dt * Dt
Pr = b00 + bt1 * Tr + bp1*Dp + b11*Tr*Dp + bt2*Tr*Tr + bp2*Dp*Dp + b12 * Dp*Tr*Tr + b21 * Dp * Dp * Tr + bp3 * Dp * Dp * Dp
print("Tr:", Tr)
print("Pr:", Pr)出力結果は以下の通り、何か間違っている。やはり先人のコードを見て検証しないとだめだ。
Tr: 72160.48901926124 Pr: 151907.33116591972
Rustで試行錯誤するのでは時間がかかりすぎるので、、比較的分かりやすそうなAliOS-Things様のドライバでPMP6988を制御してみた。結果以下と表示された。温度は大体合ってそうなのだが、気圧が99926ってのはどんな単位なんだろうか。
pressure:99926.062500 temp:33.898438
CALC_INTフラグがあって、これは多分整数演算を優先するモードと思うのだが、これを０にすると結果が少し変わって以下
pressure:100009.906250 temp:33.976463
1hPa = 100Pa らしく、上記の単位がPaだとすると、今の気圧は1000hPa(今日は晴天なので、、そんなものかも)。気象情報のサイトで確認すると、1009hPaとのことで、ちょっと誤差がある。めちゃくちゃ精度が高いと売りなのに、、いいのか？？ 標高が関係している？？*2 細かい精度はおいといて、だいたい合っている。次はこのPythonコードをばらして、どんな計算をしているのか、その時の係数は何か？を確認すると、自分のコードのどこが間違っているのかわかるはず。。
自分の実装ミス
仕様書から読み解けなかったのだが、係数を生成する時、単にビットを連結していただけだったのだが、この係数は符号付きの値なのであった。だから、符号付き20bitと解釈しないといけない。(単純にunsigned の20bitと考えてのだがそれは間違い）
サンプルソースで、12bit左シフトして、signed型にしてから右12bitシフトしてるのだが、なぜわざわざシフトして戻しているのか最初分からなかった（ソースは下記）。これは係数を12bit左シフトすることで、20bitの最上位を32bit長の最上位bit(MSB)までずらして符号の位置に持って行ってからsigned 32bit長に変換して、符号付き整数に型変換した後、右シフトして元の値に戻しているのであった。
qmp6988.qmp6988_cali.COE_a0 =
(QMP6988_S32_t)(((a_data_uint8_tr[18] << SHIFT_LEFT_12_POSITION) |
(a_data_uint8_tr[19] << SHIFT_LEFT_4_POSITION) |
(a_data_uint8_tr[24] & 0x0f))
<< 12);
qmp6988.qmp6988_cali.COE_a0 = qmp6988.qmp6988_cali.COE_a0 >> 12;MicroPyton版ドライバでも同様にunsigned->signedにするためシフトして戻している
COE_a0 = self.int32(((a_data_u8r[18] << SHIFT_LEFT_12_POSITION) | (a_data_u8r[19] << SHIFT_LEFT_4_POSITION) | (a_data_u8r[24] & 0x0f)) << 12)
COE_a0 = self.int32(COE_a0 >> SHIFT_LEFT_12_POSITION)
def int32(self, dat):
#return int(dat)
if dat > (1 << 31):
return dat - (1 << 32)
else:
return dat係数は8bit + 8bitを結合して16bit長にする演算もあり、ソースを確認するとそちらも同様に符号付き16bitと解釈する必要があるようであった。その辺りを修正したら、サンプル実装と同じ値になるのではと期待。
■追記
16bit長の補正用係数、20bit長の補正用係数いずれも符号付き数値と解釈してサンプルコードを作り直した。コードは以下
#!/usr/bin/python3
conversion_factor = {
'a1': (-6.30E-03, 4.30E-04),
'a2': (-1.90E-11, 1.20E-10),
'bt1': ( 1.00E-01, 9.10E-02),
'bt2': ( 1.20E-08, 1.20E-06),
'bp1': ( 3.30E-02, 1.90E-02),
'b11': ( 2.10E-07, 1.40E-07),
'bp2': ( -6.30E-10, 3.50E-10),
'b12': ( 2.90E-13, 7.60E-13),
'b21': ( 2.10E-15, 1.20E-14),
'bp3': ( 1.30E-16, 7.90E-17),
}
cf = conversion_factor
PRESS_TXD2 = 159
PRESS_TXD1 = 114
PRESS_TXD0 = 230
TEMP_TXD2 = 105
TEMP_TXD1 = 53
TEMP_TXD0 = 224
COE_b00_a0_ex = 0x81
COE_a2_0 = 0x73
COE_a2_1 = 0xec
COE_a1_0 = 0xc4
COE_a1_1 = 0xe4
COE_a0_0 = 0xc6
COE_a0_1 = 0xf7
COE_bp3_0 = 0xd4
COE_bp3_1 = 0xf6
COE_b21_0 = 0x09
COE_b21_1 = 0xee
COE_b12_0 = 0xec
COE_b12_1 = 0x1f
COE_bp2_0 = 0xac
COE_bp2_1 = 0x04
COE_b11_0 = 0x9f
COE_b11_1 = 0x0d
COE_bp1_0 = 0xc6
COE_bp1_1 = 0x12
COE_bt2_0 = 0xcd
COE_bt2_1 = 0x05
COE_bt1_0 = 0x16
COE_bt1_1 = 0x13
COE_b00_0 = 0xc1
COE_b00_1 = 0x5b
COE_b00_ex = COE_b00_a0_ex >> 4
COE_a0_ex = COE_b00_a0_ex & 0x0F
import pdb
def cnv2s16(val):
#print(hex(val))
if val & 0x80_00 == 0:
pass
else:
val = val - 0x1_00_00
#pdb.set_trace()
return val
def cnv2s20(val):
#print(hex(val))
if val & 0x8_00_00 == 0:
pass
else:
val = val - 0x10_00_00
#pdb.set_trace()
return val
a1 = cf['a1'][0] + cf[ 'a1'][1]*cnv2s16((COE_a1_1<<8) + COE_a1_0) / 32767
a2 = cf['a2'][0] + cf[ 'a2'][1]*cnv2s16((COE_a2_1<<8) + COE_a2_0) / 32767
bt1 = cf['bt1'][0] + cf['bt1'][1]*cnv2s16((COE_bt1_1<<8) + COE_bt1_0) / 32767
bt2 = cf['bt2'][0] + cf['bt2'][1]*cnv2s16((COE_bt2_1<<8) + COE_bt2_0) / 32767
bp1 = cf['bp1'][0] + cf['bp1'][1]*cnv2s16((COE_bp1_1<<8) + COE_bp1_0) / 32767
b11 = cf['b11'][0] + cf['b11'][1]*cnv2s16((COE_b11_1<<8) + COE_b11_0) / 32767
bp2 = cf['bp2'][0] + cf['bp2'][1]*cnv2s16((COE_bp2_1<<8) + COE_bp2_0) / 32767
b12 = cf['b12'][0] + cf['b12'][1]*cnv2s16((COE_b12_1<<8) + COE_b12_0) / 32767
b21 = cf['b21'][0] + cf['b21'][1]*cnv2s16((COE_b21_1<<8) + COE_b21_0) / 32767
bp3 = cf['bp3'][0] + cf['bp3'][1]*cnv2s16((COE_bp3_1<<8) + COE_bp3_0) / 32767
a0 = cnv2s20((COE_a0_1 << 12) + (COE_a0_0 << 4) + COE_a0_ex) / 16
b00 = cnv2s20((COE_b00_1 << 12) + (COE_b00_0 << 4) + COE_b00_ex) / 16
Dt = (TEMP_TXD2 << 16) + (TEMP_TXD1 << 8) + TEMP_TXD0 - 2**23
Dp = (PRESS_TXD2 << 16) + ( PRESS_TXD1 << 8) + PRESS_TXD0 - 2**23
print("------------------------")
print("Dt", Dt)
print("Dp", Dp)
print("------------------------")
Tr = a0 + a1 * Dt + a2 * Dt * Dt
Pr = b00 + bt1 * Tr + bp1*Dp + b11*Tr*Dp + bt2*Tr*Tr + bp2*Dp*Dp + b12 * Dp*Tr*Tr + b21 * Dp * Dp * Tr + bp3 * Dp * Dp * Dp
print("Tr:", Tr)
print("Tr:", Tr/256.0)
print("Pr:", Pr)
print("Pr:", Pr/100.0, "hPa")
実行結果は以下となり、計測タイミングの違い、丸め誤差？で結果に差があるが、サンプルコードによる出力とかなり近い値になっている。
------------------------ Dt -1493536 Dp 2061030 ------------------------ Tr: 7356.718938025345 Tr: 28.737183351661503 Pr: 100029.38430679489 Pr: 1000.293843067949 hPa
Pythonによる計算式はできたと判断して、これをRust版に置き換えると・・
■追記
かつて発生して悩まされた？ wait usb downloadのメッセージで止まってしまう件が再発した、これはbootモードがおかしくなっているためと推測。GPIO9がbootモードの選択を担当しており、誤ってLになっていたりするとReset投入後Downloadモードかなにかになってしまう（たしか。仕様書で再確認必要だが）。このような状況になったら、GPIO9を確実にHに固定するため、もしオープンだったらPullUp等を介して3.3Vに引き上げると通常通りにBootする可能性が高い。自分の場合は、GPIO9をH固定にしたら、起動の問題は解消した。一度正常に起動すると、GPIO9のH固定をはずしてオープンにしても普通に起動するのであった。ただ、、危ういのは危ういので、PullUpしておくのが無難かも
~/lang/rust/esp32-c3/qmp6988$ cargo espflash monitor
[2024-07-07T14:44:28Z INFO ] 🚀 A new version of cargo-espflash is available: v3.1.0
[2024-07-07T14:44:28Z INFO ] Serial port: '/dev/ttyACM0'
[2024-07-07T14:44:28Z INFO ] Connecting...
[2024-07-07T14:44:29Z INFO ] Using flash stub
Commands:
CTRL+R Reset chip
CTRL+C Exit
ESP-ROM:esp32c3-api1-20210207
Build:Feb 7 2021
rst:0x15 (USB_UART_CHIP_RESET),boot:0x0 (USB_BOOT)
Saved PC:0x4038063e
wait usb download■追記
Rustによるドライバは以下となった（係数計算手抜き版）
短い理由：補正用の係数計算はQMP6988のレジスタから値を取ってきて計算式を適用すべきなのだが、レジスタはいつも同じ値のようなので、あらかじめPython等で計算した定数で代用する*3。
pub struct QMP6988 {
}
impl QMP6988 {
const a0 :f32 = -2105.9375;
const a1 :f32 = -0.006391493270668661;
const a2 :f32 = -3.732941679128391e-11;
const b00:f32 = 23489.5;
const bt1:f32 = 0.11356932279427473;
const bt2:f32 = 6.638398388622699e-08;
const bp1:f32 = 0.03578676717429121;
const b11:f32 = 2.248985259559923e-07;
const bp2:f32 = -6.172249519333475e-10;
const b12:f32 = 4.795419171727653e-13;
const b21:f32 = 4.15744499038667e-16;
const bp3:f32 = 1.2433906064027832e-16;
pub fn calc_temp_press(PRESS_TXD2:u8, PRESS_TXD1:u8, PRESS_TXD0:u8, TEMP_TXD2:u8, TEMP_TXD1:u8, TEMP_TXD0:u8) -> (f32, f32){
let Dt = (((TEMP_TXD2 as u32) << 16) + ((TEMP_TXD1 as u32) << 8) + TEMP_TXD0 as u32) as f32 - 2_u64.pow(23) as f32;
let Dp = (((PRESS_TXD2 as u32) << 16) + (( PRESS_TXD1 as u32) << 8) + PRESS_TXD0 as u32) as f32 - 2_u64.pow(23) as f32;
let Tr = QMP6988::a0 + QMP6988::a1 * Dt + QMP6988::a2 * Dt * Dt;
let Pr = QMP6988::b00 + QMP6988::bt1 * Tr + QMP6988::bp1 * Dp + QMP6988::b11 * Tr * Dp + QMP6988::bt2 * Tr * Tr + QMP6988::bp2 * Dp * Dp + QMP6988::b12 * Dp * Tr * Tr + QMP6988::b21 * Dp * Dp * Tr + QMP6988::bp3 * Dp * Dp * Dp;
(Tr/256.0, Pr/100.0)
}
}インスタンス変数？使っていないのでStruct型にしなくてもいい感じですが、、ゆくゆくは係数も正しく計算したいので、一旦Struct型で実装。理想的にはSoftI2CドライバのインスタンスはQMP6988に持たせたい。
■参考URL
原因特定につながると期待できるリファレンス実装(Arduino)
M5Unit-ENV/src/QMP6988.cpp at master · m5stack/M5Unit-ENV · GitHub
Python版ドライバ（比較的分かりやすい実装）
AliOS-Things/haas_lib_bundles/python/libraries/qmp6988/qmp6988.py at master · alibaba/AliOS-Things · GitHub
仕様書
QMP6988 Datasheet PDF - Digital Barometric Pressure Sensor
