# 肺音について
生体信号は大量にデータを入手できるとは限らず、生体情報の正解が不明な信号の例しか集められないことも多い。一方で、生体信号には何らかの典型的な特徴を持っていることが多く、その特徴の組み合わせや重ね合わせで生体情報が表現されていると考えられる。

聴診は呼吸器系の異常を早期発見したり病態を評価する基本的な診察手技の一つである。聴診音は数理モデル化できるほど発生・伝播メカニズムが単純ではなく特徴を捉え難い。また、聴診する医師や看護師には聴診音から健康状態や病状を推測する知覚情報処理能力が求められる。聴診を身につけるためには音の特徴を捉えるための修練が必要なのである。

この肺聴診を支援するパターン認識とその実現のために有効であると考えられる信号処理を課題とし考察する。
***
私たちは工学部の人間なので、病気かどうかの診断装置を作るわけではない。
あくまで医学的な可能性を示唆する・補助するような情報を提供する。
***
### 目的
聴診音を構成する個々の肺音に分類するだけでなく、その特徴となる波形やパワースペクトルをデータとして記録、可視化することで医師や看護師が聴診をする際の技術的支援を目的とする。また、それだけでなく、人間では見つけることができないような肺音特徴の発見も期待できる。

digraph hierarchy {
nodesep=0.5 // increases the separation between nodes

node [color=Red,fontname=Courier,shape=box]
edge [color=black, style=dashed]

聴診音->{肺音 雑音（ノイズ）}
肺音->{呼吸音 副雑音}
副雑音->{胸膜摩擦音や肺血管性雑音等 ラ音}
ラ音->{連続性ラ音 断続性ラ音}
連続性ラ音->{類鼾音 笛声音 スクウォーク}
断続性ラ音->{水泡音 捻髪音}
{rank=same;類鼾音 笛声音 スクウォーク 水泡音 捻髪音}
}

### 背景
聴取できる種々の肺音は今の実臨床では、「たまたま、ある瞬間に、ある医者が、聴いたときの、その医者の耳の調子で、決まってしまっているのが現状」というように病理学的、聴覚的な比喩で命名されている。肺音を信号処理するために、電子デバイスで身体から音を計測する電子聴診器によって肺音の記録や伝送が可能である。また、肺聴診音は呼吸音・副雑音・ノイズの混合音と見なすことができ、この聴診音のノイズ除去や各信号への分離といった信号処理は医師や看護師の聴診技術の向上に貢献するだけでなく、コンピュータによる聴診支援システムの開発を検討できるようになるため有益である。肺音認識は与えられた入力信号の任意の時間区間に複数のクラスラベルを付与する多重ラベリングの識別問題であると言える。聴診音は呼吸音とノイズが同時に聴取されたり、様々な特性の肺音が時間領域・周波数領域で重なって聴取される場合もあり、純粋にウェーブレット変換、フーリエ変換によって聴診音をフィルタリングするだけでは識別的な特徴を得られることができないため難しい課題とされている。さらに課題はまだ存在する。肺音の波形やスペクトルは聴診器や計測方法、疾患によって大きく異なるため各肺音の信号モデルを設計することが容易でない。

### アプローチ
これまで、肺音の信号処理技術では"肺音の特徴抽出"、"異常の検出"、"識別"が試みられている。また、音声的、画像処理的なアプローチも検討されてきた。しかし、これらの信号処理では聴診音が異なる種類の肺音の混合信号であることが考慮されていないため、各種肺音の特徴を明確に抽出・分離できていない。
上記の課題を考慮し、聴診音信号をいくつかの少数の信号成分の組み合わせとして効率的に表現できる性質（スパース性）を利用して異なる特性を持つ肺音成分に分離する。

start
:入力信号;
if (ロバスト主成分分析(RPCA)[スペクトログラムを見る]u-netとか) then (スパース行列)
:連続性ラ音のスペクトログラム;
split
:類鼾音;
kill
split again
:笛声音;
kill
split again
:スクウォーク;
kill
endsplit
else (低ランク行列)
:連続性ラ音を除外した行列;
if (フーリエ・ウェーブレット分離[時間波形を見る]fistaとか) then (スパース行列)
:断続性ラ音;
split
:水泡音;
kill
split again
:捻髪音;
kill
endsplit
else (低ランク行列)
:呼吸音;
kill

単なる呼吸音だとしても呼吸の始まり・立ち上がりに音が強く録音されていると、ウェーブレット成分に誤検出することが多い。病気によって聴診箇所が変わるので呼吸音とラ音の強弱が変化してしまう。それによって検出結果が変わることも考慮する必要がある。

また、このような解析の手順では断続性ラ音を解析する前に連続性ラ音を解析してしまうと、2段階の推論という構造になるため、断続性ラ音の解析結果に信憑性が欠けてしまうことも問題点の一つである。

## 呼吸音(respiratory sounds)
主に呼吸運動によって生じ、胸部より聴診し得る音の総称。
気道・肺胞を換気する気流の音としての呼吸音(breath sounds)と、異常音である副雑音(adventitious sounds)に分類される。
特に副雑音はラ音とその他の副雑音に分類され、中でもラ音は連続性のものと断続性のものに分類される。

### 呼吸音(breath sounds)
---
#### 肺胞呼吸音(vesicular breath sounds)
正常な抹消肺に接する大部分の胸壁上、特に肺底部でよく聴取される。
音は柔らかく、吸気にはほぼ一定の大きさで聞かれるが、呼気では弱く聴取される。
肺胞呼吸音は肺胞換気の指標とされる。
#### 気管支呼吸音(bronchial breath sounds)
正常では胸骨上部の狭い範囲でのみ聞かれる。
音が大きく、好調な成分を持つ呼吸音で、吸気より呼気で大きく、吸気と呼気の間に明らかなポーズがあるように聴取される。
肺胞呼吸音が聞かれるべき部位で聴取されると肺炎や無気肺を疑う。
#### 気管支肺胞呼吸音(bronchovesicular breath sounds)
胸骨周囲、肺尖部、肩甲間部にて聴取される。
肺胞呼吸音と気管支呼吸音の中間と言われるが、その定義はやや曖昧らしい。
吸気・呼気ともに肺胞呼吸音より明瞭に聴取でき、吸気と呼気の間にポーズはない。
肺胞呼吸音が聞かれるべき部位で聴取されると病的と考えられる。
#### 気管呼吸音(trasheal breath sounds)
顎部気管上で聴取される(つまり肺とはあまり関係ないと思われる)。
音が大きく、粗い感じの呼吸音であり、吸気よりも呼気における音が大きく聴取され、吸気と呼気の間に明らかなポーズがある。

### 副雑音(adventitious sounds)
---
#### 連続性ラ音(continuous pulmonary adventitious sounds)

250[ms]以上持続するラ音。声や管楽器のように、安定して繰り返される振動波形をもつ音である。異なった周波数の音が混在した多音声の場合が多く見られ、類鼾音・笛声音、スクウォークに分類される。気流と気道壁と何らかの相互作用で振動が発生し、それが連続性ラ音として聴取される。比較的大きな音で呼気に多く聴取される。

* **```類鼾音(rhonchus)```**
周波数200[Hz]以下の低音性連続性ラ音。
* **```笛声音(wheeze)```**
周波数400[Hz]以上の高音性連続性ラ音。
* **```スクウォーク```**
吸気時に聞かれる、100[ms]前後の連続性ラ音。

**```連続性ラ音の特徴抽出```**
連続性ラ音は「基本周波数とその倍音のみから構成される音」という特徴を持っており、その特徴は時間波形よりもスペクトログラムで表現すると非常に特徴的な形で現れる。
したがって、連続性ラ音はスペクトログラムから特徴抽出することが適切である。
肺聴診音の信号を短時間フーリエ変換し変換後の信号をロバスト主成分分析によって連続性ラ音であるスパース行列とそれ以外の音として低ランク行列に分離することができる。
あくまで、連続性ラ音のスペクトログラムがわかっているわけではなく、特徴をもとにRPCAすると連続性ラ音のようなスペクトログラムが得られるという仕組みである。

#### 断続性ラ音(crackle)
持続性の短い不連続なラ音。
何かが破裂・爆発するような非楽音様の音である。
呼気時に閉塞した細い気道が吸気により再開放することが音源となるが、断続性ラ音は呼気にも吸気にもよく聞かれるものがある。
クラックル(crackle)とも呼ばれパルス状の波形として現れる。
スペクトログラムでは広い周波数帯の縦線上の瞬時周波数成分として現れる。

* **```水泡音```**
ブツブツと粗い感じで、大きな周波数の低い断続性ラ音。
細かい断続性ラ音に比べて持続時間は長く、10[ms]あるいはそれ以上であることが多い。
多くは各呼吸ごとでラ音の個数、出現タイミングに再現性がある。
* **```捻髪音```**
細かく、音が小さく、周波数が高い特徴があり、呼気に聴くことが多く、持続時間は5[ms]以下。
パルス状の波形が群発するという性質から周期、発生数、パワー、出現位置、振幅から定義される。
まれに、健常人でも聴取されることがあるので注意。
中でも代表的な特徴量である「波形」「発生数」「周波数」に着目し、その定義と算出方法を明確にしたい。

**```断続性ラ音の特徴抽出```**
断続性ラ音の特徴量の多くは典型的なパルス状の波形を多く含み、少数のウェーブレット成分で特徴づけられる。断続性ラ音はスペクトログラムでは縦線上の瞬時周波数成分として現れるが、時間周期的な特徴2CDやゼロクロス、IDW、LDWといった時間波形でしか捉えられないような特徴を持つため、時間波形的な特徴を見て抽出する。

#### その他の副雑音
* **```胸膜摩擦音```**
吸気、呼気富に出現する断続的な異常音である。肺の内部から発生する断続性ラ音と鑑別が困難なことがあるが、胸膜摩擦音は吸気、呼気いずれにも音が均等に出現するという特徴がある。
* **```Hamman's sign```**
心収縮中期に聴取されるクリック音。
心臓の音。

## なんかわかっておきたいこと
### 数学的知識について
**```フーリエ変換```**
<strong>時間tの関数から各周波数ωの関数への変換</strong>(<a href="https://www.yukisako.xyz/entry/fourier-transform">参照</a>)
心電図や音声波形など、時間に伴う複雑な波形はすべて単純なsin波が重なって形成されているだけ。単純なsin波にそれぞれ分けて考えることができれば、微積分が可能になるうえ解析がしやすくなる。そのsin波の概形を知るために必要な角周波数を基準に変換するのがフーリエ変換ということになる。
>フーリエ変換するとなぜ$\sin$波の概形を知ることができるかというと、波形を下の式のように定義しているため、各周波数は概形を知るための必要不可欠なパラメータとなり得るから。

***
$$
r \sin \omega t
$$
***
振幅rと角周波数ωに注目してωに対する振幅rの値を離散的に取り出す。例えば5sin2tだったら振幅5、各周波数2と読み取ることができる。

角周波数 $\omega$ [rad/s]は単位円において１秒間になす角度。
***
$$
\theta[rad]=\omega[rad/s]・t[/s]
$$
***
$\theta$ [rad]の成り立ちは円周/半径（ $2πr/r=2π$ ）弧度法という。
この時点で以下のように表現が可能になる。
***
$$
f(t) = \sin\theta = \sin\omega t
$$
***
1秒間に $f$ 回振動させる波を考えると、1秒間に $f$ 回振動させる場合は1秒間に $2πf$ 進めば良いことになるので
***
$$
f(t) = \sin\theta = \sin\omega t = \sin2\pi ft
$$
***
と表現することができるようになった。


複雑な時間波形でも時間波形ごとに分けて考え、更にそれをスペクトルに変換することで音の特徴を掴むことができる便利なツールがフーリエ変換。

**```離散フーリエ変換```**
時間領域から周波数領域への変換によって無限の数列に変わることが多く、実際の信号データとして扱うことができなくなるので、離散フーリエ変換を使って、有限数列として変換することで実際のデータを解析できるようになる。

[**```短時間フーリエ変換```**](https://www.jstage.jst.go.jp/article/jasj/72/12/72_764/_pdf)
フーリエ変換は周期信号を仮定しているので、時々刻々と変化する音楽などの信号には不適。そうした場合に用いるのが短時間フーリエ変換(STFT)。

$x(t)$を短時間の時間区間に分割して、その周波数構造を分析するために $0 \geq t \geq N-1$ だけで0でない値を持ち、それ以外では0となるような関数 $w_a (t)$ を $x(t)$ に乗じることを考える。
このときに用いられる $w_a(t)$ は分析窓と呼ばれ、Hamming窓、Hanning窓、Blackman窓などが挙げられる。
短時間フーリエ変換では、
***
$$
x_m(t-mS) = w_a(t-mS)x(t)
$$
***
というように分析窓関数 $w_a(t-mS)$ を $mS$ だけシフトして乗じることで離散時間信号 $x(t)$ を切り出し、切り出した時間区間の先頭に、時間原点を移動した有限長信号を定義する。
* $m$ : 有限長信号のインデックス(番号)
* $S$ : 有限長信号のシフト量
* $N$ : 有限長信号の長さ

を使って新たに $n = t-mS$ とおいて、有限長信号に離散フーリエ変換(DFT)を適用することで、短時間フーリエ変換

***
$$
X(m,k) = \Sigma^{N-1}_{n=0} x_m(n) e^{-j2 \pi kn / N}
$$
***

が定義される。

有限長信号の長さ(フレーム長)を長くすれば時間分解能は低く、周波数分解能は高くなり、逆に短くすれば時間分解能は高く、周波数分解能は低くなる。

**```離散コサイン変換```**
離散フーリエ変換は実数に対して複素数を返す $\exp$ を使うため実数列に対しても複素数列で値を扱わなければならない。
そこで、離散コサイン変換を使って $\cos$ を使った実数列での取扱いやすいものにする。

**```ウェーブレット変換```**
周波数解析手法の一つで、基底関数「ウェーブレット関数」を用いる。
小さなパルス信号の波を拡大・縮小・平行移動して足し合わせることで与えられた入力波形を表現する。
**時間的変化の特徴**と**周波数成分の混じり方**を知るために用いられ，時間と周波数にかかわる信号情報を同時に抽出することができる。詳しいことは<a href="https://www.cqpub.co.jp/hanbai/books/30/30961/30961_9syo.pdf">この資料</a>を参照のこと。
>**成り立ち**
>人の声の性質（声色，声紋）をはじめとして，楽器の音色や画像，自然界のあらゆる信号はすべて固有の周波数成分の混じり方を有している．こういった周波数成分の混じり方を解析できるのが，フーリエ変換である。
>しかし，フーリエ変換では信号のもつ性質をすべて周波数情報に置き換えてしまうため，時間に関連する情報が失われてしまう．言い換えれば，時間とともに刻々と変動する信号であっても，フーリエ変換では周波数成分だけしか得られないので，大ざっぱな信号情報しか見いだせない。
>短時間の窓関数に当てはめて時間を考慮した短時間フーリエ変換が考案されたが、解析の時間軸が固定されていることがネックとなり、時間と周波数を同時に解析するには物足りなかった。
>そこでウェーブレット変換が考案された。周波数の低い成分に対しては“解析する時間幅を長くする”，高い成分に対しては“逆に短くする”というアイデアのもと、解析する時間幅を周波数に反比例させることで多くの信号解析が合理的に行えるようになった。

石油の埋蔵場所を特定するため，人工地震波信号の精密な時間（埋蔵場所の特定）と周波数（埋蔵資源の分別，石油かどうか）の情報を同時解析するために考案された歴史がある。変換画像の効率的な符号化方式が開発されたため画像圧縮方式であるJPEG 2000に利用されるようになった。

**```離散ウェーブレット変換```**
離散ウェーブレット変換は、元信号を高周波成分と低周波成分に分解し、分解された低周波成分をまた高周波成分と低周波成分に分解するという処理を繰り返し行うことと等価である。多重解像度解析とも呼ばれる。

**```特異値分解（Singular Value Decomposition）```**
任意の行列 $\mathbf A$ (m×n)に対して $\mathbf A = \mathbf U \bf{\Sigma} \mathbf V$ と分解すること。
$\mathbf U$ : (m×m)の直交行列
$\mathbf V$ : (n × n)の直交行列
$\bf{\Sigma}$ : (m× n)の「非対角成分は0，対角成分は非負で大きさの順に並んだ行列」
この $\bf{\Sigma}$ の対角成分がそのまま大きい順に特異値となる。

例えば、顧客と商品の関係を2部グラフ(bigraph)にしたものを想定する。
<img src="https://i.imgur.com/qreYibO.png" width="600" alt="凸集合"><img src="https://i.imgur.com/QUJnrYC.png" width="600">
この顧客を $\mathbf U$ 、商品を $\mathbf V$ としたとき、「ある顧客は特定の商品が好きである商品が嫌い」や「この商品はこういう傾向のある顧客に好かれやすい」といった見方が可能である。
そしてこの傾向はそれぞれが特別な特徴を持っているわけではなく、似通った傾向のものが多い。同じような傾向は同じ特徴としてまとめるのが当然だろう。ということでまとめると、各々がお互いで表現できないような特徴(線形独立)を持った少数の特徴にまとめられる。この特徴をそれぞれベクトルで表すと、そのベクトルは互いに線形独立なベクトルであり、線形独立なベクトルの数だけ特異値の数が存在する。

**```特異値```**
数学的にいうと、行列 $\bf A$ とその転置行列 $\bf A^{T}$ の積の固有値の非負の平方根のこと。行列版固有値。
上の例でいうと、「商品が顧客に好まれる傾向」を表す。
これは、音の例だが線形独立な複数の行列でデータが表せるとして、その $\kappa$ 倍ずつされているこの値が特異値と呼ばれる。

**```直交行列```**
各成分が実数の正方行列、各列ベクトルは互いに垂直でその各ベクトルの大きさは1である。

**```固有値```**
ベクトルが線形変換された際に、ベクトルの方向は変わらずにスカラー倍されているだけのベクトルを固有ベクトルといい、スカラー倍の値のことを固有値という。まあ、つまり他の何からも影響されない独自のベクトルがあってその方向への引き伸ばし率という解釈。
そのベクトルがある点に向かって近づいていくと、どのような方向へ修正されていくのかということを表したベクトルが固有ベクトル、その値が固有値。**(怪しい。後で修正するかも)**

**```転置行列```**
行と列を入れ替えた行列。一般に $\bf A$ と転置行列 $\bf A^{T}$ の積をとるということは、行と列がひっくり返っているので自分自身との積和を求めていること(つまり内積を求めていること)になる。自身の転置行列との積を取ると、それぞれのベクトルが方向性を失ったスカラー量となるため大きさを表す。

**```内積```**
そもそも内積ってなんだったっけ？という人(私)のためにイメージしやすいものを。
内積とは一方がもう一方に対して行った仕事・貢献度のようなもの。垂直だったらもう片方に対して何もしていない同然、無駄ってことですね。

**```低ランク性```**
行列の中の線形独立なベクトルの本数が少ない行列のこと。

**```スパース性```**
ベクトルが少数の非ゼロ成分で構成されるという性質。

[**```最適(小)化問題```**](https://tech.curama.jp/entry/2018/04/20/120000")
最適化問題とは特定の集合上で定義された関数について値が最大・最小となる状態を解析する問題。最適化問題に関して、局所的には最小となる値を見つけられたとしても、それが大域的には最適な最小となるのは凸集合でない限り保証できない。

**```凸集合```**
ユークリッド空間における物体が凸(convex)であるということは、その集合に含まれる任意の２点を結ぶ線分がその集合の中に含まれる状態。

<img src="https://i.imgur.com/2KP44QL.png" width="200" alt="凸集合"><img src="https://i.imgur.com/EW6ym4I.png" width="200" alt="非凸集合">

**```凸最適(小)化問題```**
凸関数の最小化問題。局所的に求めた最小値が大域的な最小値と一致するという性質を持つ。凸集合を関数にした場合、複数の極小値は持たず一つの極小値が最小値と一致する。

[**```線型計画法```**](https://ja.wikipedia.org/wiki/%E7%B7%9A%E5%9E%8B%E8%A8%88%E7%94%BB%E6%B3%95)
いくつかの一次不等式および一次等式を満たす変数において、ある一次式を最大化または最小化する問題。
### グラフについて
**```波形信号```**
横軸に時間を、縦軸に音の大きさを示すグラフ。

**```パワースペクトル```**
ある時間の波形をフーリエ解析して横軸に周波数を縦軸に音のパワーを示すグラフ。

**```サウンドスペクトログラム```**
その時点 $t$ または $\omega$ での周波数 $f(t)$ または $f(\omega)$ を表すグラフ。スペクトログラムが大きいということは高周波数成分が存在するということである。時間信号を短時間フーリエ変換することで得られる。また、逆フーリエ変換によって時間信号に戻すこともできる。
音の時間的な変化、音色、高さ、大きさを同時に読み取ることができる。これならたしかにスペクトログラムにしたほうが解析が用意になる。

[**```標準偏差```**](https://blog.apar.jp/data-analysis/3390/)
データの、平均値からの散らばり具合を指す指標。
データが平均値の周辺にどのくらいの広がりや散らばりを持っているか。
標準偏差は、平均値と合わせて見ることによって、データを正しく把握することができる。

### 最適化問題
広く見ると、日常生活のあらゆる場面に最適化が浸透している。
災害時の復旧スケジュール計画、公共施設の配置計画、エネルギー需要バランスの制御など社会的な課題での最適化の重要性が高く、その最適な解決手段を取るための考え方。
* **数理モデル** : 多岐にわたる分野の課題を対象とし、対象となる問題を数式で記述したもの。
* **数理最適化** : 数理モデルに対して数理的な計算手法で最善策を求めること。
* **定式化** : 数理モデルを作成すること。
* **目的関数** : 数理最適化問題には目標や目的が存在し、その目的をどれだけ達成しているかという指標のこと。
* **制約関数** : 数理最適化問題の制約(条件)を表す関数。subject to。
* **最適解** : 最適化問題の解のこと。「与えられた制約条件をすべて満たし、目的関数 $f(x)$ の値が最小あるいは最大になるような $x$ の値」
* **緩和問題** : 実時間を最適化問題に定式化しても問題の複雑さや規模によって、そのまま実時間で解けないことがしばしば起こる。その場合に、与えられた最適化問題を比較的簡単な問題に緩めて解くことで元の問題の最適な目的関数値を見積もる手がかりを得るアプローチ。
* **連続緩和** : 組合せ最適化問題の各変数の整数条件を外して緩和する。
* **ラグランジュ緩和** : 制約条件に対する違反の度合いをペナルティとして目的関数に組み入れて制約条件を取り除く緩和法。

**BP**
BP(基底追跡)は $l_0$ ノルムを用いた最適化問題である貪欲法を $l_1$ ノルムを用いた最適化に近似したものである。
***
$$
\min_{\mathbf x}||\mathbf{x}||_1 \quad subject \: to \quad \mathbf{A}\mathbf{x} = \mathbf{b}
$$
***
**BPDN**
BPはBPDN(基底追跡ノイズ除去)の特別な場合である。BPDNとLASSOはLASSO回帰と等価である。ただし、 $\delta \geqq 0$ である。
***
$$
\min_{\mathbf x}||\mathbf{x}||_1 \quad subject \: to \quad ||\mathbf{b} - \mathbf{A}\mathbf{x}||_2 \leqq \delta
$$
***
**LASSO($l_2$最小化問題)**
$l_1$ ノルムに上界を設けた $l_2$ 最小化問題。定数 $r > 0$ である。

***
$$
\min_x ||\mathbf{b}-\mathbf{Ax}||_2 \quad subject \: to \quad ||\mathbf{x}||_1 \leq r
$$
***

**LASSO回帰($l_1$正則化最小二乗問題)**
LASSOはLASSO回帰の式に置き換えることができる。解は一致する。

***
$$
\min_x \frac{1}{2}||\mathbf{b}-\mathbf{Ax}||^2_2 + \lambda||\mathbf{x}||_1
$$
***

>最小化問題に対して公式を作ることはできないが、スパース解が組み合わせ問題であるため凸最適化問題の解としてスパース解を数値的に求める方法をスパース解法と呼ぶ。

### アルゴリズム
この最小問題を解くアルゴリズムとして以下のような手法(凸緩和)がある。
任意の凸集合の局所的最適解は大域的最適解であるため凸最適化問題に関して局所的最適解を求めるということは大域的最適解を求めることと等価である。

**ISTA**
LASSO回帰の最適解を求めるスパース解法はISTA(Iterative Soft shirnkage-thresholding)と呼ばれる、閾値処理のアルゴリズムである。

***
$$
\mathbf{x}^{(k+1)} := soft \left( \mathbf{x}^{(k)} + \frac{1}{L} \mathbf{A}^{\mathrm{T}} \left( \mathbf{b} - \mathbf{A} \mathbf{x}^{(k)} \right),\ \frac{\lambda}{L} \right)
$$
***

>ISTAの$l_1$ノルム、$l_2$ノルムのバランスを取るパラメータが$\lambda$なのだが、肺音に当てはめてみると、少し違和感がある。例えば、非ゼロ成分が少ないスパース成分を求める際に$l_1$ノルムはある程度小さいので、代わりに$l_2$ノルムが大きくなる(つまり忠実に再現しなくなる)。逆も然り。この２つのバランスを保っている$\lambda$を使ったISTAは実は肺音解析に向いていないのではという可能性があることが分かった。

**交互方向乗数法(ADMM)**
線形正則化項をもち、線形和最適化問題を解く反復アルゴリズムである。線形和の最適化問題ということなので、複数の最適解を求める問題ということ。
***
$$
\mathbf{x}^{(\mathbf{k}+1)} := \arg \min_x f(\mathbf{x})+\frac{\rho}{2}||\mathbf{y}^{\mathbf{(k)}}+\mathbf{F}\mathbf{x}+\mathbf{G}\mathbf{z}^{\mathbf{(k)}}-\mathbf{c}||^{2}_{2}\\
\mathbf{z}^{(\mathbf{k}+1)} := \arg \min_z g(\mathbf{z})+\frac{\rho}{2}||\mathbf{y}^{\mathbf{(k)}}+\mathbf{F}\mathbf{x}^{(\mathbf{k}+1)}+\mathbf{G}\mathbf{z}-\mathbf{c}||^{2}_{2}\\
\mathbf{y}^{(\mathbf{k}+1)} := \mathbf{y}^{(\mathbf{k})}+\mathbf{F}\mathbf{x}^{(\mathbf{k}+1)}+\mathbf{G}\mathbf{z}-\mathbf{c}
$$
***
凸関数 $f(\mathbf{x})$ と定数$\rho>0$が与えられたとき、点 $\mathbf{v}$ からの $l_2$ 正則化最小化問題の最適解を求める関数
***
$$
prox_{f,\rho}(\mathbf{v}) := \arg \min_x f(\mathbf{x})+\frac{\rho}{2}||\mathbf{x}-\mathbf{v}||^{2}_{2}
$$
***
を**近接写像**(proximity operator)という。

どんな最適化問題を解くときに使うかというと、
***
$$
\min_{(\mathbf{x},\mathbf{z})}f(\mathbf{x})+g(\mathbf{z}) \quad subject \: to \quad \mathbf{F}\mathbf{x}+\mathbf{G}\mathbf{z}=\mathbf{c}
$$
***
BPをADMMから反復計算のアルゴリズムを導出すると以下のようになる。
***
$$
\min_{(\mathbf{x},\mathbf{z})}f(\mathbf{x})+||\mathbf{z}||_1 \quad subject \: to \quad \mathbf{x}=\mathbf{z}
$$
***
### 肺音について
**```サンプリングレート```**
1秒間に標本化するデータの数。単位で表すと[データの個数/s]。これによって時間とデータ点の関係を知ることができる。グラフ化するときに便利。

## どうやって肺音を解析するのか
スパースモデリングは少数の非ゼロ成分からデータを表現するため、ノイズに対して頑健な手法であり、微弱な生体信号の解析に適している。定常的な生体現象(呼吸、心拍、etc...)を観測した信号ならば、類似した特徴が反復して出現したり同じ構成要素が含まれていると考えられる。この性質を線形従属性(低ランク性)として捉える。しかし、生体信号には例外的な構成要素(病異常信号、外的ノイズ、etc...)も含まれることがある。この例外的な信号を除外することができれば、線形従属な信号を取り出し、それ以外の信号と分離することが可能である。このように信号分離するためにロバスト主成分分析(RPCA)を用いて算出する。

上図のグラフにも少し書いているとおり、RPCAを使って聴診音から連続性ラ音のスペクトログラムとそれ以外の行列に分離できる。このとき信号のスペクトログラムを見て分離している。分離後のそれ以外の行列に対して、フーリエ・ウェーブレット成分に分離する処理を施すことで断続性ラ音を呼吸音から抽出している。このとき信号の時間波形を見て分離する。連続性・断続性それぞれの特徴を捉えやすい状態で解析している。

念頭に置いておきたいのは、どのような病気で肺のどの部分が悪いかという条件も考えて、
***
「結局、連続性・断続性ラ音はそれぞれどのような方法で特徴抽出できるのか。」そういうことを聞きたいのです。私は。
***

## 連続性ラ音の特徴抽出
肺音のスペクトログラムを行列 $\bf D$ と見なしたとき、呼吸音と断続性ラ音の混合音には「低ランク性」、連続性ラ音には「スパース性」の性質がある。つまり、呼吸音・断続性ラ音・連続性ラ音からなる肺音は低ランク行列 $\bf L$ とスパース行列 $\bf S$ の和と見なすことができる。実際の聴診音はノイズが入るのでノイズ $\bf E$ も考慮する。基本的にスペクトログラムに変換しないと解析ができない。

以下の式のように凸最適化問題として考える。
***
$$
\min_{\mathbf{L,S,E}}||\mathbf{L}||_*+\lambda_1||\mathbf{S}||_1+\lambda_2||\mathbf{E}||^2_F \quad subject \: to \quad \mathbf{D} = \mathbf{L} + \mathbf{S} + \mathbf{E}
$$
***

まず初めに理解しておきたいのは、スペクトログラム行列を低ランク行列とスパース行列に分離することで肺聴診音から連続性ラ音を取り出すことができるということである。肺音のスペクトログラムの行列 $\bf D$ を低ランク行列 $\bf L$ ,スパース行列 $\bf S$ ,ノイズ行列 $\bf E$ に分離することを考えると,凸最適化問題として定式化できる。

低ランク性とスパース性とノイズのバランスを取るパラメータ $\lambda 1$ と $\lambda 2$ があり、 $\lambda 2$ を $\infty$ として考えるとノイズ行列は $\bf E = 0$ と見なせる。よって、ノイズ正則化項がない最適化問題（ロバスト主成分分析の凸最小化問題）となる。
つまり、肺音スペクトログラムを分離する上式の最適化問題は <u>ロバスト主成分分析</u> (下式)である。
***
$$
\min_{\mathbf{L,S}}||\mathbf{L}||_*+\lambda_1||\mathbf{S}||_1 \quad subject \: to \quad \mathbf{D} = \mathbf{L} + \mathbf{S}
$$
***

### [**ロバスト主成分分析(RPCA)**](https://www.airc.aist.go.jp/achievements/ja/p-011.html)
与えられた行列を低ランク行列とスパース行列の和の形になるように分解する。大きな外れ値に対して頑健な低ランク行列近似法。もともと低ランク行列でありながら外れ値によって低ランク性を失ってしまった行列を低ランク行列と外れ値行列に分離することができる。ロバスト主成分分析の計算はバッチ算法である。

>ロバスト主成分分析の求解について**交互方向乗数法**や**閾値処理に基づく手法**による算法が実用的である。

### 西辻皐勢の研究[2018]
時刻 $t_j$ におけるスペクトログラム $\mathbf D$ の列ベクトル $\mathbf d^{(j)}$ が与えられたときに $\mathbf{d}^{(j)} = \mathbf{l}^{(j)} + \mathbf{s}^{(j)}$ のように、それまでの時刻の線形従属成分と線形従属になるような列ベクトル $\mathbf{l}^{(j)}$ とスパース成分の列ベクトル $\mathbf{s}^{(j)}$ を算出する。ただし,それまでの時刻の線形従属成分の上位 $r$ 本の主成分を列に持つ基底行列 $\mathbf{U}^{(j−1)}$ と,その特異 $\mathbf{\kappa}^{(j−1)}$ が与えられているものとする
>ロバスト主成分分析はバッチ算法であるため、肺音のデータ聴取と分離した音を聴くことを同時に実行できない。肺音の分離はこれを可能にするオンライン算法であることが望ましい。これをつかって、肺聴診音から連続性ラ音を分離する。

$\mathbf{l}^{(j)}$ は $\mathbf{U}^{(j)}$ に線形従属である。 $\mathbf{U}^{(j)}$ が $\mathbf{U}^{(j-1)}$ で近似できるならば、 $\mathbf{l}^{(-1)}$ によって $\mathbf{l}^{(j)}$ を近似できる。

このような問題は再構成問題と呼ばれ、再構成した際の誤差を最小化する問題となる。 $\mathbf{s}$ のスパース正則化のもとで再構成の $l_2$ 誤差を最小化する問題は、凸関数でスパース解を推定する最小化問題LASSO回帰と等価である。

>最終的にこのオンライン算法(オンラインRPCA)の中でLASSO回帰をISTAで実行しているという形になる。

### 尾道拓海の研究[2022]
**```U-Netを使った特徴抽出```**
先行研究であるRPCAによる連続性ラ音の分離では、肺音データに対して毎回反復計算を行う必要がある。また、その計算は低ランク性とスパース性にのみに基づいて分離される。しかし、実際の聴診音には患者による個人差や多数の聴診音に共通する特徴が含まれると考えられる。そこで、深層学習による肺音の学習を提案する。深層学習により多数の聴診音を学習することで、各肺音に共通する特徴や、呼吸器や身体と相関する特徴の学習が可能であると考えられる。しかし、聴診音は各肺音が混合したものが聴取されるため、正解データを準備することは困難である。そこで、肺音の特徴を表す事前知識として、RPCAの最小化問題を目的関数として与えることで教師なし学習が可能になる。
ネットワークに与えるRPCAの目的関数を以下の式に示す。
***
$$
\min_{\Theta}||\mathbf{D-S_{\Theta}(D)}\|_*+\lambda\|\mathbf{S_{\Theta}(D)} \|_1
\quad subject \: to \quad
\mathbf{D = L + S}
$$
***
${\Theta}$ はネットワーク内の学習可能なパラメータの集合を指す。

深層学習を行うDNNモデルには、入出力が肺音の信号であることと中間層で特徴の学習が可能であることからend-to-endでの学習が可能なDual Frame U-Netを用いる。[Dual Frame U-Net](https://arxiv.org/pdf/1708.08333.pdf)とは、畳み込み層による対照的なエンコーダ・デコーダ型の構造を持つ。通常のU-Netの構造に加え、デコーダでアップサンプリングする直前の特徴量から、エンコーダにおけるプーリング直後の特徴量を差し引く。通常、プーリングでは周波数特性の微細な構造が失われ、粗い構造がU-Netの深い層まで特徴抽出に使われる。Dual Frame U-Netはこの不平等を避ける。したがって、連続性ラ音の分離においては、微細な倍音構造を持つ連続性ラ音の分離に適している。
現状では、少数のデータの学習でも汎化性を確認できている。


## crackle(断続性ラ音)の特徴抽出
断続性ラ音は局所的に発生し、個々の発生区間がまばらであるため波形の最大ピークを機械的に検出したり発生数を機械的に計数することが難しい。呼吸音の発生と重なる場合にピークを持たなかったり、確実に断続性ラ音のみを分離できなければならないという条件が検出を難しくしている。

### 波形の算出
>上記のような課題を解決して断続性ラ音の特徴抽出するために振幅が非ゼロ地点を開始、ゼロの地点を終了として断続性ラ音の発生区間を特定する。

1. 分離したウェーブレット成分の波形に対し、ソフト閾値処理を使って無音区間を作ることで断続性ラ音の発生区間を特定する。
2. 特定した発生区間の波形に対して最大の絶対値を持つピークを最大ピークと呼ぶ。
3. 最大ピークの直前直後のゼロクロスを特定し、それからLDWを算出する。

**ゼロクロス**
信号の値が0を通過する事。音の信号は通常0を中心としてプラス側とマイナス側に振動しているため、符号が変わる際に0を通過する瞬間がある。この値が0を通過する事を「ゼロクロス」、通過する瞬間の事を「ゼロクロス点」と呼ぶ。

### 発生数の算出
>波形を算出したときと同様に、振幅が非ゼロ地点を開始、ゼロの地点を終了として断続性ラ音の発生区間を特定する。

1. 最大ピーク、2CDを算出。
2. 個々の断続性ラ音の区間算出。
3. 捻髪音/水泡音の判定。
$2CD < 10ms$ の最大ピークの総数から捻髪音の発生数を算出
$2CD \geq 10ms$ の最大ピークの総数から水泡音の発生数を算出

<参考>
crackleの定義は論文によって異なっているためどのような線引きをすればよいのかも曖昧である。

### 周波数の算出
>呼気吸気の区間が明確である聴診音に対し、低周波成分である心音をカットするために$60Hz$のハイパスフィルタを前処理としてかけ、以下のとおり算出する。

1. 短時間フーリエ変換により聴診音の各区間のパワースペクトルを算出。
2. 各区間のパワースペクトルの総和を取り、それから平均値を算出。
3. パワー最大の周波数を $F_{max}$ とし、周波数方向に積分して総面積 $25\%$ 、 $50\%$ 、 $75\%$ のところをそれぞれ $F25$ 、 $F50$ 、 $F75$ とする。

[**パワースペクトル**](https://www.onosokki.co.jp/HP-WK/c_support/faq/fft_common/fft_spectrum_6.htm)
対象とする波形データに対して、データがその周波数をどれだけ含んでいるか（パワー）を周波数毎に表したもの。

### タイミングの算出
>現在、吸気区間を判定する方法がないため、吸気区間が既知であることを前提に算出する。

1. 1吸気を3分割し、分割された区間を早い順にEarly、Mid、Lateとラベル付けする。
2. 断続性ラ音の発生数を計数。
3. Early、Mid、Late区間内の発生数をタイミングの個数として計数。

### 振幅の算出
>分離したウェーブレット成分に対して算出する。

1. 最大ピーク、2CDを算出。
2. 個々の断続性ラ音の区間算出。
3. 最大振幅と、そこからゼロクロスをまたいで一つ前のピークの振幅との絶対値和を取ることで $A_1$ 、 $A_2$ 、 $A_3$ を求められる（ $A$ はおそらく振幅(Amplitude)のこと。どういう定義か書いていないのでいまいちわかんないけど）
4. $\frac{A_2}{A_1}$ 、 $\frac{A_3}{A_1}$ を求めることで振幅がわかる。

>なぜ、これらの肺音特徴量が必要なのか。ということをはっきりさせたい。

肺音の特徴を算出するには様々な特徴量を求める必要がある。
波形・発生数・周波数・タイミング・振幅の5つを算出するために肺音特徴を知っておく必要がある。
**2CD**
two-cycle deflection。音が始まったとされる時点から2周期分の幅のこと。
**IDW**
initial deflection width。初めの半周期の幅を表す。
**LDW**
largest deflection width。最も幅が大きいとされる半周期の幅のこと。


これらを求めることでそれぞれの肺音特徴が算出されることになる。
### crackleの解析研究
肺音解析の研究は、肺音に明確な定義が存在しないことから、何が正解なのかを可視化したデータというものは存在しない。現段階では、同じ臨床データに対して、「医師の診断」と「機械的な解析結果」のすり合わせをしながら検証しなければならない。
しかし、「機械的に解析する」ことは人間の目や耳では捉えられないような異常音や、病状に関係するような特徴波形を発見する可能性がある。
医師の診断を基盤とし、解析によってより詳しい肺音の特徴を発見できれば人間が見つけることのできる特徴以上の医学的情報を提供できるのでは？
#### 西辻皐勢の研究[2018]
crackle(断続性ラ音)は典型的なパルス状の波形を持つ音である。
そのパルス状の波形であるという事前知識を用いて、聴診音をパルス状の波形を持つウェーブレット成分とフーリエ成分(厳密には離散コサイン成分)に分類し、そのウェーブレット成分からcrackleを見つけ出そうという手続き。

パルス信号は少数の非ゼロ成分から構成されているため、FISTAという最適化問題を解くアルゴリズムで分離している。スパースモデリングを使った手法となっている。

#### 山田頼弥の研究[2021]
**```unrolled ISTAを使った特徴抽出```**
スパースモデリングを使ったアルゴリズムのパラメータ・基底行列といった値の修正を施せるように辞書学習可能にしたものを提案。

pytorchを使った深層学習で肺音の解析をするために西辻さんのコードをpythonに書き直し、学習できるようにpytorch版に修正。
***
* スパースモデリングを用いて聴診音から呼吸音と断続性ラ音を抽出する技術がある[酒井ら, 2013]。

* スパース解と基底から非ゼロ成分がある程度小さくなるようにデータを表現する。

* この最小解を求めるための反復アルゴリズムがある。

* 計算をブロックに分け、一つのブロックがノードとなる計算グラフになる。

***

#### 若杉勇弥の研究[2021]

## 【付録】結局どの問題をどう解けばよいのか(SPM基礎勉強会資料より)
#### $l_0$
解が極めてスパースであるならば,その非ゼロ成分を次々と特定する貪欲法がお勧めである。実
際、凸緩和の解法よりもずっと少ない手間数で解が得られることが多い。しかし、途中で台の選択に失敗するとそこで解への収束が破綻する。そのような失敗は基底$\mathbf{A}$の性質に依存する。
#### $l_1$
凸緩和の解法は、貪欲法よりも計算時間を要することが多いが、任意の初期値からスパース
解へ収束する反復計算になっているので、極端にスパースではないスパース解を求める応用問題に向いている。実際、画像や音はフーリエまたはウェーブレット基底で表現すると、ある程度のスパース性が現れる。また、貪欲法で失敗しやすい基底 $\mathbf{A}$ や低次元の $\mathbf{b}$ に対しても凸緩和の解法で適切なスパース解が得られることがある。スパース解が $\mathbf{Ax} = \mathbf{b}$ を厳密に満たすならば OMP や BP で良いが、ノイズの混入を考慮すると BPDN 、 LASSO 、 $l_1$ -LS などが実用的であろう。ただし、凸緩和では定数を利用者が設定する必要があり、しばしば手探りで見つける作業を要する.

## 特徴抽出が難しい事例
**断続性ラ音が独立せずに、連続して発生している聴診例。**
２つ以上の断続性ラ音が重なった部分は最大ピークの数が一つにまとめられるため、最大ピークの推定を誤り、波形や振幅の特徴を見誤ってしまう。
**実際の聴診音に対して断続性ラ音のパワーが小さいこと。**
断続性ラ音のパワーが強い聴診音を分離するとパルス状の成分が多く抽出されるが、パワーが弱い聴診音を分離するとウェーブレット成分に呼吸音の成分が誤抽出される。断続性ラ音の特徴を算出することが困難になる。
**断続性ラ音のみの信号を完全に抽出するのが困難であること。**
断続性ラ音のみの成分をウェーブレット成分として完全に抽出できるわけではない。