MATLAB体温計モドキ

**本記事で紹介するのは体温計「モドキ」です。ご自身の体調管理には正規の体温計を使用してください**
# 1 概要
App Designerの練習のため、予測式体温計モドキを作ったのでレポートします。プログラムはすべて[こちら](https://github.com/HppyCtrlEngnrng/MATLAB_Thermometer_Modoki)にアップロードしています。
## 予測式体温計の概要
一般家庭に流通している体温計は、温度センサを腋に挟んで温度を計測しています。腋をぴったりと閉じることにより、温度センサの示度は測定したい体内の温度(深部体温)に近づいていき、やがて定常値に収束します。この定常値を体温(腋窩温)とするわけですが、腋窩温の収束には10分程度の時間がかかります(図1)。
![thermal_dynamics.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211913/07cb8910-7c85-bde0-e61f-45eb0f43f7b9.png)
図1 筆者の温度特性

この計測時間を短縮するため、事前設計された温度変化モデルに基づき短時間の計測データから定常値を予測するのが予測式体温計です。

## 作成した環境の機能仕様
今回は上記の予測式体温計の原理を模擬した簡易的な体温計モドキをArduino+MATLABで作成していきます。実現したい機能としては以下の通りとしました。

1. サーミスタを用いた温度センサユニットを有し、腋に挟んだサーミスタの検出温度(腋窩温度)を連続計測できること
2. 1で計測された腋窩温度データと、温度変化モデルにより平衡温度を予測できること
3. 2の**予測アルゴリズムと温度変化モデルはユーザが変更できること**

# 2 構成
図2に作成した環境の機能構成を示します。本環境は、Arduino Nanoベースのサーミスタ温度計測ユニットとデータ処理用のMATLABがインストールされたPCで構成されています。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211913/b1161eaa-ab70-97e6-bcf8-46c49aedd038.png)
図2 構成図

## 2.1 温度計測ユニットと機能
温度計測ユニットの外観図および回路図を図3に示します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211913/e1559242-ed7b-3197-6cfd-ade3c40025c5.png)
図3 温度計測ユニットのハードウェア
### 2.1.1 温度→電圧変換
サーミスタは抵抗変化型のセンサで、今回は単純な分圧回路によりこの抵抗変化を検出します。サーミスタの温度$\theta_T\ [\mathrm{K}]$と抵抗値$R_{T}\ [\Omega]$の関係は、ある基準温度$\theta_0\ [\mathrm{K}]$およびそのときの抵抗値$R_0\ [\Omega]$、そして定数$B$を用いて

$$
\theta_{T}=\frac{1}{\frac{1}{B}\ln\left(\frac{R_{T}}{R_0}\right)+\frac{1}{\theta_0}}
$$
と表されます。分圧則より、$R_{th}$はArduinoからの供給電圧$V_i\ [\mathrm{V}]$と基準抵抗$R_{f}\ [\Omega]$、検出電圧$V_o\ [\mathrm{V}]$により

$$
R_{th}=\frac{V_i-V_o}{V_o}R_f
$$
と求まるため、$V_o$の計測により$\theta_T$が計算されます。

### 2.1.2 計測値エンコード
Arduino NanoのADCの分解能は10ビットです。通常であればこれを2バイトに分割しシリアル通信すればよいのですが、生データに観測雑音の影響が大きかったことから、4 msで取得したADC出力を25回加算してから100 ms周期で送信することにします。これにより送信データの値域は0～1023×25(0b0000 0000 0000 0000～0b0110 0011 1110 0111)となります。これを上下7ビットに分割し、末尾に0b1111 1111を付与した3バイトを送信します(下位側を左シフトしてますがたぶん不要)。これにより、MATLAB側で終端子を255として簡単に受信処理を記述できます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211913/e872c105-5a9a-2933-030b-6470b78891e4.png)

図4 送信データのエンコーディング

### 2.1.3 シリアル送信
Serial.writeで送るだけです。エンコード処理を含めた関数は以下のようになります。

``` C:[Arduino]エンコード処理
byte snd_buf[3];

void send_mes() {
  snd_buf[0] = (byte)((mes_val & 0x7F) << 1);
  snd_buf[1] = (byte)(mes_val >> 7);
  // snd_buf[2] = 0xFF;

  Serial.write(snd_buf,3);
}
```

## 2.2 MATLAB側の機能
### 2.2.1 シリアル入力
serialportオブジェクトを作成し、温度計測ユニットから送られたデータを受信します。configureTerminator関数で終端子を255に指定することで、readline関数による受信待機を行います。

``` MATLAB:[MATLAB]シリアル通信設定
s = serialport(comport_prev, app.serial_baud);
configureTerminator(s, 255);
```

### 2.2.2 計測値デコード
readline関数で読みだされた受信データを温度データに変換します。まずサーミスタ回路の電圧値を復元するために以下の処理を行います。
1. readlineの戻り値はstringであるため、char array→double arrayに変換しArduinoの送信データを復元する
2. 2バイトに分割された送信データを結合し、計測値の25回加算値を得る
3. ADCの分解能を補正して25で割り電圧値を得る
一行にまとめると以下の処理となります。

``` MATLAB:[MATLAB]デコード処理
double(char(readline(s))) * [2^-1; 2^7] * 5/1023/25
```
つぎに電圧値を温度へ変換するために、2.1.1の式をもとに以下の関数を使用します。

``` MATLAB:[MATLAB]電圧→温度変換
function Th = conv_thermister(app, Vo, B, R0, Rf, T0, Vi)
    Rth = Rf*(Vi-Vo)./Vo;
    Th = 1./(1/B*log(Rth/R0)+1/T0)  - 273.15 ;
end
```

### 2.2.3 温度変化モデル同定
2.2.2の手順を用いて温度変化を記録(最大10分)、ユーザ指定のアルゴリズム、温度変化モデルでパラメータ同定を行う機能です。App Designerの画面構成は図5の通りです。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211913/1af6ee17-4347-c0d6-0e09-df721546e110.png)
図5 モデル同定画面

Arduinoが接続されたCOMポートを指定し、サーミスタを腋に挟んでから[計測開始]ボタンを押下すると、100 ms周期の温度が記録されます(図6)。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211913/69668444-557e-c281-ddcc-9b0ae2e83ada.png)
図6 モデル同定用温度計測

計測されたデータは[データ保存]でmatファイルにエクスポートでき、また過去に計測したデータを[データ読み込み]でロードできます。同定を実行するには、まずデータが計測されている状態で[アルゴリズムを選択]します。アルゴリズムはclassとして定義され、App側とのI/Fとして計測データからユーザのモデルパラメータと計測データへのフィッティングカーブを出力するidentメソッドを持ちます。例えば、温度変化モデルを深部体温(定数)→腋窩の表面温度→サーミスタ検出温度の熱伝導と考え、各々時定数の異なる線形1次モデルとした図7の場合、identメソッドは計測データからパラメータ$\theta_B,a_A,a_T$を推定します。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211913/c253c6ec-d561-ec43-64db-21d35f6ef895.png)
図7 温度変化モデルの例

[同定実行]を押下すると指定されたアルゴリズム(クラス)のidentメソッドがコールされ、計測データにフィッティングカーブが重畳されます(図8)。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/211913/30414687-c738-f565-0db6-2d827a785ff6.png)
図8 同定結果の表示

### 2.2.4 体温測定
2.2.3で得られたモデルを用い、リアルタイム計測された体温データから平衡点を予測出力する機能です。体温測定を実行するには、サーミスタを腋に挟んでから[体温予測開始]をクリックします。同定の場合と同様、アルゴリズムのクラスに定義されたestimメソッドが1サンプルごとにコールされ、予測値を出力します。予測アルゴリズムは、例えば図7のモデルに対して深部体温をステップ状の外乱とみなした外乱オブザーバを用いることができます。動作例は以下のツイートをご覧ください。

<blockquote class="twitter-tweet"><p lang="ja" dir="ltr">App Designerのれんしう<br><br>サーミスタがあったのでMATLABで予測式体温計モドキを作ってみた。 <a href="https://t.co/8OlHoRRUCi">pic.twitter.com/8OlHoRRUCi</a></p>&mdash; 幸福の制御工学 (@HppyCtrlEngnrng) <a href="https://twitter.com/HppyCtrlEngnrng/status/1249392523920830464?ref_src=twsrc%5Etfw">April 12, 2020</a></blockquote> <script async src="https://platform.twitter.com/widgets.js" charset="utf-8"></script>