十字韌帶重建術後膝部異常聲音之訊號分析

(1)

行政院國家科學委員會專題研究計畫成果報告

十字韌帶重建術後膝部異常聲音之訊號分析

計畫類別：個別型計畫計畫編號： NSC91-2314-B-002-254- 執行期間： 91 年 08 月 01 日至 92 年 07 月 31 日執行單位：國立臺灣大學醫學院骨科計畫主持人：江清泉報告類型：精簡報告處理方式：本計畫可公開查詢

中華民國 92 年 10 月 16 日

(2)

前十字韌帶重建術後膝部異常振動信號之數學模式與頻譜分析

Mathematical Modeling and Spectrum Analysis of Abnormal Knee Clicks

in the Anterior Cruciate Ligament Reconstruction Patients

Ching-Chuan Jiang , Ju-Hong Lee , Yung-Chung Luo

摘要

在骨科的臨床上發現，一旦人體膝關節發生病變時，在運動時膝關節便會產生異常的聲音，亦即為膝關節在擺動下會產生振動信號。 Vibration Arthrometry（VAM）即是藉由分析此一振動信號來診斷膝關節之病變。由於 VAM 是一種非侵襲性的檢查工具，因此即發展潛力。本論文即是應用 VAM 在量測到之膝關節髕骨振動信號上，並利用自迴歸模型（Autoregressive Model）與倒頻譜模型（Cepstral Model）經由可適性區段（adaptive segmentation）後所建立之信號數學模式，與臨床發現對照後，找出膝部前十字韌帶（Anterior Cruciate Ligament；ACL）重建術後病變患者與退化性膝關節炎（degenerative osteoarthristis；OA）患者之特徵參數。在本研究完成後，希望能夠將所建立之專家診斷系統普遍地應用在骨科臨床診斷上，提供醫師們一具非侵襲性而準確的診斷工具，並將其納入全民健保給付的參考標準中，俾能在解除病患痛苦的同時避免醫療資源的浪費，以增進社會廣大民眾的福祉。

簡介

膝關節（knee joint）是人體中最大的活動關節，也是滑液性關節，亦是人體所有關節中最常受傷的關節[1][54]。因此前方膝痛（anterior knee pain）可說是骨科門診中最常見到的問題之一，而其發生之主要原因大多是與髕－股膝關（patellofemoral joint）有關。且根據統計，膝關節的外科手術占了全部關節外科手術極高的比例。隨著台灣社會型態趨於老年化，退化性關節炎（degenerative osteoarthritis；OA ）也愈來愈常見；一般，在 60 歲以上的老人，女性有 25％，男性有 15％會有厲害的膝關節退化性關節炎。其中更有相當比例的病患是兩腿都需要手術。除了年紀造成的退化性關節炎外，尚有少部分病患因為痛風，創傷等等的原因會造成膝關節的退化[53]。在臨床上骨科醫師發現一種很特殊的外源性膝關節振動信號：膝部前十字韌帶（Anterior Cruciate Ligament；ACL）重建術後，有些病人膝部會產生連續且綿密的振動信號，令病人相當困擾。膝部前十字韌帶受傷是年輕人最常見的運動傷害，發生的年齡大多在 40 歲以下。在台灣最常見的是打籃球上籃時突然感到膝關節有東西拉斷，接著出現膝關節腫脹疼痛；慢慢的，病人會發現不能跑、跳，做急轉彎時，膝關節好像會脫臼似的。如不處理，股四頭肌會逐漸萎縮，大腿腿

(3)

圍變小。對病人最大的困擾是無法再從事跑跳的運動。治療 ACL 斷裂，目前最有效的方法是用自體髕骨兩端帶骨韌帶行韌帶重建手術。術後效果很好，能使膝關節恢復穩定性，讓病人重拾運動的樂趣。但是有些病人術後會出現當膝關節活動時病人可以感到或觸覺到膝關節一連串綿密的聲音，讓病人困擾萬分。如能偵測此聲音加以分析診斷，便可成為骨科醫生臨床診斷時另一種重要的輔助工具。目前診斷膝關節病變的工具大致上可區分為三類：第一類為侵襲性（invasive）的檢查，例如關節鏡（arthroscopy）。關節鏡雖然能夠清楚地看出關節內部組織的狀況，但是以關節鏡來觀察關節軟骨的磨損情況，無論就時間或金錢上而言都是不經濟的，而且在手術進行前需先行麻醉，因此具有一定的危險性。第二類為非侵襲性（noninvasive）的檢查，例如；電腦斷層掃描（computer tomography

scanner；CT）、核磁共振造影（magnetic resonance imaging；MRI）、X-ray imaging

等。它們能診斷出較嚴重的關節軟骨磨損，但對於診斷初期的軟骨磨損並無法提供有效之幫助。最後一類為半侵襲性（semi-invasive）的檢查，例如關節腔攝影。此種技術與 X-ray 與 CT 一樣，病患需接受放射線照射，而且無法對其軟組織成像，因此對於診斷初期軟骨病變的功效亦是有限。相較於上述檢查膝關節病變之工具，利用關節病變所產生的聲音來診斷膝關節的狀況，將是一種具有非侵襲性、簡單而價廉的一種臨床檢查工具。若能利用此種方法，將各種膝關節病變所產生的振動信號聲音加以收集並分析，以建立專家診斷系統，將有助於提昇此方法的實用性。

研究動機

在骨科的臨床上發現，當膝關節發生病變時，在運動時會產生異常的聲音。 vibration arthrometry（VAM）是藉由分析膝關節的振動信號來診斷膝關節之病變，經由已發表的文獻報告可以發現 VAM 是一種極具發展潛力的非侵襲性診斷工具。本研究的主要目的，即是利用由膝關節髕骨表面所擷取到之振動信號，並應用信號處理方面的技術於所量測到的振動信號上，用數學模型來加以分析時域（time domain）及頻域（frequency domain）下之特徵參數，以找出足以區分出膝關節病變的特徵參數，進而找出其與膝關節疾病的關係，達到分辨病症以提供骨科醫師臨床診斷之參考。

本研究係將加速規置於髕骨外面中心點的皮膚表面，以量測膝關節在快速（67°/Sec）與慢速（2°/Sec）擺動下的振動信號[31-32]。在慢速擺動下的振動信號，即是所謂的髕股顫振（Physiological Patellofemoral Crepitus；PPC），至於快速擺動下的振動信號稱為 Vibration Arthrometry signal（VAM）信號。本研究將會利用線性預測（linear prediction）演算法的技術建立其數學模型，從時域（time

domain）及頻域（frequency domain）中，嘗試找出可以區分出：正常者（Normal）、

ACL 重建術後患者、退化性關節炎（OA）重建術前患者之特徵參數。本論文所使用的數學演算法為自迴歸模型（Autoregressive Model）及倒頻譜模型（Cepstral

(4)

Model）。

Ⅲ 材料與方法

人體之膝關節於正常情況下的伸直-彎曲擺動即會有振動信號之產生，而此一振動信號則可能會反映出膝關節內部組織的情況。一般而言，膝關節的病變可能會對膝關節之內部組織造成某種程度的影響，因此，可預知的是，從病患膝關節所量測到的振動信號其特性將會和正常人有所不同。也因為不同膝關節的病變對於膝關節內部的組織可能會造成某種程度之影響，因此所量測到的膝關節振動信號之特性也應該會不一樣。關節振動測量術（vibration arthrometry; VAM）[3] 基本上即是利用此一原理，藉由量測膝關節於彎曲擺動時所產生之振動信號，並應用信號處理的技術在所量測到的振動信號上，建立數學模式，以找出足以區分膝關節病變的特徵參數。

髕股顫振（physiological patellofemoral crepitus；PPC）是一種由正常膝關節

在慢速（小於 50/秒）主動擺動下所產生在髕骨之振動信號，它是由 Beverland 於 1986 年利用關節振動量測術（vibration arthrometry）所首要發現的[33]。且有研究人員量測了 250 個正常者的膝關節，並根據此結果統計出約有 99％的正常者的膝關節可以量測到 PPC 信號[33]。因此，一般咸認為 PPC 信號是由髕骨與股骨之關節軟骨在相對慢速運動下的黏滑（stick slip）現象所造成，而此一信號將可反映出關節軟骨之完好情況。曾有研究藉由活體實驗證明了 PPC 信號確實是由關節軟骨在滑動下之相互磨擦所產生[31]，該研究並指出 PPC 信號會受到關節軟骨的完整性、接觸面積與壓力、擺動速度、髕骨重量以及軟骨附近的軟組織等因素之影響。且依據由 PPC 信號在時域中所得之特徵參數，將退化性關節炎重建術前區分成三種型態，而此三種型態在骨科臨床上分別對應到三種不同之軟骨磨損情況。之後又有研究藉由將 PPC 信號數學模型化，而在頻域中找出可更精確地區分出三種退化性關節炎重建術前之型態的特徵參數[29-30][32-33]。此外，亦有文獻指出[26]，當相鄰之 PPC 信號的間隔之變異數變小時，則其頻譜之峰值會向高頻移動。然而該研究僅僅只對 PPC 信號的數學模型做模擬，並沒有實際量測取得人體膝關節之 PPC 信號作分析。本研究將對膝關節於快速擺動以及慢速擺動下所產生的振動信號加以研究分析，其中膝關節在快速擺動下所產生的振動信號稱之為 VAM 信號，而膝關節在慢速擺動（小於 5 /sec）下所產生的振動信號則被稱為髕股顫振（physiological patellofemoral crepitus；PPC）。為了使每一次所量測紀錄到的膝關節振動信號，其研究分析之結果能夠在一個公平合理的基礎上進行比較，因此，在進行膝關節振動信號量測實驗時，必須固定其擺動之速度。根據文獻上的記載[3]，67 /sec 此一速度最適合人體的膝關節作持續性的伸直-彎曲擺動，因此規定快速擺動的速度為 67 /sec 。至於慢速擺動的速度則是依循之前相關的研究為 2 /sec[28][31-33]。在本研究中，首先針對 VAM 信號進行量測，而後為 PPC 信號之量測。對於所量測得到之膝關節振動信號，本研究將運用線性預測（linear

(5)

prediction）之頻譜估側技術建立其數學模式，從頻域（frequency domain）中嘗試找出足以區分出正常者、前十字韌帶重建術後病變患者與退化性膝關節炎患者其相互群組間之特徵參數。由於自迴歸模型（autoregressive model;AR）係數可以得到整個信號的動態及頻譜訊息，因此具有相當高之解析度；而且自迴歸模型較其它時序（time-series）模型常用於隨機（stochastic）信號的模擬，其受歡迎的原因是在於係數之數學演算法是線性的，故它亦具有低運算複雜度之優點[37][39]，因此在本研究中將採用此一數學模型。自迴歸模型係數相當於數位語音處理中的線性預測係數，在語音辨識上，倒頻譜模型係數（Cepstral model coefficients）經證明比自迴歸模型係數效果好，並從文獻[41][50]中得知，倒頻譜模型係數在特徵空間（feature space）內比自迴歸模型係數更具較高之可區別性，且在做表面肌電圖型樣辨識（ surface electromyography pattern recognition）時，倒頻譜模型係數的平均辨識率至少高於自迴歸模型係數 5%，因此在本研究中亦將應用此一數學模式於膝關節振動信號之處理。

A 研究對象

本論文研究之對象為前十字韌帶重建術後病變患者以及退化性膝關節炎患者。在前十字韌帶重建術後病變患者方面，收集的病例都是重建術後膝部會產生連續且綿密的振動信號，令病人相當困擾之患者，此患者一共收集了 40 個病例，其中男性 23 例，女性 17 例，平均年齡 35.4 歲。在退化性膝關節炎患者方面，所收集的病例都是必須進行人工膝關節置換之病患，並且在術前量取其膝關節振動信號，此患者一共收集了 52 個病例，其中男性 6 例，女性 46 例，平均年齡 67.2 歲。此外，也收集了 40 例經過骨科醫師判定其膝關節為正常之自願者的膝關節振動信號以作為本研究之對照組，其中包括男性 26 例，女性 14 例，平均年齡 28.7 歲。

B 實驗設備

本研究所使用的儀器設備主要為一套資料擷取系統，我們稱之為關節振動量測儀。其組成的儀器包含加速規（ accelerometer ）、電子量角器

（electrogoniometer）、放大器（amplifier）、濾波器（filter）、類比數位轉換器（A/D

converter）以及個人電腦（personal computer）等，如圖 3.1。茲將各組成儀器的作用以及相關之規格說明如下： 1.加速規：這是本系統的感測元件(sensor)，負責將振動信號轉換成電信號後傳送給放大器。加速規可用來測量物體的加速度，因為一般表面振動的位移量都很小，不容易測量，所以大多以量測物體表面的加速度來代表物體的振動情形。壓電式加速規(piezoelectric accelerometer)是透過一壓電晶體，將振動信號轉

(6)

換成電信號，如圖 3.2，壓電晶體在受到一軸向力時會產生形變，造成晶體兩端的電壓差，相反的在晶體兩端加上不同的電壓差，將使晶體產生形變。本實驗所採用的加速規為一壓電式加速規（PCB PIEZOTRONICS, NEW YORK, U.S.A ； Model ： 352A22 ）。此型加速規所具有的頻寬最寬 (3dB 頻寬為 1Hz~16kHz)，且因其形狀扁平、體積小而適合固定於人體膝關節的表面。其尺寸大小為 3.3(mm) × 6.35(mm) × 9.14(mm)，重量為 0.48 公克，±3dB 頻寬為 1Hz~16kHz，靈敏度為 10mV/g。詳細的規格請參見附錄。 2.電子量角器：其主要元件為一圓形線性可變電阻，利用支架固定於膝關節外側轉軸處，當膝關節彎曲時會帶動可變電阻之懸鈕隨之轉動，如此隨著角度的改變，將使可變電阻產生相對應之電阻值，若在電阻兩端加上一定電壓源，則可獲得與角度變化成正比的電壓值。本實驗所採用之可變電阻其規格為 20k Ohm(Ω)，線性度為 1%，電壓與角度之關係為 25mv/degree。 3.放大器與濾波器：其作用為將加速規所感測到的振動信號加以放大及濾波以利區分及信號的檢視與紀錄。本實驗所採用的放大器為與加速規同廠牌的四通道放大器（Model：482A16），工作寬頻為 0.22Hz~100kHz，最大輸出為±10v，信號對雜訊比(SNR)大於 72dB，輸出阻抗小於 50Ω。本實驗所設定使用的增益

（gain）為 100，在此增益下的信雜比（signal to noise ratio, SNR）大於 72dB。除了上述的放大器以外，本實驗另外還在類比數位轉換器（A/D converter) 前加裝一濾波器，以避免高頻雜訊在類比數位轉換器作取樣（sample）時造成干

擾效應（aliasing effect）。本實驗所採用的濾波器為 AAF-3 抑制干擾（anti-alias）

濾波器（Costa Mesa, U.S.A），10kHz 的低通濾波器，此型濾波器能夠使用軟體程式來控制截止頻率。

4.類比數位轉換器：其作用乃是將所量測到的類比信號轉換成數位形式，使電腦可加以儲存與運算。本實驗係採用 Bakker BE490 類比數位轉換器（Bakker

electronic, Atlanta, Georgia3003, U.S.A），它是一個具有 8 通道（channel）的 A/D

card，可允許的最大取樣速率（sampling rate）為 1MHz，解析度（resolution）為 12bits/sample，並具有 256KB 的記憶體以及取樣並暫存（sample and hold）的功能，亦即對每個通道所取樣的時間都是同步的。在本實驗中我們僅使用 A/D converter 其中 4 個通道，其中第一個通道的信號是來自電子量角器，第二至第四個通道的信號則來自加速規，每個通道的取樣速率為 50kHz，經過轉換後的數位信號則直接以檔案的型式存入硬碟。其餘規格請參見附錄說明。

5.個人電腦：其功能是儲存所量測到的信號，並對信號作分析處理。本實驗係使用與 IBM（International Business Machine）相容的 Pentium-200 機型。

在本研究所使用的儀器除上述的關節振動量測儀外，本實驗在量取信號時尚需一個能夠輔助病人的膝關節做等速度擺動的儀器。我們所使用的是等速肌力測

(7)

試儀（Biodex isokinetic dynamometer, Biodex； New York, U.S.A），該儀器的旋轉臂可依照使用者所設定的角度範圍與擺動速度作等速度擺動。

C 實驗方法

首先讓病患坐於等速肌力測試儀（Biodex isokinetic dynamometer; Biodex,

New York, USA）上，此儀器的功用乃是使患者的膝關節可以以67o/sec的等角

速度來運動，為了量測 VAM 信號，將一個 PCB（PCB DIEZOTRONICS, New York, USA）加速規（model: 352A22）放置在膝關節內、外側脛骨粗隆(medial tibial bony protuberance)處的皮膚表面上，如圖 3.3，並以橡皮帶圈覆的方式固定加速規，

膝關節從 flextion( O

90 )到 extention(0 )再回到O 90 來回擺動，並由 personal O

computer（PC）進行記錄，在記錄之前，經由 PCB 加速規輸出的信號需先經 PCB （PCB DIEZOTRONICS, New York, USA）放大器（model: 482A16）將信號予以放大，之後需進行濾波，本實驗使用 AAF-3 anti-alias filter(Costa Mesa, USA)； 10kHz 的低通濾波器，濾波後的類比信號須經 analog to digital converter 的數位化處理後，信號才能進入 PC 中予以儲存，本實驗使用 Bakker BE490 A/D converter（Bakker electronic, Atlana, USA）對 PPC 信號作 50k Hz 的取樣速率及 12 bits/sample 的量化，電子量角器被置於膝關節的外側，藉由電壓的變化來表示角度的變化。茲將實驗步驟分述如下[31]： 1.記錄受測者基本資料。 2.受測者坐於 Biodex 的座椅上，此時雙腳懸空膝關節自然彎曲至 90 度。 3.將加速規和量角器固定妥當。 4.設定等速肌力測試儀旋轉臂的角度，由膝關節完全伸直到膝關節彎曲 90 度。 5.設定 Biodex 旋轉臂的移動速度為 67o/sec。 6.打開關節振動量測儀，並確定信號線接妥，儀器設定符合前節所述。 7.告知受測者，請其小腿前緣緊貼著旋轉臂的護墊，做等速的伸直彎曲運動。 8.啟動旋轉臂，在受測者完成數次膝關節彎曲-伸直-彎曲的動作、熟悉膝關節的運動速度後，才啟動電腦以取樣速率之 50kHz 來記錄一個膝關節運動週期的信號。在步驟 8 中，待受測者完成數次的彎曲-伸直-彎曲擺動動作後，才開始記錄振動信號，其目的是在於減少因受測者本身的肌肉收縮所造成的干擾信號[31]。在量測 PPC 信號時，除了將等速肌力測試儀之旋臂的擺動速度設定為 sec / 2o 的等角速度運動，及將類比數位轉換器之取樣頻率設定為 5kHz 外，其餘實驗方法完全相同。

D 數學演算法

由於膝關節振動信號的產生，被模擬成一種隨機過程，因此可用時序模型模擬之。自迴歸模型假設信號是線性過程產生，而且信號的二次動差（second

(8)

monent）是穩定的（stationary），雖然膝關節振動信號既非線性亦非穩定，但經過可適性區段 (adaptive segmentation)之方式處理後，卻可視為符合這些條件，而且自迴歸模型係數的演算法，較其他非線性模型簡單快速，因此常被用來模擬生物醫學信號。應用自迴歸模型轉換後可得倒頻譜模型，此兩種模型轉換後的係數，將用來當作區分膝關節病變之特徵參數。此外，經由模型化後之特徵參數，將可求得信號的功率頻譜密度（Power spectral density；PSD）以及在不同頻帶範圍內的信號

主極點功率比（Spectral power ratio of dominant pole）、最大峰值功率比（Spectral

power ratio of maximum peak）等之重要特性依據，最後由單方變異數分析（one-way-anova）的方法，來判斷最適合區別膝關節病變的特徵參數。除了利用這些特性資料來加以區別發生病變時之膝關節的振動信號特性外，這些特性資料未來亦可提供作為研判膝關節生理組織狀況的應用。 1.自迴歸模型（Autoregressive Model）頻譜估測技術中最常用到的時序建立模型之方法是自迴歸模型，因為自迴歸模型係數之精確估計值可以用解一組線性方程式求得。自迴歸模型是假設系統之轉移函數（transfer function）的極點（pole）都在單位圓（unit circle）內且轉移

函數之零點（zero）不存在，因此自迴歸模型又稱為全極點模型（all pole model）。

自迴歸模型基本上是一種線性預測（linear prediction），它是利用一組已量測到的資料值以線性組合方式來預測下一時間會量測到的資料值。假設

( )

n x ,x(n−1),x(n−2),....為靜止的（stationnary）隨機信號，且x

( )

n 的現在值為其過去值與某個輸入驅動信號u

( )

n 的線性組合，其表示式如下[38-39]：

( )

n a

( ) (

k x n k

) ( )

u n x P k P − + − =

∑

=1 (3.1）其中 a_P

( )

k ,k=1,2,3,...,P 為自迴歸模型的係數或預測的係數 (prediction coefficients)，P 為自迴歸模型的階數（order），u

( )

n 為輸入驅動信號。在一般的情況下（包括本研究）u

( )

n 是無法得知的，通常假設u

( )

n 為一白色雜訊(white

noise)，其平均值（mean）為零，功率頻譜密度（power spectral density）為一常數。目前常用來演算自迴歸模型係數的方法，一種是批次式（batch）的，就是將一段資料一起處理，這種方法在模型的階數未定時較適合，可用來嘗試不同階數的模型。另一種是連續式（sequential）的，先設定模型的階數，當每一筆資料進來後，都要重新計算自迴歸模型係數，這種方法適用於信號的頻譜會隨時間變化的情況。在本研究中之膝關節振動信號屬於連續式的，故預先設定階數為 25。圖 3.4 為自迴歸模型在時域（time domain）中的圖示。若將（3.1）式之等式兩邊

(9)

做 Z 轉換 ,則

∑

= − + = = _P k k P k z a z U z X z H 1 ) ( 1 1 ) ( ) ( ) ( (3.2) 此為自迴歸模型在頻域（frequency domain）中的表示法，如圖 3.5。而(3.1)式亦可視為一個線性濾波器（linear filter），且其系統之轉移函數為(3.2)式。在自迴歸模型中，假設輸入驅動信號u

( )

n 是未知的，例如在本研究中的膝關節振動信號即是。因此可將x

( )

n 只由過去值的線性組合來得到近似值，表示式如下：

( )

_∑

( ) (

)

=

−

=

P k P

k

x

n

k

a

n

x

1 ^ （3.3）在此選擇與自迴歸模型的階數相同。由（3.2）式、（3.3）式可得到介於實際值x

( )

n 與估計值

xˆ

( )

n

間的預測誤差(prediction error)為：

( ) ( ) ( ) ( )

∑

( ) (

)

=

−

+

=

−

=

P k P

k

x

n

k

a

n

x

n

x

n

x

n

e

1 ^ （3.4）欲求解自迴歸模型的係數：

a

_P

( )

k

，乃是基於下列之條件： minimize

{

( )

}

2

n

e

E

_p

=

_（3.5）

其中E 稱為均方預測誤差(mean squared prediction error)。由（3.5）式可得： _p

          ₊ ₋     ₊ ₋ =

∑

= = P P P k P P E x n a k x n k x n a x n E 1 * * * 1 ) ( ) ( ) ( ) ( ) ( ) ( l l l    ₊ ₋ ₊ ₋ =

∑

= = P P k P P x n x n a k x n k a n x n x n x E 1 1 * * * * ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( l l l     ∑ = ∑= − − + p k p n x k n x p a k p a 1 1 ) ( * ) ( ) ( * ) ( l l l

{

}

_∑

{

}

_∑

{

}

= = − + − + = P P k P P E x n x n a k E x n k x n a n x n x E 1 1 * * * * ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( l l l

{

}

∑∑

= = − − + P k P P P k a E x n k x n a 1 1 * * ) ( ) ( ) ( ) ( l l l

(10)

∑∑

∑

= = = = − + − + + = P k P xx P P P k xx P P xx P xx a r a k r k a k a r k r 1 1 * 1 1 * ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 0 ( l l l l l l （3.6）其中

{

}

{

(

)

(

)

}

)

(

)

(

)

(

* *

m

n

x

n

x

E

n

x

m

n

x

E

m

r

_xx

−

=

+

=

{

( ) ( )

}

) ( m E x n x* n m r_xx − = +

{

}

) ( ) ( ) ( * * * m r n x m n x E xx = + = （3.7）且*為取共軛數之意定義：

[

]

T p xx xx p r r p r ≡ (1) . . . . ( ) ₁_* （3.8）

[

]

T p P p P a a p a ≡ (1) . . . . ( ) ₁_* （3.9）因此（3.6）式又可以寫成： p p H p p H p p H p xx p r a r r a a R a E = (0)+ + + ₋₁ （3.10）其中 1 − p R =                 − + − + − − ) 0 ( . . . . ) 1 ( . . . . . . . . . . . . ) 2 ( . . . ) 0 ( ) 1 ( ) 1 ( . . . ) 1 ( ) 0 ( xx xx xx xx xx xx xx xx r p r p r r r p r r r （3.11）所以E 可以被表為_p a 之二次多項式，基於求解_p a 的條件，因此存在一組_p a_p 使得E 有極小值，也就是令_p p p a E ∂ ∂ =0 ，將會得到一組最佳解。其中_p

[

]

T P p 0 0 . . . 0 1* 0 = 因此將（3.10）式對a 微分，得 _p p p a E ∂ ∂ =rP +R_P−1a_P =0_P

(11)

則 a_p=−R_p−1₋₁r_p （3.12）將（3.12）式帶入（3.10）式中，可整理得 Ep =rxx(0)+rpHa_p （3.13）由（3.12）式及（3.13）式，可得到 Yule-Walker 方程式如下所示                     =                                         + − − − 0 . . . 0 ) ( . . ) 2 ( ) 1 ( 1 . ) 0 ( . . . . ) ( . . . . . . . . . . . . . . . . . . ) 1 ( . . . ) 0 ( ) 1 ( ) ( . . . ) 1 ( ) 0 ( _P P P P xx xx xx xx xx xx xx xx E p a a a r p r p r r r p r r r （3.14）或                     =               − 0 . . . 0 1 ) 0 ( 1 p P p p H p xx E a R r r r （3.15）

基於的運算效率，選擇以 Levinson Durbin 演算法來求解 Yule-Walker 方程式 [38][52]，可得自迴歸模型的係數。以上之求解過程在學理上稱為向前預測理論（forward prediction），而其之所以被稱為向前的原因在於目前的樣本值是由之前 p 個樣本值的線性組合而成。至於向後預測理論（backward prediction）的信號x

( )

n 以下式來表示[37-39] [52]：

∑

= + + − = P k b P k x n k u n a n x 1 ) ( ) ( ) ( ) ( （3.16）其中ab_P(k)為向後預測的係數，它與向前預測係數a_P(k)的關係為

a

_Pb

(

k

)

=

(

a

_P

(

k

)

*

k

=

1 L

P

（3.17）同樣地，與向前預測的 Yule-Walker 方程式的推導觀念相同，我們可以推得向後預測的 Yule-Walker 方程式如下表示：

(12)

                    =                     −                     + − − − b P b P b P xx xx xx xx xx xx xx xx E p a p a r p r p r r r p r r r . . . 0 0 1 . . . ) 1 ( ) ( . ) 0 ( . . . . ) ( . . . . . . . . . . . . . . . . . . ) 1 ( . . . ) 0 ( ) 1 ( ) ( . . . ) 1 ( ) 0 ( （3.18）將（3.14）式及（3.18）式分別改寫成：                       ∆ =                                             + − − − − − +1 0 . . . 0 0 ) ( . . ) 2 ( ) 1 ( 1 ) 0 ( . . . ) ( ) 1 ( ) 1 ( . . . . ) ( . . . . . . . . . . . . . . . . . . ) ( . . . ) 0 ( ) 1 ( ) 1 ( . . . ) 1 ( ) 0 ( P P P p p xx xx xx xx xx xx xx xx xx xx xx E p a a a r p r p r r p r p r r r p r r r （3.19）                      ∆ =                       −                       + − − − − − ₊ b P P b P b P xx xx xx xx xx xx xx xx xx xx xx E p a p a r p r p r r p r p r r r p r r r . . . 0 0 1 . . . ) 1 ( ) ( 0 ) 0 ( . . . ) ( ) 1 ( ) 1 ( . . . . ) ( . . . . . . . . . . . . . . . . . . ) ( . . . ) 0 ( ) 1 ( ) 1 ( . . . ) 1 ( ) 0 ( * 1 （3.20）其中，

∑

= + = + + + − ∆ P k xx P xx P r p a k r p k 1 1 ( 1) ( ) ( 1 )， b P P E E = ，組合（3.19）式及（3.20）式；亦即（3.19）式＋α （3.20）式可得                      ∆ +                       ∆ =                                             +                       ₊ + + b P P P P b P b P P P P E E a p a p a a R . . . . 0 . . . . 0 1 ) 1 ( . . . ) ( 0 0 ) ( . . . ) 1 ( 1 * 1 1 1 α α _（3.21）將（3.21）式與下式（3.22）式作比較：

(13)

                      =                       + + + + + + 0 . . . . 0 ) 1 ( ) ( . . . ) 1 ( 1 ₁ 1 1 1 1 P P P P P E p a p a a R （3.22）並由（3.21）式及（3.22）式可以得到下列參數間的關係：從第一行: α =a_P₊₁(p+1) （3.23） P P b P P E E α α + ∆ = + ∆ = + + 1 1 0 （3.24）從最後一行: * 1 1 + + = P + ∆P P E E α （3.25）從其它行: a_P₊₁(k)=a_P(k)+αab_p(p+1−k) k =1_Lp （3.26） α稱為反射係數（reflection coefficients），由（3.23）式至（3.26）式可以分別推得： P P P P E p a E ) 1 ( 1 1 + − = − = ∆ + + α （3.27） P P P E p a 1 1( 1) + + ∆ − = + = α （3.28）

(

P

)

P P P P P E p a E E E E 2 1 2 1 ) 1 ( 1 ) 1 ( * ) ( + − = − = − + = + + α α α （3.29）因此便可求得解自迴歸模型係數

a

P

(k

)

之 Levinson-Durbin 演算法：初始條件： 1 = k ) 0 ( ) 1 ( ) 1 ( 1 xx xx r r a =− （3.30） ) 1 ( 1 =rxx ∆ （3.31）

( ) ( )

[

( ) ( )

]

( )

(

1 1

)

( )

0 0 1 1 0 ) 1 ( ) 0 ( 2 1 * 1 1 * 1 1 1 xx xx xx xx r a r a a r a r E − = − + = ∆ + = （3.32） p k =2_L

∑

= + = + + + − ∆ P k xx P xx P r p a k r p k 1 1 ( 1) ( ) ( 1 ) （3.33）

(14)

(

)

( ) (

)

p p k xx p xx P P P E k p r k a p r E p a

∑

= + + − + + + − = ∆ − = = + 1 1 1 1 1 ) 1 ( α （3.34） p k k p a p a k a k a_P₊₁( )= _P( )+ _P₊₁( +1) *_P( +1− ) =1_L （3.35）

(

)

(

P

)

P P P P P P P

E

p

a

E

2 1 2 * * 1 1

)

1 (

1

1 )

(

+

−

=

−

=

−

+

=

∆

+

=

+ + +

α

（3.36）而 Levinson-Durbin 演算法之過程摘要，可參照圖 3.6。使用自迴歸模型的技術在每一個區域穩態（locally stationary）的區段（segment），並求出其自迴歸模型的係數。得到自迴歸模型的係數後，則以自迴

歸模型近似信號；即x∧(n)的功率頻譜密度(power spectrum density)，P∧

( )

f 便可求

得[37] [39-40] [51]。亦可由（3.2）式得出其表示式為： 2 1 2 ) ( 1 1 ) ( ˆ

∑

= − + = P k k f j P k e a f P π （3.37）由估測信號x∧(n)的功率頻譜密度P∧

( )

f 及信號x

( )

n 的功率頻譜密度P

( )

f ，便

可求得預測信號的全部誤差（total error of prediction）；E，如下表示[40]:

∑

− = ∧ + = m m i i i f P f P m E ) ( ) ( 1 2 1 （3.38）其中 i=-m,...,-1,0,1,...,m,係均勻分佈在[-π,π]之間，2m+1則為計算功率頻譜密度之所有點數。（3.38）式在討論可適性區段方法的章節中是一項重要的參數。 2.倒頻譜模型（Cepstral Model）一個信號 x(n)的倒頻譜原來是要將該信號的頻譜取對數後，再經反傅氏轉換（inverse Fourier transform）才可求得(圖 3.7)。如下式表示：

( )

ω π ω π π ω d e e X c_n

∫

j j n − = log 2 1 （3.39）但當其轉移函數的極點都在單位圓內時，倒頻譜模型的係數可由自迴歸模型係數遞迴演算而得，避開兩次傅氏轉換之繁瑣運算，以加快信號處理的速度而節省時

(15)

間[27][41][46-48][50]。故由 n n nz c z C z H − ∞ =

∑

= ≡ 1 ) ( ] ) ( [ ln （3.40）其中 H(z)為 x(n)的自迴歸模型的系統轉移函數，根據之前自迴歸模型（3.2）式的結果可得到

∑

∞ = − = − =             + 1 1 1 1 ln n n n p k k k z c z a （3.41）將上式等號兩邊對 z -1微分 ) ( 1 1 ln 1 1 1 1

∑

∞ = − − = − − =             + n n n p k k k z c dz d z a dz d （3.42）則 1 1 1 1 1 ) 1 /( ) ( − + ∞ = = − = + −

∑

+ = − n n n p k k k p k k kz a z nc z ka （3.43）或 1 1 1 1 1 ) 1 ( − + ∞ = = − = + −

∑

= + − n n n p k k k p k k kz a z nc z ka （3.44）比較係數

( )

n p n k c a c p n n k c a a c a c c k n n k k n k n n k k n n > − − = ≤ < − − − = − = = − = − − =

∑

1 1 1 0 1 1 1 1 1 0 （3.45） n=1~p即為倒頻譜模型的係數。得到這些倒頻譜模型的係數後，可將其視為原始信號經由自迴歸模型後的系統轉移函數，在倒頻譜領域（Qfrecy domain）的表示式；由於已經模型化，不再是隨機信號，直接執行傅氏轉換，即可求得信號在倒頻譜的功率頻譜密度，針對每一個區段的信號皆可求得其倒頻譜係數與功率頻譜密度，並計算所有區段的最大峰值功率密度比；可以得出其功率頻譜密度(power spectrum density)的表示式如下： 2 2 1 ˆ ( ) p j f n n n P f c e− π = =

∑

（3.46）

(16)

3.可適性區段（adaptive segmentation）此種方法已經被成功的應用在許多非穩態之生物醫學信號的分析處理上 [3][22][27]，而本研究中所量測到的膝關節振動信號其本質為非穩態的隨機程序，且在本研究中所使用之頻譜估測技術皆是針對穩態的信號，所以在應用這些方法之前，必須先將原始的信號以可適性區段的方法 segment 成許多區段，其中每一區段皆可視為區域穩態（locally stationary）。本研究中的可適性區段之方法乃是參考[3]的作法。首先選取一個參考窗戶（reference window）長度為 N 點，並計算此 N 點信號的功率頻譜密度，則自迴歸模型和倒頻譜模型的功率頻譜密度可以由（3.37）式、（3.46）式求得。因此，預測誤差（prediction error）便可由（3.38）式求得，從參考窗戶獲得之自迴歸模型係數被用來估側下一個 N 點信號。如果下一個 N 點信號的預測誤差小於一

個閥值（threshold），則此 N 點被加入目前的參考區段（reference segment）中，

此過程將持續到預測誤差超過閥值為止，以下為可適性區段的步驟[3]： Step1:選取首 N 點信號，並求解此 N 點信號之自迴歸模型模型係數。 Step2:選取下一個 N 點信號，並求出其自迴歸模型係數，使用此自迴歸模型係數我們可以求得預測信號誤差E ，如（3.38）所示： o

( )

∑

− = ∧ + = m m i i o i f P f P m E 1 2 1 0 （3.47）

( )

f P and Po

( )

f ∧ 分別表示下一個 N 點信號之功率頻譜密度及自迴歸模型信號之功率頻譜密度。

Step3:使用 Step 1 中所獲得之自迴歸模型係數去估測 Step 2 中之 N 點信號，並

求出預測信號的全部誤差

E

₁ :

( )

1 1 1 2 1 m i i m i P f E m ₌₋ _P∧ _f = +

∑

（3.48）

( )

f P1 ∧ 為 Step 1 中以自迴歸模型係數求得之功率頻譜密度。 Step4:定義 1 0 E E = β 。如果

β

大於某一個選定的閥值，則視作 Step1 與

Step2 中的信號具有相似的穩態（similar stationary）特性；亦即此二個連續 N 點的信號具有相似穩態特性，因此可將信號的區段點數繼續延伸增

加 N 點，並將運算回到 Step2，若否，即

β

小於這個選定的閥值，則此區

段之邊界（boundary）找到，並且運算回至 Step 1 重新開始下一個信號的邊界尋找。延續之前的研究，所選定之閥值為 0.85[32-33]。

(17)

上述的

β

值可以由如下表示[38]：

∑

− = − = = = m m i i i m m i i i f P f P f P f P E E ) ( ˆ ) ( ) ( ˆ ) ( 1 0 1 0 β （3.49）即

β

可由P( f)、Pˆ₀(f)、Pˆ₁(f)求得。

E 特徵參數

1.根均方值(Root Mean Square；RMS)

首先可直接計算與能量有關的參數：VAM 信號的根均方值，其為將週期時間內之每點信號大小的平方和求其總和，再取平均值後再開根號。振動信號 x(n) 之 RMS，其數學表示式如下[42]:

( )

2 1 1 2 1 RMS     =

∑

= N n n x N （3.50）其中 x(n)為振動信號每個取樣點的值，N 為信號的總點數，且 RMS 之因次為(公尺/秒2 )。 2.內部群組距離（intraclass distance）內部群組距離係表示一個非穩態隨機信號經由可適性區段後，所得到之全部區域性穩態區段，彼此間自迴歸模型係數之緊密程度。令 {a_j

( )

1 ,a_j

( )

2 ,....,a_j

( )

p }為第 j 個區段的自迴歸模型係數，則其所有區段的自迴歸模型係數之平均值可為：

( )

∑

( )

= = nT j j T k a n a 1 1 l （3.51）其中k=1,2,...., p ,而n 為所有區段的數目。內部群組距離即表示每一個區段的自_T 迴歸模型係數與平均之自迴歸模型係數值間的緊密程度，亦即信號之所有區段間的差異性越大，內部群組距離也越大，定義如下[3]：

∑ ∑

(

( ) ( )

)

= =      ₋ = nT j p k j T k a k a n D 1 2 1 1 2 1 （3.52）此值即可表示每個信號中之所有區段間的自迴歸模型係數的離散程度，且其為一無因次的量。

(18)

主極點之意義為使功率頻譜密度最大之相對應頻率，且此頻率亦必須為系統轉移函數的極點之一。在自迴歸模型中，信號經可適性區段之方法後，一個區段中之主極點的功率為 Ω_k =P_k

( )

f_d_p =Max{P_k

( )

f_{( )}_j },k =1,...,n_T;j=1,....,25 （3.53）其中 nT表示一個信號中之區段總數，Ω 為第 k 個區段主極點所對應之功率，k

( )

f P_k 為第 k 個區段的功率頻譜密度， f_{( )}_j 為第 j 個極點對應之頻率， f_d_p為主極點相對應之頻率。則定義主極點功率比為：

∑

= ∈ < < dp ≡ nT k k G k k f f f R 1 2 1 Ω Ω （3.54）其中 G 表示在 nT個區段中，主極點頻率 fdp介於 f1與 f2之間的區段的集合。由於不同的信號其能量在頻譜之分布情況亦不相同，而主極點功率比可反映出分布在某一個頻帶內之能量所佔的比例，因此具有分析的價值。

4.最大峰值功率比（spectral power ratio of maximum peak）

最大峰值之意義為使功率頻譜密度最大之相對應頻率。在倒頻譜模型中，信號經可適性區段之方法後，一個區段中之最大峰值的功率為 Ω_k =P_k(f_m_p)=Max{P_k

( )

f },k =1,...,n_T （3.55）其中 nT表示一個信號中之區段總數，Ω 為第 k 個區段最大峰值之功率，k Pk

( )

f 為第 k 個區段的功率頻譜密度， f_mp為最大峰值相對應之頻率。則定義最大峰值功率比為：

∑

= ∈ < < mp ≡ nT k k G k k f f f R 1 2 1 Ω Ω （3.56）其中 G 表示在 nT個區段中，最大峰值頻率 fmp介於 f1與 f2之間的區段的集合。由於不同的信號其能量在頻譜之分布情況亦不相同，而最大峰值功率比可反映出分布在某一個頻帶內之能量所佔的比例，因此具有分析的價值。至於 PPC 信號之特徵參數，除了主極點功率比及最大峰值功率比方面，需再重新尋找適合之頻帶外，其餘部分與 VAM 完全相同。

F 單方變異數分析（one-way ANOVA）

單方變異數分析是一項統計學上常用的理論，它可以用來檢定多組樣本所

(19)

分別代表之母群體的平均值彼此之間是否有差異性存在。假設現有 k 組樣本

{

x1(1),x1(2),...,x1(n1)

}

，

{

x2(1),x2(2),...,x2(n2)

}

，…，

{

xk(1),xk(2),...,xk(nk)

}

，則第 i 組樣本之平均值為：

∑

= −

=

ni j i i i

x

j

n

x

1

)

(

1

（3.57）而全部樣本之平均值可表示為：

N

x

n

x

k i i i

∑

= − −

=

1 （3.58）其中

∑

= = k i i n N 1

為樣本總數。組內均方（within groups mean squre）的定義為：

( )

2 1 1

k

N

x

j

x

MSW

k i n j i i i

−









₋

≡

∑∑

= = − （3.59）其自由度（degree of freedom）為 N-k。組間均方（between groups mean squre）的定義為： 2 1

1 −











 −

≡

∑

= − −

k

x

n

MSB

k i i i （3.60）而其自由度為 k-1。計算

MSW

MSB

F

=

_（3.61）並且查閱 F 分配（F-distribution）之機率表，若 F 所對應之機率值；p，小於某一個由使用者所自行選定之 p 值，則表示各群組彼此之間有差異。在一般的情況下，只要 F 所對應之機率值小於 0.05，即可宣告各群組之間有差異性。使用變異數分析的結論，只能知道各群組間之平均值並非完全相同，然而這樣的結論，並無法確定是否其中特定的二群組有差異性。倘若群組的數目大於 2，通常都是先對所有的群組進行單方變異數分析之檢定。若檢定的結果顯示各群組之間確實有差異的存在，則再進行兩兩檢定，以確定究竟是哪些群組之間有差異[43-45]。

Ⅳ VAM 信號之分析結果

關節振動量測術(Vibration arthrometry；VAM)是藉由分析膝關節的振動信號

(20)

來診斷膝關節之病變，經由已發表的文獻報告可以發現 VAM 是一種極具發展潛力的非侵襲性診斷工具。且從過去的研究結果得知，VAM 能有效地診斷出膝關節的髕骨與內、外側半月板之一些病變情況[28-33][49][55]。本章將利用自迴歸模型與倒頻譜模型兩種信號處理的數學模式來實驗分析膝關節振動信號。本章研究的主要對象即為前十字韌帶重建術後病變患者；參考圖 4.2。而圖 4.3 為關節鏡直視下之退化性膝關節炎病人髕股關節，圖 4.1 為關節鏡直視下之正常髕股關節。經由對髕骨 VAM 振動信號的數學模式與頻譜分析，以期能找出區分正常

者（Normal）、前十字韌帶重建術後病變患者（ACL）、退化性膝關節炎患者（OA）

之特徵參數。圖 4.4、圖 4.5、圖 4.6 依序分別為正常者、前十字韌帶重建術後病變患者、退化性膝關節炎患者之 VAM 信號圖，由這三圖之時域波形，即可明顯看出其中的差異。即正常者的振幅較小，且上下擺動較為均勻，而前十字韌帶重建術後病變患者與退化性膝關節炎患者則有週期性振幅較大之突波出現，且此二圖之突波週期並不相同。因此由圖 4.4、圖 4.5、圖 4.6 可作為膝關節病變患者 VAM 信號圖之初步比對診斷。以下將針對所量測到的 VAM 信號進行自迴歸模型與倒頻譜模型之分析研究。

A 自迴歸模型數學演算法之結果

1 根均方值表 4.1 列出了 40 位正常者其 VAM 信號的根均方值以及其平均值與標準差，表 4.2 列出了 40 位前十字韌帶重建術後病變患者其 VAM 信號的根均方值以及其平均值與標準差，表 4.3 列出了 52 位退化性膝關節炎患者其 VAM 信號的根均方值以及其平均值與標準差。正常者 40 例 VAM 信號的根均方值其平均值為 0.22574，標準差為 0.09888，前十字韌帶重建術後病變患者 40 例 VAM 信號的根均方值其平均值為 0.36579，標準差為 0.27075，退化性膝關節炎患者 52 例 VAM 信號的根均方值其平均值為 0.53870，標準差為 0.54011。因此可得知退化性膝關節炎患者的根均方值之平均值最大，前十字韌帶重建術後病變患者次之，正常者最小。將此三個群組之“根均方值”進行 one-way ANOVA 的分析，得其 F 值為 8.0716（表 4.13），經查表後發現其所 F 值對應之機率值 p<0.05（表 4.13），因此認為三個群組之間具有差異性。接下來再進行兩兩群組之間的檢定，所得到的 F 值和 P 值，如表 4.13。由表 4.13 可得知，除了（前十字韌帶重建術後病變患者與退化性膝關節炎患者），其 F 值所對應之機率值 p>0.05 不具差異性外；其餘兩兩群組之間，其 F 值所對應之機率值皆小於 0.05，具有差異性。因此由“根均方值”這個特徵參數進行 one-way ANOVA 的分析時，只有（ACL、 OA）二群不可區分，其餘不論三群或二群皆可區分。

2 內部群組距離

(21)

準差，表 4.5 列出了 40 位前十字韌帶重建術後病變患者其 VAM 信號的內部群組距離以及其平均值與標準差，表 4.6 列出了 52 位退化性膝關節炎患者其 VAM 信號的內部群組距離以及其平均值與標準差。正常者 40 例 VAM 信號的內部群組距離其平均值為 0.80283，標準差為 0.04124，前十字韌帶重建術後病變患者 40 例 VAM 信號的內部群組距離其平均值為 0.78373，標準差為 0.06347，退化性膝關節炎患者 52 例 VAM 信號的內部群組距離其平均值為 0.76949，標準差為 0.04233。因此可得知正常者的內部群組距離之平均值最大，前十字韌帶重建術後病變患者次之，退化性膝關節炎患者最小。將此三個群組之“內部群組距離” 進行 one-way ANOVA 的分析，得其 F 值為 5.1582（表 4.14），經查表後發現其所 F 值對應之機率值 p<0.05（表 4.14），因此認為三個群組之間具有差異性。接下來再進行兩兩群組之間的檢定，所得到的 F 值和 P 值，如表 4.14，由表 4.14 可得知，除了（正常者與前十字韌帶重建術後病變患者）及（前十字韌帶重建術後病變患者與退化性膝關節炎患者），其 F 值所對應之機率值 p>0.05 不具差異性外；其餘兩兩群組之間，其 F 值所對應之機率值皆小於 0.05，具有差異性。因此由“內部群組距離”這個特徵參數進行 one-way ANOVA 的分析時，只有（Normal、ACL）與（ACL、OA）二群不可區分，其餘不論三群或二群皆可區分。 3 主極點功率比圖 4.7、圖 4.8、圖 4.9 依序分別為 40 例正常者、40 例前十字韌帶重建術後病變患者、52 例退化性膝關節炎患者其 VAM 信號主極點功率比之平均值分佈圖。由圖中可發現三者的主極點幾乎都集中在 10Hz 以下，但正常者最多，前十字韌帶重建術後病變患者次之、退化性膝關節炎患者最少。且退化性膝關節炎患者在頻帶[550Hz~1000Hz] 及 1000Hz 以上之主極點分佈較前十字韌帶重建術後病變患者來的多。在量取膝關節振動信號的實驗中，所使用之取樣率為 50kHz，為了避免干擾效應（aliasing effect），所以在 A/D convertor 前加了一個 10kHz 的前置濾波器。因此對每個區域穩態區段的頻譜估測之有效頻寛為 0~10kHz。故在 0~10kHz 中的任一組頻率範圍[f1,f2]間，可求得正常者、前十字韌帶重建術後病變患者和退化性膝關節炎患者的主極點功率比，然後將此三群的主極點功率比進行 one-way ANOVA 分析，如果得到之 F 值很高且相對應的機率值 p<0.05，表示在[f1,f2]頻帶間的主極點功率比，可當作有效區分的特徵參數，即表示此群組間具有差異性。在本論文研究中，使用 10Hz 當作最小的區別頻段範圍。所以在 0Hz~10kHz，計有[0Hz,10Hz]，[0Hz,20Hz] ，[0Hz,30Hz] ，……，[0Hz,9990Hz] ，[0Hz,10kHz] 共 1000 種組合；同理，在 10Hz~10kHz，計有[10Hz,20Hz]，[10Hz,30Hz]， [10Hz,40Hz]，……，[10Hz,9990Hz]，[10Hz,10kHz]共 999 種組合，依此類推；因此在 0Hz 至 10kHz 之區段中共有 500,500 種組合，亦即有 500,500 個主極點功率比之特徵參數可以使用。如同第三章所敘述：不同的特徵參數之特性與如何產生膝關節振動信號的生理組織有密切之關聯；即不同的膝關節病變會對應到不

(22)

同之特徵參數。故必須利用 one-way Anova 的分析，來尋求最適合之特徵參數。在 500,500 個組合中，發現[0Hz ~500Hz]這個區段最有意義，因為（正常者, 前十字韌帶重建術後病變患者,退化性膝關節炎患者）三群在[0Hz ~500Hz]這個區段，進行 one-way Anova 分析後，得到 F 值為 19.3665 且其對應之機率值 p<0.05 （表 4.17）。即表示在[0Hz ~500Hz]這個區段間之主極點功率比，可以用來當作可區分之特徵參數。在這個區段內，正常者之平均功率為 99.5559，標準差為 0.74330（表 4.19），前十字韌帶重建術後病變患者之平均功率為 85.5993，標準差為 24.1810（表 4.20），退化性膝關節炎患者之平均功率為 71.6293，標準差為 26.6455（表 4.21），皆含有相當高之主極點功率比，表示此三群不同之膝關節振動信號，主要的資訊（information）皆集中在中頻以下。三群在中頻區間可區分，並不代表任意二群在中頻區間可區分。因此，仍必需進行兩兩群組之間的檢定。經 one-way Anova 自動搜尋最佳頻帶後發現，任意二群在[0Hz ~500Hz]間之主極點功率比皆具有差異性；亦即其 F 值所對應之機率值皆小於 0.05（表 4.17）。因此，不論三群或二群，皆可由[0Hz ~500Hz]間之主極點功率比來當作區別群組的特徵參數；此最佳頻帶為中頻帶區間。

B 倒頻譜模型數學演算法之結果

1 根均方值因根均方值為時域中之特徵參數，不會因為數學演算法之改變而有所不同，故其分析之結果與 4.2.1 完全相同。在此將其分析之結果整理在表 4.7、表 4.8、表 4.9、表 4.15 中。 2 內部群組距離表 4.10 列出了 40 位正常者其 VAM 信號的內部群組距離以及其平均值與標準差，表 4.11 列出了 40 位前十字韌帶重建術後病變患者其 VAM 信號的內部群組距離以及其平均值與標準差，表 4.12 列出了 52 位退化性膝關節炎患者其 VAM 信號的內部群組距離以及其平均值與標準差。正常者 40 例 VAM 信號的內部群組距離其平均值為 0.47925，標準差為 0.02367，前十字韌帶重建術後病變患者 40 例 VAM 信號的內部群組距離其平均值為 0.50009，標準差為 0.02756，退化性膝關節炎患者 52 例 VAM 信號的內部群組距離其平均值為 0.51653，標準差為 0.02612。因此可得知退化性膝關節炎患者的內部群組距離之平均值最大，前十字韌帶重建術後病變患次之，正常者最小。將此三個群組之“內部群組距離”進行 one-way ANOVA 的分析，得其 F 值為 23.4863（表 4.16），經查表後發現其所 F 值對應之機率值 p<0.05（表 4.16），因此認為三個群組之間具有差異性。接下來再進行兩兩群組之間的檢定，所得到的 F 值和 P 值，如表 4.16。由表 4.16 可得知，任何兩兩群組之間，其 F 值所對應之機率值皆小於 0.05，因此都具有差異性。故由“內部群組距離”這個特徵參數進行 one-way ANOVA 的分析時，不論三群或二群皆可區分。

(23)

3 最大峰值功率比圖 4.10、圖 4.11、圖 4.12 依序分別為 40 例正常者、40 例前十字韌帶重建術後病變患者、52 例退化性膝關節炎患者其 VAM 信號最大峰值功率比之平均值分佈圖。由圖中可發現三者的最大峰值幾乎都集中在 10Hz 以下，但正常者最多，前十字韌帶重建術後病變患者次之、退化性膝關節炎患者最少。且退化性膝關節炎患者在頻帶[550 Hz ~1000 Hz] 之最大峰值分佈較前十字韌帶重建術後病變患者來的多。在 0Hz 至 10k 之區段中共有 500,500 個最大峰值功率比之特徵參數可以使用。故必須利用 one-way Anova 的分析，來尋求最適合之特徵參數。在 500,500 個組合中，發現[0Hz ~150Hz]這個區段最有意義，因為（正常者,前十字韌帶重建術後病變患者,退化性膝關節炎患者）三群在[0Hz ~150Hz]這個區段，進行 one-way Anova 分析後，得到 F 值為 4.5241 且其對應之機率值 p<0.05（表 4.18）。即表示在[0Hz ~150Hz]這個區段間之最大峰值功率比，可以用來當作可區分之特徵參數。在這個區段內，正常者之的平均功率為 83.0990，標準差為 9.1899（表 4.22），前十字韌帶重建術後病變患者的平均功率為 78.8602，標準差為 10.1493（表 4.23），退化性膝關節炎患者的平均功率為 76.2465，標準差為 12.3681（表 4.24），皆含有相當高之最大峰值功率比，表示此三群不同之膝關節振動信號，主要的資訊（information）皆集中在低頻。三群在低頻區間可區分，並不代表任意二群在低頻區間可區分。因此，仍必需進行兩兩群組之間的檢定。經 one-way Anova 自動搜尋最佳頻帶後發現；（正常者,前十字韌帶重建術後病變患者）在[1180Hz ~1490Hz]間之最大峰值功率比具有差異性，（正常者,退化性膝關節炎患者）在[0Hz ~150Hz] 間之最大峰值功率比具有差異性，（前十字韌帶重建術後病變患者,退化性膝關節炎患者）在[960Hz ~1270Hz] 間之最大峰值功率比具有差異性；且其任意二群的 F 值所對應之機率值皆小於 0.05（表 4.18），並將各個二群組可區別頻帶之最大峰值功率比的平均值與標準差分別列於表 4.25、表 4.26、表 4.22、表 4.24、表 4.27、表 4.28。因此，不論三群或二群，皆可找到最佳頻帶之最大峰值功率比當作區別群組的特徵參數。其中（正常者,前十字韌帶重建術後病變患者,退化性膝關節炎患者）三群及（正常者,退化性膝關節炎患者）二群之最佳頻帶為低頻帶區間；而（正常者,前十字韌帶重建術後病變患者）二群與（前十字韌帶重建術後病變患者,退化性膝關節炎患者）二群之最佳頻帶為高頻帶區間。

PPC 信號之分析結果

髕股顫振（physiological patellofemoral crepitus；PPC）是一種由正常膝關節

在慢速（小於 50/秒）主動擺動下所產生在髕骨之振動信號，而此一信號將可反

映出關節軟骨之完好情況。

本研究在第四章中已成功的利用膝關節在快速擺動下所產生的振動信號

(24)

（ACL）與退化性膝關節炎患者（OA）的特徵參數。在本章的研究中，仍將採用和第四章相同的信號處理之方法，利用自迴歸模型和倒頻譜模型的分析方法，找出足以區分出正常者、前十字韌帶重建術後病變患者與退化性膝關節炎患者的特徵參數。圖 5.1、圖 5.2、圖 5.3 依序分別為正常者、前十字韌帶重建術後病變患者、退化性膝關節炎患者之 PPC 信號圖，由這三圖之時域波形，即可明顯看出其中的差異。即正常者的振幅於上下擺動時較為均勻綿密，而前十字韌帶重建術後病變患者與退化性膝關節炎患者則有不規則性之突波出現，且前十字韌帶重建術後病變患者之突波振幅較大。因此由圖 5.1、圖 5.2、圖 5.3 可作為膝關節病變患者 PPC 信號圖之初步比對診斷。以下將針對所量測到的 PPC 信號進行自迴歸模型與倒頻譜模型之分析研究。

A 自迴歸模型數學演算法之結果

1 根均方值表 5.1 列出了 40 位正常者其 PPC 信號的根均方值以及其平均值與標準差，表 5.2 列出了 40 位前十字韌帶重建術後病變患者其 PPC 信號的根均方值以及其平均值與標準差，表 5.3 列出了 52 位退化性膝關節炎患者其 PPC 信號的根均方值以及其平均值與標準差。正常者 40 例 PPC 信號的根均方值其平均值為 0.12689，標準差為 0.08161，前十字韌帶重建術後病變患者 40 例 PPC 信號的根均方值其平均值為 0.09997，標準差為 0.05497，退化性膝關節炎患者 52 例 PPC 信號的根均方值其平均值為 0.07271，標準差為 0.07233。因此可得知正常者的根均方值之平均值最大，前十字韌帶重建術後病變患者次之，正常者退化性膝關節炎患者最小。將此三個群組之“根均方值”進行 one-way ANOVA 的分析，得其 F 值為 6.6828（表 5.13），經查表後發現其所 F 值對應之機率值 p<0.05（表 5.13），因此認為三個群組之間具有差異性。接下來再進行兩兩群組之間的檢定，所得到的 F 值和 P 值，如表 5.13。由表 5.13 可得知，除了（正常者與退化性膝關節炎患者），其 F 值所對應之機率值 p<0.05 具有差異性外；其餘兩兩群組之間，其 F 值所對應之機率值皆大於 0.05，不具有差異性。因此由“根均方值”這個特徵參數進行 one-way ANOVA 的分析時，只有三群及（Normal、OA）二群可區分，其餘二群皆不可區分。 2 內部群組距離表 5.4 列出了 40 位正常者其 PPC 信號的內部群組距離以及其平均值與標準差，表 5.5 列出了 40 位前十字韌帶重建術後病變患者其 PPC 信號的內部群組距離以及其平均值與標準差，表 5.6 列出了 52 位退化性膝關節炎患者其 PPC 信號的內部群組距離以及其平均值與標準差。正常者 40 例 PPC 信號的內部群組距離其平均值為 1.77250，標準差為 0.56168，前十字韌帶重建術後病變患者 40 例 PPC 信號的內部群組距離其平均值為 1.40345，標準差為 0.30302，退化性膝關節炎患者 52 例 PPC 信號的內部群組距離其平均值為 1.32100，標準差為 0.23504。因此

(25)

可得知正常者的內部群組距離之平均值最大，前十字韌帶重建術後病變患者次之，退化性膝關節炎患者最小。將此三個群組之“內部群組距離”進行 one-way ANOVA 的分析，得其 F 值為 17.1408（表 5.14），經查表後發現其所 F 值對應之機率值 p<0.05（表 5.14），因此認為三個群組之間具有差異性。接下來再進行兩兩群組之間的檢定，所得到的 F 值和 P 值，如表 5.14，由表 5.14 可得知，除了（前十字韌帶重建術後病變患者與退化性膝關節炎患者），其 F 值所對應之機率值 p>0.05 不具差異性外；其餘兩兩群組之間，其 F 值所對應之機率值皆小於 0.05，具有差異性。因此由“內部群組距離”這個特徵參數進行 one-way ANOVA 的分析時，只有（ACL、OA）二群不可區分，其餘不論三群或二群皆可區分。 3 主極點功率比圖 5.4、圖 5.5、圖 5.6 依序分別為 40 例正常者、40 例前十字韌帶重建術後病變患者、52 例退化性膝關節炎患者其 PPC 信號主極點功率比之平均值分佈圖。由圖中可發現正常者與前十字韌帶重建術後病變患者的主極點幾乎都集中在 10Hz 以下及[100Hz~450Hz]間，但正常者在 10Hz 以下明顯分佈較多；而退化性膝關節炎患者的主極點分佈則是由低頻向高頻逐漸遞減。且退化性膝關節炎患者在頻帶[10Hz~100Hz] 及[450Hz~550Hz] 間之主極點分佈較前十字韌帶重建術後病變患者來的多。在量取膝關節振動信號的實驗中，所使用之取樣率為 5kHz，為了避免干擾效應（aliasing effect），所以在 A/D convertor 前加了一個 2kHz 的前置濾波器。因此對每個區域穩態區段的頻譜估測之有效頻寛為 0~2kHz。故在 0~2kHz 中的任一組頻率範圍[f1,f2]間，可求得正常者、前十字韌帶重建術後病變患者和退化性膝關節炎患者的主極點功率比，然後將此三群的主極點功率比進行 one-way ANOVA 分析，如果得到之 F 值很高且相對應的機率值 p<0.05，表示在[f1,f2]頻帶間的主極點功率比，可當作有效區分的特徵參數，即表示此群組間具有差異性。在本論文研究中，使用 10Hz 當作最小的區別頻段範圍。所以在 0Hz~2kHz，計有[0Hz,10Hz]，[0Hz,20Hz] ，[0Hz,30Hz] ，……，[0Hz,1990Hz] ，[0Hz,2kHz] 共 200 種組合；同理，在 10Hz~2kHz ，計有 [10Hz,20Hz] ， [10Hz,30Hz] ， [10Hz,40Hz]，……，[10Hz,1990Hz]，[10Hz,2kHz]共 199 種組合，依此類推；因此在 0Hz 至 2kHz 之區段中共有 20,100 種組合，亦即有 20,100 個主極點功率比之特徵參數可以使用。在 20,100 個組合中，發現[0Hz ~490Hz]這個區段最有意義，因為（正常者, 前十字韌帶重建術後病變患者,退化性膝關節炎患者）三群在[0Hz ~490Hz]這個區段，進行 one-way Anova 分析後，得到 F 值為 20.8873 且其對應之機率值 p<0.05 （表 5.17）。即表示在[0Hz ~490Hz]這個區段間之主極點功率比，可以用來當作可區分之特徵參數。在這個區段內，正常者之平均功率為 98.9768，標準差為 0.84620（表 5.19），前十字韌帶重建術後病變患者之平均功率為 97.9025，標準差為 3.60680（表 5.20），退化性膝關節炎患者之平均功率為 86.9145，標準差為 14.8193（表 5.21），皆含有相當高之主極點功率比，表示此三群不同之膝關節振

十字韌帶重建術後膝部異常聲音之訊號分析

行政院國家科學委員會專題研究計畫 成果報告