第35回 回路とシステムワークショップ
The 35th Workshop on Circuits and Systems
in Kitakyushu & on online, August 19-20, 2022

レーダセンサと信号源分離アルゴリズムを用いた運転者の心拍推定

Estimation of Driver’s Heart Rate Using Radar Sensors
and Source Separation Algorithms
村田佳斗 †

Keito MURATA†

1

北村大地 †齋藤諒 ‡† 香川高等専門学校‡ 株式会社村田製作所

植木大地 ‡

Daichi KITAMURA†
Ryo SAITO‡†National Institute of Technology, Kagawa College‡Murata Manufacturing Co., Ltd.

Daichi UEKI‡

はじめに自動車の運転中に運転者が発作や意識喪失等に見舞わ
Ch. 12
Ch. 3
Ch. 4
Ch.

れることは致命的な状況となる．そのため，運転中に運転者の状態を何らかの方法で管理することは重要課題の一つとなっている．この課題に取り組むために，本稿ではFig. 1 (a) に示すようなレーダ非接触型生体センサアレイ（以後，レーダセンサと呼ぶ）を用いて，運転中の運転者の心拍を測定及び推定するシステム（以後，振動測定系と呼ぶ）を取り扱う．
Fig. 1 (b) に振動測定系の概観を示す．シートの背もたれ内部にレーダセンサを埋め込むことで運転者の背部の微小変位が測定できるため，心拍由来の体表面変位から心拍を推定することを目指す．なお，運転中の系全体の振動を模擬するために，シートは振動台の上に設置されている．このような振動測定系では，心拍由来の体表面変位以外にも系全体（振動台）の振動，運転者の体動，呼吸による体表面変位等の成分が同時に観測されてしまう．
本稿では，前述の振動測定系のレーダセンサで得られる観測信号に対してブラインド信号源分離（blind source(a) Beams from radar sensor(b) Vibration system

Fig. 1: Vibration measurment system with radar sensorattached to back of driver’s seat.

1.2 Hz の正弦波で上下方向に振動させた．レーダセンサはシートの背もたれ内部に埋め込まれているため，Fig. 1 (a)のように被験者の背部の微小な体表面変位の時間変化が測定される．また，レーダセンサは 4 つの異なる指向性で人体に照射しており，同時に近傍 4 点の体表面変位を4 ch の信号として得ることができる．このレーダセンサのサンプリング周波数は 40 Hz である．なお，被験者には接触型心電センサ（以後，ECG センサと呼ぶ）を胸部に貼り付けて測定しており，このセンサから得られる信号を解析して計算された心拍を参照値に利用する．

separation: BSS）[1, 2] を適用し，心拍由来の体表面変位を抽出することで，心拍の推定精度の向上を目指す．特に，音響信号処理分野で代表的な BSS である独立ベクトル分 2.2 レーダセンサの観測スペクトログラム
析（independent vector analysis: IVA）[3, 4] 及び独立低ランク行列分析（independent low-rank matrix analysis:
レーダセンサの観測信号（Ch. 1）及び ECG センサのILRMA）[5, 6] を適用し，心拍の推定精度を比較する．
観測信号のスペクトログラムを Fig. 2 に示す．レーダセンサの観測信号には心拍由来の成分以外にも，振動台の2 振動測定系と観測信号の概要振動成分及び呼吸に起因する成分等が混在していること
2.1 振動測定系の測定条件が確認できる．一方で，ECG センサの観測信号にはほとFig. 1 (b) の振動測定系では，被験者が座った状態で測 んど心拍由来の成分のみ含まれている．本稿では，レー定を行った．測定は全体で 420 s 間行い，途中の 60 s から ダセンサの観測信号に BSS を適用し，ECG センサの観360 s の時間区間でのみ，振動台を振幅 10 mm，周波数 測信号のような心拍成分のみを抽出することを目指す．

－ 96 －Vibration
system Observed signal (1 ch) Heartbeat

Demixing Observed Estimatedsignalmatrix
signal

ECG signal Heartbeat(a) Mixing model

Breathingmotion

Fig. 3: Mixing and demixing models (M = N = 2).

Fig. 2: Spectrograms of observed signal obtained byradar sensor (left) and ECG signal (right).

おいて全周波数ビンの成分をまとめたベクトルであり，次式で定義される．

sjn = [s1jn , s2jn , · · · , sijn , · · · , sIjn ]T ∈ CI

BSS 及び心拍推定アルゴリズム

3
3.1(b) Demixing model

BSS の定式化

y jn = [y1jn , y2jn , · · · , yijn , · · · , yIjn ]T ∈ CI

本章では，本稿で適用する時間周波数領域 BSS の定式化を示し，IVA 及び ILRMA のアルゴリズムを説明する．さらに，本稿で適用する心拍推定アルゴリズムについても述べる．BSS では，心拍，呼吸，振動台等の複数の信号源が混合されて観測される状況を仮定する．この信号源を s̃[l] =[s̃1 [l], · · · , s̃n [l], · · · , s̃N [l]]T ∈ RN と定義する．ここで，
l = 1, 2, · · · , L 及び n = 1, 2, · · · , N はそれぞれ離散時間及び信号源のインデクスである．同様に，観測信号を x̃[l] =[x̃1 [l], · · · , x̃m [l], · · · , x̃M [l]]T ∈ RM ，BSS を適用して得られる分離信号を ỹ[l] = [ỹ1 [l], · · · , ỹn [l], · · · , ỹN [l]]T ∈RN とそれぞれ定義する．ここで，m = 1, 2, · · · , M は観測チャネルのインデクスである．時間周波数領域 BSS では，時間信号に短時間 Fourier 変換（short-time Fouriertransform: STFT）を適用して得られる時間周波数領域の信号を取り扱う．s̃[l]，x̃[l]，及び ỹ[l] の時間周波数領域の信号をそれぞれ次式で定義する．
sij = [sij1 , sij2 , · · · , sijn , · · · , sijN ]T ∈ CN
xij = [xij1 , xij2 , · · · , xijm , · · · , xijM ] ∈ CT

yij = [yij1 , yij2 , · · · , yijn , · · · , yijN ] ∈ CTMN(7)信号源数がチャネル数と等しく（N = M ）かつ Ai が可逆な場合は，分離行列 Wi = [wi1 wi2 · · ·C

N ×M

wiN ]H ∈が存在し，分離信号を次式で表せる．

yij = Wi xij(8)

ここで，·H は行列及びベクトルのエルミート転置を表す．
この分離系を Fig. 3 (b) に示す．時間周波数領域 BSS は，混合行列 Ai が未知の条件下で，観測信号 xij のみから分離行列 Wi（≈ A−1i ）を推定し，信号源 sij の推定値である分離信号 yij を得ることを目的とする．
本稿では，レーダセンサから M = 4 チャネルの観測信号が得られるため，混合している信号源も N = 4 つと仮定して BSS を適用する．このとき，BSS で得られる 4 つの分離信号の内，どの信号が心拍信号を多く含む推定心拍信号であるかは最適化アルゴリズムの初期値に依存し(2)択する方法については今後の課題とし，実験では全て手(3)動で ECG センサの信号に近いスペクトログラムを持つ分離信号を推定心拍信号として選択している．

た，式 (1)–(3) の各信号においては，時間周波数行列（スI×J

3.2IVA

本稿で適用する 1 つ目の時間周波数領域 BSS は IVA [3,， 4] と呼ばれる．IVA は「信号源間の統計的独立性」と「各信号源の全周波数成分の共起性 [3]」の両方を仮定して，

Xm ∈ C，及び Yn ∈ Cとして定義しておく．
時間周波数領域 BSS では，複数の信号源が混合して観測される現象を，周波数毎の時不変な瞬時混合でモデル化する．即ち，観測信号は周波数毎の混合行列 Ai ∈ CM ×Nを用いて次式で表せると仮定する．
I×J(6)て変化する．本稿では，分離信号から推定心拍信号を選れ周波数ビン及び時間フレームのインデクスを表す．まペクトログラム）としての表記もそれぞれ Sn ∈ CI(1)

ここで，i = 1, 2, · · · , I 及び j = 1, 2, · · · , J はそれぞI×J(5)

xjm = [x1jm , x2jm , · · · , xijm , · · · , xIjm ] ∈ C
T

観測信号のみから分離行列 Wi を推定する．各信号源の全周波数成分の共起性とは，Fig. 4 (a) に示すように，分離信号の周波数ベクトル y jn 内の全成分が同じ時刻 j で大きな値をとることを仮定したものである．例えば，Fig. 2の心拍成分や振動台の成分には明確な調波構造が見られ，

xij = Ai sij(4)

基本周波数とその整数倍の成分が共起していることから，この混合系を Fig. 3 (a) に示す．ここで，図中の sjn ，xjn ， Fig. 4 (a) の IVA の信号源モデルがある程度成立するこ及び y jn はそれぞれ信号源，観測信号，及び分離信号にとが分かる．

－ 97 －

Frequency

源と排他的に近似されるため，信号源分離が促進される．

ce

前述のように，心拍成分や振動台の成分は時間的に連続Time(b) ILRMASour

Time(a) IVASource

Frequency

繰り返す．この NMF によるモデル化の際に，他の信号ILRMA における分離行列 Wi の推定は，次式の最小化

IVA における分離行列 Wi の推定は，次式の最小化問題として定式化される [3, 4]．
XX
Minimize − 2J
log | det Wi | +∥y jn ∥2(9)i

モデルは IVA の信号源モデルよりも，さらに緻密に各信号源の時間周波数構造を表現できることが予想される．

Fig. 4: Source models assumed in IVA and ILRMA.

Wな調波構造が現れるため，Fig. 4 (b) の ILRMA の信号源問題として定式化される [5, 6]．

Minimize − 2J
W,T,VX

log | det Wi | +X

D |Yn |.2 |Tn Vnni(13)N
ここで，T = {Tn }N
n=1 及び V = {Vn }n=1 は NMF 信号j,n

ここで，W = {Wi }Ii=1 は最適化変数の集合であり，∥ · ∥2I×K
源モデルに関する最適化変数の集合を表し，Tn ∈ R≥0は L2 ノルムを表す．式 (9) の最小化問題は，解析的な解は及び Vn ∈ RK×Jは NMF における非負基底行列及び非≥0

未発見であるが，補助関数法 [7] 及び反復射影法（iterative

負係数行列である．さらに，D(·|·) は 2 つの引数行列の要projection: IP）に基づく効率的な反復最適化アルゴリズ
ムが補助関数 IVA [4] として提案されている．この手法では，分離行列 Wi に何らかの初期値を与えた上で，次式の計算を収束するまで反復する．
11XqP
Gin =xij xH(10)ij
J j
H x |2|w
in ij
i

素間の板倉斎藤擬距離 [8] の総和であり，K > 0 は各信

win ← (Wi Gin )−1 en(11)1Hwin ← win (winGin win )− 2(12)ここで，en ∈ {0, 1}N は n 番目の要素のみが 1，他要素号源のパワースペクトログラムを低ランク近似する際のNMF の基底ベクトル数である．式 (13) の最小化問題は，IVA と同様に反復最適化アルゴリズムで解くことができる．ILRMA の反復最適化アルゴリズムは，分離行列 Wi及び NMF の行列 Tn 及び Vn に何らかの初期値を与えた上で，次式の計算を収束するまで反復する．
vPuP−2H x |2 v
u j |win
ijkjn (k′ tik′ n vk′ jn )ttikn ← tiknPP
−1j vkjn (k′ tik′ n vk′ jn )sPが 0 のベクトルである．なお，この反復最適化アルゴリ
ズムは，1 回の更新で目的関数 (9) の値が減少する又は変vkjn ← vkjn

(14)
H2i |win xij | tikn (P動しないこと（単調非増加）が理論的に保証されている．

i tikn (PPk′ tik′ n vk′ jn )−1k′ tik′ n vk′ jn )−2(15)
3.3

ILRMA
Uin =

本稿で適用する 2 つ目の時間周波数領域 BSS は IL-RMA [5, 6] と呼ばれる．ILRMA は「信号源間の統計的独立性」と「各信号源の時間周波数構造の低ランク性」の両方を仮定して，観測信号のみから分離行列 Wi を推定する．
各信号源の時間周波数構造の低ランク性とは，Fig. 4 (b)に示すように，分離信号のパワースペクトログラム |Yn |.2が低ランク行列で良く近似できる（同じスペクトルパターンの繰り返しを多く含む）ことを仮定したものである．ここで，行列に対する絶対値記号及びドット付き指数乗はそれぞれ要素毎の絶対値及び指数乗を表す．この低ランクモ
デルには，非負値行列因子分解（nonnegative matrix factorization: NMF）[8] を用いている．Fig. 5 は，ILRMAによる BSS の原理を示した図である．ILRMA は，周波数毎の分離行列 Wi で分離された信号源のパワースペク
トログラム |Yn |.2 を NMF で近似的にモデル化し，そのモデルを用いて再度分離行列を更新するという手続きを

1X
JjP

1xij xH
ijtvikn
kjnkwin ← (Wi Uin )−1 enHwin ← win (winUin win )

(16)
(17)

− 21(18)

ここで，tikn 及び vkjn はそれぞれ Tn 及び Vn の要素であり，k = 1, 2, · · · , K は基底ベクトルのインデクスである．
この反復最適化アルゴリズムも，IVA と同様に，目的関数 (13) の単調非増加性が保証されている．

3.4

心拍推定アルゴリズム

本節では，BSS で推定された推定心拍信号に適用する心拍推定アルゴリズムについて説明する．BSS で得られる時間領域の推定心拍信号を (ỹHR [l])Ll=1 と表記する．ECG
センサ等で得られる心拍信号から心拍を推定する手法は，

R 波（心拍信号に現れるスパイク状の波形）と次の R 波までの間隔（R-R interval: RRI）を用いることが一般的－ 98 －Fig. 5: BSS principle of ILRMA.

Emphasize
heartbeat

ログラムである2 ．また，IVA を適用して得られる 4 つのBand passfilter

分離信号のスペクトログラムを Fig. 8 に示す．この分離信号の中の推定心拍信号が，Fig. 2 の右側に示した ECG
センサで得られる心拍信号のスペクトログラムにどの程度近いかが重要となる．IVA の分離結果を見ると，3 番目

Contourtracking ofspectral peaksの分離信号に心拍の高調波成分が強く見られる（3.5 Hz，Heart rateestimates

Fig. 6: Process ﬂow of heart rate estimation algorithm.
である．しかしながら，レーダセンサと BSS で得られる推定心拍信号には明確な R 波が含まれないため，RRI の計算が困難である．そこで，本稿では Fig. 6 のような心拍推定アルゴリズムを適用する．心拍推定アルゴリズム
ではまず，信号の調波構造を強調しつつ低周波成分を除去し，BPF によって心拍成分を取り出す．そして，得られた心拍成分に STFT を適用してスペクトログラムに変換し，心拍スペクトログラムの最大ピークとなる周波数を時間フレーム毎に求めることで，推定心拍を得る．
フィルタを適用しない場合の BSS 及び心拍推定実験

44.1

実験条件

本章では，まず 2 章に示した観測波形にそのまま IVAを適用して得られる分離信号のスペクトログラム及び推5 Hz，及び 6 Hz 付近）一方で，0.5 Hz 付近には，呼吸による体動に起因する成分が非常に強く残留しており，心拍以外の成分を完全には除去できていない．振動台の成分の多くは他の分離信号に現れているが，1.2 Hz の基本周波数成分は 3 番目の分離信号に多く残留してしまっていることが確認できる．
Fig. 8 の 3 番目の分離信号に心拍推定アルゴリズムを適用した結果を Fig. 9 に示す．ECG センサから得られる信号（Fig. 2 の右側）から計算された心拍の参照値も重ねて示している．この結果より，振動台の影響が無い 0–60 s及び 360–420 s の区間は，推定心拍が参照値と概ね一致していることが分かる．一方，振動台が影響する 60–360 sの区間は多くの推定誤差を含む結果となった．
以上の結果より，現状の方法では呼吸由来の成分が全ての分離信号に強く残留してしまうことが明らかとなった．また，振動台の成分も完全には除去できず，心拍推定の結果に悪影響を及ぼしたことが推測される．次章の実
験では，BSS の前処理としてハイパスフィルタを適用することで呼吸の成分を除去した場合の実験結果を示す．

定心拍信号に 3.4 節の心拍推定アルゴリズムを適用した結果を示す．IVA を適用するために，観測信号に適用する5フィルタを適用した場合の BSS 及び心拍推定実験

STFT は Fourier 変換長（窓長）を 1.6 s，窓関数をハミ 5.1 前処理として適用するフィルタの条件
ング窓，窓のシフト長を 0.1 s に設定した1 ．IVA におけ
本章では，観測信号に対して呼吸の基本周波数成分る分離行列 Wi の初期値は，全ての周波数に対して単位 （0.5 Hz 付近）を除去するハイパスフィルタを BSS の前行列とした．また，IVA の反復最適化アルゴリズムの反 処理として適用する．本実験では，カットオフ周波数を復回数は 100 回に設定した．
1.5 Hz，次数を 170 に設定した有限インパルス応答（ﬁnite
impulse response: FIR）フィルタを用いた．
4.2 実験結果
5.2 実験条件
Fig. 7 はレーダセンサから得られる観測信号のスペクト観測信号に適用する STFT の条件及び IVA の条件は 4.1
1
この BSS のための STFT の設定は，本稿の結果で掲載しているスペクトログラムを得る際の STFT の設定とは大きく異なる．後者の設定は成分の見やすさのために，全て窓長を 25.6 s，シフト長を 0.8 s に設定している．

2Fig. 7 の 1 番目のチャネルのスペクトログラムは Fig. 2 の左側と同じものである．

－ 99 －Fig. 7: Four-channel spectrograms of observed signal obtained by radar sensor.

Fig. 8: Spectrograms of separated signals estimated byIVA.
Reference(ECG signal)

節の記載内容と同様である．また，本実験では ILRMA を適用した際の結果も示す．ILRMA では，IVA と同様に分Heart rate [bpm]離行列 Wi の初期値を単位行列とし，NMF の行列 Tn 及
び Vn の初期値は区間 (0, 1) の一様乱数とした．また，各信号源の基底ベクトル数を K = 3 とし，反復最適化アルゴリズムの反復回数は 100 回に設定した．

5.3実験結果

BSS の前処理にハイパスフィルタを適用した場合の IVA及び ILRMA で得られる分離信号のスペクトログラムをそれぞれ Figs. 10 及び 11 に示す．いずれも 3 番目の分離信号に心拍の成分が多く含まれており，3 番目が推定心拍信号と考えられる．まず IVA の結果（Fig. 10）では，推定心拍信号に振動台の成分がほとんど含まれておらず，心拍の高調波が 3.5 Hz，5 Hz，及び 6 Hz 付近に強く現れている．さらに，前処理で 1.5 Hz 以下の成分を除去したことから，0.5 Hz 付近の呼吸の成分も大幅に減衰している．次に ILRMA の結果（Fig. 11）では，推定心拍信号に対して IVA とほぼ同様のことがいえる他，心拍の高調波が IVA の場合よりも大きなエネルギーで残留していることが分かる．従って，いずれの BSS においても，レーダセンサの観測信号から心拍由来の成分を高精度に抽出できる可能性が示唆された．
Figs. 10 及び 11 の 3 番目の分離信号に対して，心拍推定アルゴリズムを適用した結果を Figs. 12 及び 13 にそれぞれ示す．これらの結果と 4 章の結果（Fig. 9）を比較すると，振動台が影響する 60–360 s の区間の推定精度が大幅に向上していることが分かる．特に，IVA よりも ILRMA

Radar sensor
and BSS

Fig. 9: Estimated (red) and reference (blue) heart rates
obtained by IVA.の方が高精度な心拍推定を実現している．ILRMA を適用した場合の心拍推定結果では，60–120 s の区間のみが参照値から外れている．この時間区間は振動台の加振の初期段階であり，被験者の身体が大きく動いてしまったことが参照値と外れた要因と推察される．
結論6

本稿では，運転者の心拍のモニタリングを目的とした振動測定系を構築し，実際にレーダセンサから得られる観測信号に対して心拍を推定する方法について検討した．
提案手法では，心拍由来の体表面変位のみを抽出するために IVA 及び ILRMA を適用し，その分離信号から心拍を推定した場合の精度を比較した．実験結果より，BSS はレーダセンサの観測信号に対しても効果があることを実証－ 100 －Fig. 10: Spectrograms of separated signal obtained by
IVA with high-pass ﬁltering.

Fig. 11: Spectrograms of separated signal obtained by
ILRMA with high-pass ﬁltering.

した．また，呼吸の成分を前処理で除去することで，BSSReference(ECG signal)

参考文献

Radar sensorand BSS

Heart rate [bpm]及び心拍の推定精度を大幅に改善できることを示した．

[1] P. Comon, “Independent component analysis, A newconcept?,” Signal Processing, vol. 36, no. 3, pp. 287–
314, 1994.
[2] H. Sawada, N. Ono, H. Kameoka, D. Kitamura, andH. Saruwatari, “A review of blind source separation
methods: Two converging routes to ILRMA originating from ICA and NMF,” APSIPA Trans. Signal andInformation Processing, vol. 8, no. e12, pp. 1–14, 2019.
[3] T. Kim, H. T. Attias, S.-Y. Lee, and T.-W. Lee, “Blindsource separation exploiting higher-order frequency dependencies,” IEEE Trans. Audio, Speech, and Language
Processing, vol. 15, no. 1, pp. 70–79, 2007.
[4] N. Ono, “Stable and fast update rules for independentvector analysis based on auxiliary function technique,”Proc. IEEE Workshop on Applications of Signal Processing to Audio and Acoustics, pp.189–192, 2011.

Fig. 12: Estimated (red) and reference (blue) heart rates
obtained by IVA with high-pass ﬁltering.

[5] D. Kitamura, N. Ono, H. Sawada, H. Kameoka,and H. Saruwatari, “Determined blind source separation unifying independent vector analysis and nonnegative matrix factorization,” IEEE/ACM Trans. Audio, Speech, and Language Processing, vol. 24, no. 9,
pp. 1626–1641, 2016.
[6] D. Kitamura, N. Ono, H. Sawada, H. Kameoka,and H. Saruwatari, “Determined blind source separation with independent low-rank matrix analysis,” Audio Source Separation, S. Makino, Ed., pp. 125–155.
Springer, Cham, 2018.
[7] D. R. Hunter and K. Lange, “A tutorial on MM algorithms,” The American Statistician, vol. 58, no. 1,
pp. 30-–37, 2004.
[8] C. Févotte, N. Bertin, J.-L. Durrieu, “Nonnegativematrix factorization with the Itakura-Saito divergence.
With application to music analysis.” Neural Computation, vol. 21, 793–830, 2009.

Fig. 13: Estimated (red) and reference (blue) heart rates
obtained by ILRMA with high-pass ﬁltering.

－ 101 －