基于MATLAB GUI編程的血流/溫度信號時頻分析軟件開發

唐元梁，賀纓

大連理工大學 能源與動力學院，遼寧 大連 116024

引言

微循環是人體血液循環系統的重要組成部分，其基本功能是進行血液和組織液的物質交換，因此健全的微循環功能是人體重要臟器功能正常運行的重要保障。許多疾病的發生，包括高血壓、糖尿病及心腦血管疾病等，均與微循環障礙有關。若能及時檢測出微循環改變的危險信號，將對微循環疾病的早期診斷和治療有極大幫助。

通過對人體皮膚血流信號進行頻譜分析，研究者發現其包含多個波動特征頻率[1-4]（圖1）。雖然這些波動并非是嚴格周期性的，但很顯然也不是隨機出現，表現出類周期性的特點。通過干預交感神經活動、內皮活動等，研究者發現這些波動特征頻率與微循環的調節活動有關，這些調節活動及其活動周期范圍如下：內皮調節（50～105 s）；神經調節（20～50 s）；肌源性調節（7～20 s）；呼吸（2～7 s），心跳（0.5～2 s）。因此，若能通過濾波提取出與各調節活動相對應的波動分量并進行數據分析，可以更加準確地判斷微循環調節功能病變的發生原因。

圖1 人體皮膚血流信號的波動特征

由于皮膚溫度受皮下微循環血液灌注的影響，因此能夠間接地反映皮下微循環的狀況，且相較于血流信號的采集（一般采用激光多普勒血流儀等），體表溫度采集更容易、測量穩定性更好。而且實驗分析顯示[5-7]，人手指尖等部位的皮膚溫度也表現出與血流信號類似的波動特征，且兩者具有一定的波動相關性。因此，綜合考慮微循環調節活動的類周期性特征、皮膚溫度與血流的波動相關性以及熱信號采集方面的優勢，Frick等[8]指出，皮膚溫度的波動特征可以用來標記微循環的調節功能。近年來微循環功能的熱波動檢測逐漸引起的研究者的關注。Smirnova等[9]發現在單側手接受冷水刺激時，健康人群的對側手皮膚溫度在內皮、神經、肌源性調節頻段的波動幅值均減小，冷刺激結束后波動幅值迅速回升至正常；而二型糖尿病、糖耐量受損人群在結束冷刺激后，3個頻段內的溫度波動幅值均沒有恢復正常，表明二型糖尿病組、糖耐量受損組的皮膚微循環調節功能均出現了障礙。研究者發現指尖局部熱刺激時健康人群的皮膚溫度在肌源、神經以及內皮調節頻段的波動幅值都有明顯的增加[10-12]，這反映了熱刺激誘導的血管舒張調節變化；而二型糖尿病患者的指尖在熱刺激前后其內皮調節頻段的溫度波動幅值并沒有顯著變化，表明后者可能存在內皮調節功能障礙。

微循環功能熱波動分析方法可以應用于相關軟、硬件產品的開發。FM Diagnostics公司與俄羅斯科學院連續介質研究所共同研發了Microtest儀器[12]，見圖2，該儀器主要包括了機身、探頭兩部分。探頭部分尺寸約20 mm×30 mm×10 mm，溫度傳感器（HRTS-5760，Honeywell International，Inc.，Morristown，NJ，USA） 安裝在探頭圓孔（直徑約8 mm）的中央，傳感器外周有4個小型電阻加熱器，探頭內填充低熱導率材料來減少環境變換的影響。經信號放大和A/D轉換后溫度測量精度可達0.001℃。將人體待檢測部位皮膚緊貼Microtest溫度傳感器，啟動儀器后可以對皮膚進行變功率加熱和皮膚溫度采集，然后通過配套的軟件自動對獲取的皮膚溫度信號進行時頻分析，并根據不同功率熱刺激下皮膚溫度波動幅值的變化情況對皮下微循環功能進行評價。Sagaidachnyi等[7,13-15]開發了TBF Converter軟件，該軟件提供了信號的Fourier分析和小波分析、PPG信號包跡線提取功能，并能根據不同模型將血流和溫度信號進行相互轉化（圖3）。

圖2 Microtest儀器

利用血流或溫度信號的波動特征對微循環調節功能進行評價時，均首先需要對信號進行頻譜分析。信號分析中，傅里葉變換是常用的方法，但經過傅里葉變換后時域信息丟失，因此傅里葉變換只在頻域里有局部分析能力，采用標準傅里葉變換不能得到信號波動的時變特征。短時傅里葉變化和小波變換均能保留時域信息，但相較于窗口寬度固定、只有單一的分辨率的短時傅里葉變化，小波變換采用可變的時間窗，即低頻信號可采用較寬的時間窗，而高頻信號采用較窄的時間窗，因此可以在時、頻兩域都具有表征信號局部特征的能力。小波變換非常適合于探測正常信號中夾帶的瞬態反常現象并展示其成分，譽為信號分析的“顯微鏡”。作為一種具有多分辨率特點的時間—尺度分析方法，小波變換是分析微循環信號（或皮膚溫度信號）不同頻率波動成分及其變化特征的非常有用的數學工具。

從目前情況來看，盡管基于熱波動分析方法對皮下微循環調節功能進行檢測的前景非常可觀，但相關的軟、硬件的發展還處于比較早期的階段，距臨床應用還有較大距離。目前已有的相關產品數量較少，功能較為簡單。如Microtest儀器對于微循環檢測的模式只有加熱刺激這一種，而且其配套軟件只能處理該儀器獲取的數據，不能夠從外部導入；TBF Converter軟件側重于信號分析和轉換，沒有信息挖掘功能，不能直接用于微循環功能的評價。盡管小波分析方法非常適合于微循環信號分析，但其實現過程比傳統的傅立葉分析更復雜一些，如小波函數的選取、尺度與頻率的轉換以及邊緣效應等問題均可能會對研究者帶來困擾。雖然MATLAB平臺提供的Wavelet Toolbox工具箱能極大地簡化用戶進行信號小波分析的過程，但該工具箱僅包含基本的小波變換功能，后續的數據處理，如不同頻段上信號波動分量的提取、波動幅值的分段統計，不同數據間的波動相關性分析，以及數據處理結果的圖像/表格顯示等仍需用戶自己編程實現。

圖3 TBF Converter軟件界面

因此，為便于使用小波分析方法對微循環血流/溫度信號進行波動分析，本文基于MATLAB平臺的guide工具箱開發了一個有圖像用戶界面的應用軟件，期望對血流/溫度信號的小波分析、微循環功能的檢測提供幫助。

1 生物信號的小波分析方法

1.1 一維連續小波變換

小波分析將信號分解成一系列小波函數的疊加，而這些小波函數都是由一個母小波函數經過平移與伸縮得來的（圖4）。母小波函數φ(t)經伸縮、平移后可得小波序列φa,b(t)：

其中a為尺度因子，b為平移因子，兩者均為實數。對信號x(t)進行連續小波變換后得到小波系數Wx(a,b)：

式中T1、T2分別為時域信號x(t)的起、止時間，Wx(a,b)為信號x(t)的小波系數，φ?為φ的共軛函數。尺度因子a與頻率f呈倒數關系：

其中fcen為母小波函數φ(t)的中心頻率，Tsam為信號x(t)的采樣周期。在信號分析領域，Morlet小波是最常用的小波函數，其具有良好的時-頻定位性質，時域下Morlet復小波函數的簡化形式為：

圖 4 一維連續小波變換圖示

1.2 信號頻譜特征分析

類似于傅立葉分析，小波分析也可以提取信號的頻譜特征。小波系數的模方表征了信號的小波功率譜，將小波功率譜進行時間平均可得到某時間段上信號的小波全譜M(a)：

小波全譜M(a)反映了在尺度a上信號的波動強度，根據頻率與尺度a的關系[式(3)]，小波全譜圖的局部峰值代表了信號的波動特征頻率（或周期）。

1.3 波動相關性分析

信號x(t)、y(t)的整體互相關性可以利用如下的互相關函數計算：

1.4 波動分量提取

利用信號的小波系數Wx(a,b)進行小波反變換，可以提取出原始信號在某一頻段上的波動分量，如對頻段[f1, f2]上信號x(t)的波動分量x1(t)進行提取：

其中[a2, a1]為頻段[f1, f2]對應的尺度因子范圍。提取出波動分量x1(t)后，可以利用其標準差來表征該頻段上信號波動幅值的大小。

1.5 信號時頻分析軟件的開發

為了便于對血流/溫度信號進行基于小波變換的時頻分析，本文通過MATLAB GUI編寫了可視化的數學程序，并生成可執行的軟件（exe文件）。本文所設計開發的軟件主要包括了文件模塊、小波分析模塊、檢測結果模塊以及相關性分析模塊4部分，見圖5。

圖5 軟件模塊

文件模塊主要功能是導入血流、溫度等信號并進行顯示；相關性分析模塊主要功能是計算導入的血流、溫度信號的互相關性；小波分析模塊可以對導入的信號進行小波變換，并提取出不同頻段的信號波動分量幅值變化信息；若導入的血流/溫度信號是血管反應性測試實驗中獲取的數據，如反應性充血、熱刺激和冷刺激等，則根據信號波動分量的幅值變化信息，可以在檢測結果模塊中統計出信號在不同實驗階段的波動幅值變化情況。

在MATLAB軟件的命令行中輸入“guide”，便可創建新的GUI。將編寫的所有GUI添加到一個工程目錄下并編譯通過后，便可以打包生成可執行程序，雙擊該程序便可運行。不過需注意的是MATLAB GUI編程所得的exe軟件的正常運行依賴于MATLAB數學函數庫，因此若計算機上未安裝MATLAB軟件，需要安裝MCR（MATLAB Compiler Runtime）編譯器。

2 軟件使用介紹

2.1 信號小波變換和波動幅值提取

軟件的小波分析模塊界面如圖6所示。其中區域①顯示了導入的血流時間序列信號，該信號為利用ADI公司的ML191激光多普勒血流儀采集，采樣周期0.001 s；信號采樣周期默認為1.0 s，小波分析尺度默認為1～100，根據信號的實際采樣周期和小波分析尺度要求，在“采樣周期”和“分析尺度”后方的文本編輯框內可修正相應的數值；點擊“小波時頻分析”按鈕，軟件對導入的溫度信號進行小波分析，并對分析結果進行顯示；區域②可以顯示某段時間內信號的小波系數圖和小波系數模圖，圖像模式可在左上方的彈出菜單選擇，時間參數可在“開始時間”和“分析時長”后方的文本編輯框內設置，圖中顯示的是區域①所示信號在15～21 min內（對應區域①中兩條藍線間的信號部分）的小波系數圖；區域③顯示了信號的小波方差圖；區域④的3條曲線由上到下分別為所輸入信號在內皮、神經和肌源性頻段的波動幅值變化，默認的尺度因子變化范圍與人體微循環3個主動調節活動的頻段相一致，如需調整則可分別在“最大尺度”和“最小尺度”方的文本編輯框內修改。

2.2 檢測結果的分析

在某些機械、溫度刺激條件下，正常皮膚微循環的調節活動強度會發生變化，而有微循環功能障礙的人群其調節活動響應可能會出現異常。這些調節反應的狀況可以通過皮膚血流/溫度信號波動幅值的變化進行評估。

本軟件在“檢測結果”模塊中，在對不同刺激實驗（熱刺激、對側手冷水刺激、反應性充血）獲取的溫度信號進行小波分析后，將提取的波動幅值進行了分段統計；根據不同實驗階段的溫度波動幅值變化并參考已有文獻中的溫度波動變化指標，可以對（對側手冷水/熱刺激下）皮下微循環狀態進行定性評價。軟件的檢測結果模塊的界面，見圖7。圖中顯示的是一人手指尖熱刺激實驗獲取的皮膚溫度信號在不同實驗階段的平均波動幅值變化情況。該結果顯示Stage 1階段各頻段的平均溫度波動幅值均高于Stage 2階段，表明微循環調節活動強度變化對熱刺激敏感，調節功能正常。

除溫度刺激實驗檢測結果分析外，該模塊還提供了反應性充血實驗數據的分析功能。利用大鼠反應性充血實驗，本研究小組發現二型糖尿病組大鼠在血管閉塞階段內皮、神經調節頻段的皮膚溫度波動幅值小于對照組大鼠且差異明顯[16]；刺激結束后，二型糖尿病組大鼠神經、肌源性頻段的皮膚溫度波動幅值大于對照組大鼠且差異明顯。該實驗結果對于評價人體糖尿病微循環障礙有一定價值，但還需大量的人體實驗進行驗證。

圖6 小波分析模塊

圖7 檢測結果模塊

2.3 信號間的相關性分析

圖8 互相關分析模塊

軟件的相關性分析模塊的界面，見圖8。其中血流、溫度圖的橫坐標為時間，單位s；小波互相關譜圖的橫坐標為波動周期，單位s。圖中紅色、綠色曲線分別為反應性充血實驗獲取的相同位置指尖皮膚血流和溫度信號，其中皮膚血流信號采用激光多普勒血流儀采集，溫度信號采用Agilent溫度記錄儀采集，精度為0.001℃。經互相關分析后求得兩信號的整體互相關系數為0.7385；黑色曲線為兩信號在1～100尺度范圍內的小波互相關譜，從整體上看兩信號在較大尺度上波動相關性較高，表明該溫度信號能夠表征血流信號的低頻脈動特征。

3 總結

利用MATLAB GUI編程，本文基于小波分析方法開發了一款用于微循環血流/溫度信號時頻分析的軟件。在小波分析模塊中，軟件提供了小波變換功能，可以對信號的頻譜特征進行分析；除了基本的小波變換功能外，本軟件還能夠提取信號在不同頻段的波動幅值變化信息，并能夠對導入的血流/溫度信號進行波動相關性分析。由于在微循環調節活動影響下人體皮膚血流信號呈現出類周期性特征，因此皮膚血流信號波動幅值變化信息能夠反映皮下微循環調節活動是否出現異常。而當血流和溫度信號有較高波動相關性時，皮膚溫度信號波動幅值的變化也能用于微循環功能的檢測。因此，結合該軟件各模塊的不同功能，可以為微循環的熱波動檢測提供較為直觀、可信的數據信息。

該軟件還需進一步完善。目前我們在檢測結果模塊中只對信號波動幅值變化進行了分段統計，還需更多量化指標以全面評價微循環。我們將進一步通過動物和臨床實驗獲取更多數據信息，并通過數據分析來逐步優化軟件，以適應將來微血管功能評價的臨床應用需求，為糖尿病足等的早期診斷提供幫助。

[參考文獻]

[1] Kvernmo HD,Stefanovska A,Bracic M,et al.Spectral analysis of the laser Doppler perfusion signal in human skin before and after exercise[J].Microvasc Res,1998,56:173-182.

[2] Bracic M,Stefanovska A.Wavelet-based analysis of human blood- flow dynamics[J].B Math Biol,1998,60(5):919-935.

[3] Kvernmo HD,Stefanovska A,Kirkeb?en KA,et al.Oscillations in the human cutaneous blood perfusion signal modified by endothelium-dependent and endothelium-independent vasodilators[J].Microvasc Res,1999,57(3):298-309.

[4] Soderstrom T,Stefanovska A,Veber M,et al.Involvement of sympathetic nerve activity in skin blood flow oscillations in humans[J].Am J Physiol Heart C,2003,284(5):1638-1646.

[5] Frick P.Wavelet-based correlations of skin temperature and blood flow oscillations[J].Cardiovasc Eng,2008,8(3):185-189.

[6] Bandrivskyy A,Bernjak A,Mcclintock P.Wavelet phase coherence analysis: Application to skin temperature and blood flow[J].Cardiovasc Engin,2004,4:89-93.

[7] Sagaidachnyi AA,Skripal AV,Fomin AV,et al.Determination of the amplitude and phase relationships between oscillations in skin temperature and photoplethysmography-measured blood flow in fingertips[J].Physiol Meas,2014,35(2):153-166.

[8] Frick P,Mizeva I,Podtaev S.Skin temperature variations as a tracer of microvessel tone skin temperature variations as a tracer of microvessel tone[J].Biomed Signal Proces,2015,21:1-7.

[9] Smirnova E,Podtaev S,Mizeva I,et al.Assessment of endothelial dysfunction in patients with impaired glucose tolerance during a cold pressor test[J].Diab Vasc Dis Res,2013,10(6):489-497.

[10] Parshakov A,Zubareva N,Podtaev S,et al.Local heating test for detection of microcirculation abnormalities in patients with diabetes-related foot complications[J].Adv Skin Wound Care,2017,30(4):158.

[11] Podtaev S,Stepanov R,Smirnova E,et al.Wavelet-analysis of skin temperature oscillations during local heating for revealing endothelial dysfunction[J].Microvasc Res,2015,97:109-114.

[12] Parshakov A,Zubareva N,Podtaev S,et al.Detection of endothelial dysfunction using skin temperature oscillations analysis during local heating in patients with peripheral arterial disease[J].Microcirculation,2016,23:406-415.

[13] Sagaidachnyi AA,Usanov DA,Skripal AV,et al.Correlation of skin temperature and blood flow oscillations[J].Proceedings of SPIE,2012,8337(4):9.

[14] Sagaidachnyi A,Usanov DA,Skripal AV,et al.Restoration of finger blood flow oscillations by means of thermal imaging[A].11thInternational Conference on Quantitative Infra Red Thermography[C].2012.

[15] Sagaidachnyi AA,Usanov DA,Skripal AV,et al.Thermographybased blood perfusion imaging in hands : spectral amplification and time shift[A].Quantitative Infrared Thermography[C].2014.

[16] 唐元梁,糖尿病微血管病變熱檢測方法的機理研究[D].大連:大連理工大學,2017.