第35回信号処理シンポジウム
2020年11月4日〜6日（リモート開催）

付け爪型センサを用いたケプストラム分析及び MUSIC 法に基づく心拍推定Heart Rate Estimation Based on Cepstrum Analysis and MUSIC Method Using
Nail Tip Sensor
梶谷奈未 †

北村大地 †† 香川高等専門学校

石井耕平 †

Nami KAJITANI†
Daichi KITAMURA†
Kohei ISHI††National Institute of Technology, Kagawa College
アブストラクト

不快な装着感が無く長期間にわたる強固な固定が可能なウェアラブルセンサとして，付け爪型センサが提案されている．このセンサでは，指先の血流の変化によって光電脈波信号が得られるため，実生活環境下での常時心拍モニタリングへの応用が期待されている．
本研究では，実生活環境下において付け爪型センサから得られた光電脈波信号に対して，音響信号処理分野で一般的なケプストラム分析及び multiple signal classiﬁcation（MUSIC）法に基づく心拍推定を行い，その推定精度について検討する．

Fig 1: Nail tip sensor system [1].
1

はじめに在宅医療では，患者の健康状態を管理するために，遠隔でモニタリングを行うことが重要である．この目的のために，患者は心拍や血圧等の情報を取得するセンサを装着し続ける必要があり，その煩雑さから生活の質の低下を招く問題がある．この問題の解決として，文献 [1] では，Fig. 1 に示すような，不快な装着感のない付け爪型の遠隔モニタリング用センサが開発されている．試作された付け爪型センサは，Fig. 2 のような小型脈拍計測回路を搭載しており，爪に固定した状態で長時間の光電脈波信号（photoplethysmogram: PPG）を測定することがで

Fig 2: Tiny PPG sensor circuit and its application [1].

きる．この PPG を解析することで心拍が推定でき，冒頭で述べた在宅医療における遠隔モニタリングへの応用がが，その推定結果は外乱の影響を強く受けてしまい，推期待されている．

定精度の劣化が生じている．また理想的には，デッドタイ

付け爪型センサで得られる PPG は，指先の血管の膨ムの設定値は心拍の真値に合わせて変動させるべきパラ張や収縮に起因する反射光の変動を測定したものである． メータであるが，その具体的な方法は提案されておらず，この信号の基本周波数は通常心拍数と一致する．従って， 常に一定の値を用いているという問題もある．
測定された PPG の基本周波数を推定することで心拍の推本稿では，付け爪型センサで得られる PPG に対して，定値が得られるが，生活環境下では，指先の運動や圧力より外乱に頑健かつ高精度な心拍の推定方法を検討する．の変化等の外乱が混在するため，頑健な心拍の推定は困この推定問題は，音響信号処理で一般的な基本周波数推定難な問題となる．文献 [1] では，PPG に対してデッドタと等価であることから，音響信号処理で広く用いられていイムを設けたゼロクロス法を適用し心拍を解析しているるケプストラム分析 [2] と多重信号分類（MUltiple SIgnal‑ 5 ‑Classiﬁcation: MUSIC）法 [3] を PPG に適用する．
Voltage [a.u.]

Observed PPG signal

2.1

特徴

文献 [1] で提案されている付け爪型センサの利点は以下の 3 点である．

• 爪には感覚神経が無いため，不快な装着感がない

10-1

Observed ECG signal

2
Voltage [a.u.]

付け爪型センサの概要

20-2
18:00

• 爪は組織が固いため，数週間固定を維持できる

22:00

02:00
06:00
10:00
Japanese standard time

• 複数の爪に装着することで多点同時計測ができ，よりPPG (04:00)

Voltage [a.u.]

0.3

動作原理

センサシステムの概要と外観図を Figs. 1 及び 2 に示す．

0Voltage [a.u.]

線送信機，及び制御用マイコンが搭載されている．脈波

0.3受信機へ送られる．さらに，インターネット回線を通じPPG (14:15)

200-1-0.3

計測回路にて観測された PPG は，無線通信により外部のPPG (07:18)

0.3

1

付け爪型センサ本体には脈波計測回路，バッテリー，無-0.3

ECG (04:00)-2ECG (07:18)

0.3

ECG (14:15)

21000-1-0.3て遠隔地にあるモニタリング端末へと情報を送ることが0123Time [s]45-0.30123Time [s]45-20123Time [s]45(b) ﬁve seconds in sleep (left), moderate activity (center), andhard exercise (right) periods可能である．

2.318:00(a) 24 hours

発展的な計測が期待できる

2.2

14:00

測定波形

Fig 3: Observed PPG (top) and ECG (bottom) signals

文献 [1] で実施された 24 時間の測定実験で，付け爪型

for (a) 24 hours and (b) ﬁve seconds.

センサから得られた PPG 及び心臓付近に取り付けた電極から得られた心電信号（electrocardiograph: ECG）を次に，(x[t])Tt=1 を短時間区間信号に分割する．短時間

Fig. 3 に示す．まず，Fig. 3 (a) で示す 24 時間全体の波形区間の信号長を L，短時間区間のシフト長を S とすると，を見ると，大部分において波形がクリッピングしているこ(x[t])Tt=1 の短時間区間信号への分割は次式で表される．とが分かる．これは，日常的な活動の動作時に生じる（心拍以外の）信号（本稿ではこれをノイズ信号と考える）が

xs [i, l] = x[l + (i − 1)S]加算されることによって，センサの測定可能範囲を超えてここで，i = 1, 2, · · · , I は短時間区間インデクス，l =しまうことが原因である．しかしながら，Fig. 3 (b) 左列

1, 2, · · · , L は短時間区間内の離散時間インデクスを表す．に示す睡眠時の拡大波形では，PPG 及び ECG のいずれの波形も心拍に由来すると思われる波形（本稿ではこれ3.2(2)ケプストラム分析に基づく心拍推定を目的信号と考える）が測定されている．ただし，ECGケプストラム分析 [2] では，まず (xs [i, l])L
l=1 に何らかのはスパイク状の目的信号が現れるのに対し，付け爪型セ窓関数 (w[l])Ll=1 を乗じて L 点離散 Fourier 変換（discreteンサで得られる PPG の目的信号は高周波成分をあまり含Fourier transform: DFT）を施し，短時間区間の複素ス

まない波形となっている．Fig. 3 (b) 中央列の通常活動下ペクトル (F [i, ω])Ω
ω=1 を得る．

や，Fig. 3 (c) 右列の激しい運動中では，運動の激しさにF [i, ω] =応じてノイズ信号が増し，心拍に由来する信号の SN 比がL∑l=1低下する．

ここで，j =
基本周波数推定法を用いた心拍推定

3
3.1

(3)

−1 であり，ω = 1, 2, · · · , Ω は離散周波数インデクスである（但し Ω = L）．実数ケプストラムΩ(C[i, q])Q
q=1 は (F [i, ω])ω=1 の絶対値の対数を逆 DFT し信号の前処理観測した PPG を (xraw [t])Tt=1 とする．ここで，t =て得られる．

1, 2, · · · , T は離散時間インデクスであり，T は観測時間C[i, q] =長である．この信号を零平均の信号 (x[t])Tt=1 に変換する．

1∑
xraw [t]T t=1
Tx[t] = xraw [t] −√][2π(ω−1)(l−1)
w[l]xs [i, l] exp −jL(1)[]Ω1 ∑2π(ω−1)(l−1)
log |F [i, ω]| exp j(4)Ω ω=1L

ここで，q = 1, 2, · · · , Q はケフレンシーと呼ばれるケプストラム領域のインデクスである（但し Q = L）．即ち，ケ‑ 6 ‑1
プストラム (C[i, q])Q
q=1 は対数スペクトルのスペクトル列 Ri が次式のように得られる．

であり，基本周波数 fo [i] とその倍音 (nfo [i])N
n=2 の周波Ri = E[xi,m xTi,m ]数をもつ正弦波から成る目的信号のスペクトル解析がで

1 ∑K×Kxi,m xTi,m ∈ R
M m=1M≈

きるため，基本周波数推定に用いられる．なお，ケプストラム (C[i, q])Q
q=1 の時間変動はケプストログラムと呼ば(8)

ここで，·T は転置を表す．

れる．
短時間区間毎の基本周波数の推定値 (fˆo [i])Ii=1 は，次式のようにケプストラム (C[i, q])Q
q=1 が最大値を取るケフレI
ンシーインデクス (q̂o [i])i=1 から求められる．

fˆo [i] = fs /(q̂o [i] − 1)q̂o [i] = arg max C[i, q]

MUSIC 法では，信号相関行列 Ri の部分空間から目的信号成分とノイズ成分を分離し，MUSIC スペクトルと呼ばれる特徴量を算出する．信号相関行列 Ri を固有値分解すると，次式となる．

(5)

Ri = Vi Λi Vi−1
λi,10λi,2
 0
Λi = ..
 ..
 .
.

(6)q

ここで，fs は信号のサンプリング周波数である．実際は，基本周波数の推定範囲の下限 fmin と上限 fmax を事前に定めたうえで，それらに対応するケフレンシーインデク

Vi =0···..
.

0..
.

0···

0

λi,Kvi,1vi,2···vi,K[スの下限 qmin = round(fs /fmax + 1) 及び上限 qmax =···





(9)(10)](11)round(fs /fmin + 1) の範囲内で (6) の最大値探索を行う． ここで，λ は R の固有値であり，λ ≥ λ ≥ · · · ≥i,kii,1i,2
ここで，round(·) は整数への丸め処理を表す．なお心拍は λのように降順に並んでいるものとする．また，v
一般的に beat per minute (bpm) で表現されるため，推定心拍は ĥo [i] = 60fˆo [i] となる．

3.3i,Ki,kは λi,k に対応する固有ベクトルである．いま，短時間区間信号 (xs [i, l])Ll=1 が次式のように表されると仮定する．

MUSIC 法に基づく心拍推定

xs [i, l] = s[i, l] + ρ[i, l]MUSIC 法 [3] は部分空間法に基づく信号の分類手法で

s[i, l] =あり，センサーアレイ信号処理で信号の到来方向推定によく用いられる．一方で，時間信号に MUSIC 法を適用し基本周波数を推定する手法 [4], [5] も提案されており，ノイズに対して頑健という特徴がある．本稿では，この MUSIC法による基本周波数推定を PPG からの心拍推定に応用する．
まず，短時間区間信号 (xs [i, l])L
l=1 中の長さ K （但しN∑(12)

An [i] sin(2πnfo [i]l + θn [i])

(13)

n=1

ここで，An [i] 及び θn [i] は n 次基本成分の振幅及び位相差Lをそれぞれ表し，(ρ[i, l])Ll=1 は (s[i, l])l=1 と無相関なノイズ信号である．即ち，観測信号は基本周波数 fo [i] とその倍音 (nfo [i])Nn=2 の周波数をもつ正弦波から成る目的信号
L(s[i, l])Ll=1 とノイズ信号 (ρ[i, l])l=1 の混合と仮定する．も∑L∑L
し各信号のパワー比が l=1 |s[i, l]|2 ≫ l=1 |ρ[i, l]|2 で

K < L）の部分信号を切り出して次式のようにベクトルあれば，式 (9) で得られる固有値のうち，上位 2N 番目まで定義する．

2Nでの固有値 (λi,n )2N
n=1 に対応する固有ベクトル (vi,n )n=1


xs [i, 1 + (m − 1)P ]

..

.


K

xi,m =  xs [i, k + (m − 1)P ] 
∈R


..
.

xs [i, K + (m − 1)P ]

が目的信号 (s[i, l])L
l=1 の部分空間を張り，残りの固有値K(λi,n )K
n=2N +1 に対応する固有ベクトル (vi,n )n=2N +1 がノ2
イズ信号 (ρ[i, l])L
l=1 の部分空間を構成する ．従って，式(7)(9) は次式のように分解できる．
(s)

ここで，m = 1, 2, · · · , M は xs [i, l] 中の部分信号インデクス，k = 1, 2, · · · , K は部分信号内の離散時間インデク
ス，P は部分信号を切り出す際のシフト長をそれぞれ表す．このベクトルより，i 番目の短時間区間の信号相関行
1 より厳密には，log |F [i, ω]| は ω に関して偶関数であるため，実数
ケプストラム C[i, q] は log |F [i, ω]| の cos 展開係数といえる．

(s)(s) −1

Ri = Ui Di Ui

λi,10

λi,2 0(s)Di = ..
 ..
 .
.
0···

(ρ)

(ρ)

(ρ) −1+ Ui Di Ui
···0
···0 .. ..
.
. 
0(14)

(15)

λi,2N

2 目的信号に対応する固有値の数が正弦波の数 N の 2 倍となるのは，式 (13) において実数正弦波の重ね合わせを仮定しているためである．
複素正弦波の重ね合わせであれば，対応する固有値の数と正弦波個数は一致する [5]．

‑ 7 ‑

(s)

Ui

[=

ui,2Nλi,2N +10···00..
.

λi,2N +2..
.

···
..
.

0..
.

0···

0λi,Kui,2N +1ui,2N +2···ui,K(ρ)
Di = (ρ)

Ui[=]···

ui,1ui,2(16)拍推定を適用し，その精度について議論する．具体的な



測定条件は文献 [1] を参照されたい．なお，今回の実験で

は心拍の真値は定かではない問題があるが，PPG の解析(17)

結果及び PPG と同時に測定した ECG の解析結果がどの程度整合するかを確認し，付け爪型センサから得られるPPG の有用性について議論する．但し，ECG による心](18)

拍測定は通常仰臥位かつ平静状態で記録するものである
ため，Fig. 3 に示すような日常的な活動下で測定された

理想的には，式 (14) の第一項及び第二項がそれぞれ目的

ECG による推定心拍も多分に誤差が含まれている点に注信号の信号相関行列及びノイズ信号の信号相関行列に対意する必要がある．

応する．

心拍推定においては，PPG 及び ECG のサンプリング

MUSIC スペクトルは，ノイズ信号の固有ベクトル(ρ)(ui,k )K
k=2N +1 （Uiの列ベクトル）を用いて次式で計算される．

周波数 fs を測定時の 1 kHz から 20 Hz にダウンサンプリングして解析した．また，心拍の推定範囲は 35–210 bpm（fmin = 0.58 Hz 及び fmax = 3.5 Hz）に制限した．短時

1

G[i, γ] = ∑K

T2
k=2N +1 |ui,k b[γ]|
12π(γ−1)



e−j K


4π(γ−1)−jK
eb[γ] =  ∈ CK


.

..

2π(K−1)(γ−1)
K
e−j(19)間信号長は L = 1024（51.2 s）とした．ケプストラム分
析では，短時間区間のシフト長を S = L/4，DFT の窓関数 (w[l])L
l=1 をハミング窓に設定した．また MUSIC 法では，短時間区間のシフト長を S = L（オーバーラップ無(20)し），式 (7) の部分信号長を K = L/2，部分信号のシフ
ト幅を P = K/32，仮定する目的信号中の正弦波の個数を N = 1 に設定した．

ここで，γ = 1, 2, · · · , Γ は離散周波数インデクスである（但し Γ = K ）．b[γ] は DFT 基底であるため，uTi,k b[γ]

4.2

推定結果

測定された 24 時間の PPG 及び ECG の解析結果をそはノイズ信号の固有ベクトル ui,k の DFT そのものであれぞれ Figs. 4 及び 5 に示す．なお，文献 [1] には各時刻る．式 (19) はノイズ信号の固有ベクトルの DFT が分母における動作の詳細があるため，合わせて参照されたい．にあることから，MUSIC スペクトル (G[i, γ])Γ
γ=1 は，目的信号 (s[i, l])L
l=1 の部分空間に対応する周波数（式 (13)における正弦波周波数 (nfo )Nn=1 ）で顕著なピークを持つ．

まず，Figs. 4 及び 5 の上端に示したスペクトログラムに着目する．これらを見ると，基本周波数とその倍音成分と思われる目的信号成分の軌跡を確認できる．しかし，活動従って，ノイズ信号に対して頑健に基本周波数の推定が状態によっては波形がクリッピングし，全周波数にエネル可能となる．

ギー現れて目的信号成分が確認できなくなる．特に，PPG短時間区間毎の基本周波数の推定値 (fˆo [i])Ii=1 は，次式のように MUSIC スペクトル (G[i, γ])Γ
γ=1 が最大値を取る
周波数インデクス (γ̂o [i])Ii=1 から求められる．

fsfˆo [i] = (γ̂o [i] − 1)Kγ̂o [i] = arg max G[i, γ]はその割合が多く，例えば 22:00–23:00 及び 08:00–12:00
はノイズ信号の成分が支配的である．ECG は，PPG よりも SN 比が良いが，それでもいくつかの時間帯で目的信号の成分が確認できなくなる．なお，14:00–15:00 付近は(21)(22)

1.3 Hz 程度に強い基本周波数成分が現れているが，この時間帯はランニング運動をしており，これは心拍ではな

γ

く体動に由来するノイズ信号の基本周波数成分と思われケプストラム分析の場合と同様に，実際は推定範囲の下る．このランニング運動中の実際の心拍は 2.0 Hz 以上と限 fmin と上限 fmax に対応する離散周波数インデクス予想され，ECG のスペクトログラムにはかすかに目的信の下限 γmin = round(Kfmin /fs + 1) 及び上限 γmax =号と思しき成分が 2.5 Hz 付近に確認できる．
次に，PPG を用いたケプストラム分析及び MUSIC 法

round(Kfmax /fs +1) の範囲内で (22) の最大値探索を行う．に基づく心拍推定 (ĥo [i])Ii=1 の結果（Fig. 4 の中央と下端）

実験

44.1を比較する．ケプストラム分析での推定心拍は値が激し実験条件く変動しているが，MUSIC 法の推定心拍は比較的穏やか

本稿では，2.3 節で示した PPG 及び ECG の測定波形 [1]である．特に 01:00–07:00 の睡眠時は MUSIC 法においてに対して，ケプストラム分析及び MUSIC 法に基づく心50–60 bpm といった尤もらしい心拍が高精度で推定でき‑ 8 ‑Frequency [Hz]Spectrogram (PPG signal)
3210

Quefrency [s]Cepstrogram (PPG signal)
3020
101Heart rate [bpm]Estimated heart rates based on cepstrum analysis (PPG signal)
200
150
100500Frequency [Hz]MUSIC spectrogram (PPG signal)
3210Heart rate [bpm]Estimated heart rates based on MUSIC method (PPG signal)
200
150
10050
0
18:00

22:00

02:00

06:00

10:00

14:00

18:00

Japanese standard time

Fig 4: Analysis results using observed PPG, where spectrogram (top raw), cepstrogram (second raw), estimatedheart rates based on cepstrum analysis (third raw), MUSIC spectrogram (fourth raw), and estimated heart rates
based on MUSIC method (bottom raw) are depicted.ていることが確認できる．ただし，MUSIC 法は通常活動

通常活動時及び運動時の推定心拍が PPG よりも高精度で時及び運動時の推定心拍が安定せず，ケプストラム分析あるが，逆に睡眠時は PPG の推定心拍よりも不安定であの方が（スパイク状の変動を無視すれば）辛うじて心拍る．MUSIC 法の推定心拍も全体的に PPG の場合より不の変動のような軌跡が確認できる．この違いは，ケプス安定である．これらの原因は，Fig. 3 (b) に示すような目トログラム（Fig. 4 の上から 2 番目）と MUSIC スペクト的信号の波形の違いによるものと思われる．ケプストラム

ログラム（Fig. 4 の上から 4 番目）でも確認でき，通常活分析と MUSIC 法はいずれも，目的信号に対して式 (13)動時及び運動時はケプストログラムのピーク値が，睡眠のような倍音構造を仮定しているため，心拍由来の信号が時は MUSIC スペクトルのピーク値が比較的精度良く現スパイク状になる ECG よりも，より正弦波に近い PPGれている．なお，いずれの推定心拍も，移動中央値フィルの方が高精度な心拍推定ができたものと推察される，

タ等を適用することでスパイク状の変動を除去でき，よ5り尤もらしい滑らかな心拍の時間変化が得られる．

まとめ

最後に，ECG に対する両手法の推定心拍（Fig. 5 の中

本稿では，付け爪型センサから得られる PPG にケプス央と下端）について述べる．ケプストラム分析の結果はトラム分析及び MUSIC 法に基づく基本周波数推定を適用‑ 9 ‑Frequency [Hz]Spectrogram (ECG signal)
3210

Quefrency [s]Cepstrogram (ECG signal)
3020
101Heart rate [bpm]Estimated heart rates based on cepstrum analysis (ECG signal)
200
150
100500Frequency [Hz]MUSIC spectrogram (ECG signal)
3210

Heart rate [bpm]Estimated heart rates based on MUSIC method (ECG signal)
200
150
100
500
18:00

22:00

02:00

06:00

10:00

14:00

18:00

Japanese standard time

Fig 5: Analysis results using observed ECG, where spectrogram (top raw), cepstrogram (second raw), estimatedheart rates based on cepstrum analysis (third raw), MUSIC spectrogram (fourth raw), and estimated heart rates
based on MUSIC method (bottom raw) are depicted.
し，心拍の推定精度について議論した．結果より，PPG を

Digital Signal Processing, IEEE Press, New York,用いた心拍推定がある程度可能であることを確認した．ま1979.

た，ケプストラム分析は通常活動時及び運動時，MUSIC[3] R. O. Schmidt, “Multiple emitter location and sig-

法は睡眠時の PPG に対して高精度な推定心拍を与えるこ

nal parameter estimation,” IEEE Trans. Antennasとが確認された．

Propag., vol. AP-34, no. 3, pp. 276–280, 1986.
[4] M. G. Christensen, S. H. Jensen, S. V. Andersen,

参考文献

and A. Jakobsson, “Subspace-based fundamental[1] K. Ishii and N. Hiraoka, “Nail tip sensor: toward reliable daylong monitoring of heart rate,”
Trans. Electrical and Electronic Eng., vol. 15, no. 6,frequency estimation,” Proc. EUSIPCO, pp. 637–640, 2004.
[5] M. G. Christensen, A. Jakobsson, and S. H. Jensen,

pp. 902–908, 2020.

“Joint high-resolution fundamental frequency and[2] IEEE Acoust., Speech, and Signal Process. Society,order estimation,” IEEE Trans. Audio, Speech, andDigital Signal Process. Committee, Programs forLang. Process., vol. 15, no. 5, pp. 1635–1644, 2007.

‑ 10 ‑