強迫振蕩技術（FOT）是一種主動式肺功能測量技術，通過給呼吸系統外加激勵信號的方式，辨識呼吸系統力學特性。FOT 技術常用的激勵信號包括單頻正弦、偽隨機和周期性脈沖三種。本文針對偽隨機多正弦組合信號存在的時域幅度過沖問題，研究了偽隨機信號的相位優化，嘗試了隨機相位組合以及時-頻域交換算法，以波峰因數評價優化效果。進一步根據振蕩單元的頻率響應曲線對低頻段（4～18 Hz）輸入信號的幅值進行了補償，使最終產生的偽隨機信號的時-頻域特性在 4～40 Hz 范圍內能夠滿足呼吸系統辨識要求。研究結果表明，時-頻域交換算法能夠有效優化偽隨機信號的相位組合，振蕩單元的幅頻特性經過低頻補償后，能夠產生滿足性能要求的激勵信號。
引用本文： 劉曉莉, 張楠, 梁洪, 張政波, 李德玉, 王衛東. 強迫振蕩肺功能的偽隨機激勵信號優化研究. 生物醫學工程學雜志, 2017, 34(5): 660-666. doi: 10.7507/1001-5515.201607064 復制
引言
強迫振蕩技術（forced oscillation technique，FOT）是一種主動測量呼吸系統力學特性的技術，其核心思想是將人體的呼吸系統視為一個“黑箱”，利用系統辨識的原理，通過外界施加不同類型的壓力振蕩，測量呼吸系統的響應，即壓力、流量，進而辨識出呼吸系統阻抗[1-3]。FOT 是一種無創技術，無需受試者主動配合，測量過程中只需要受試者平靜地進行自主呼吸，即可非常靈敏地測試出呼吸系統力學特性參數，特別適合于老年人、兒童使用[4]；它可以分辨吸氣相和呼氣相阻力變化，對于辨別慢性阻塞性肺疾病和哮喘具有重要意義[5-8]；另外，它還有助于對吸煙患者呼吸系統病變的辨識，可以幫助臨床早期診斷由吸煙導致的呼吸功能改變[9]。因而，FOT 在呼吸系統疾病診斷和指導用藥方面具有廣闊的應用前景[3]。
FOT 通過壓力振蕩源施加連續的激勵信號以檢測系統的響應[3]。FOT 常采用的激勵信號有三種，分別為單頻正弦信號、偽隨機信號（多正弦組合）和周期脈沖信號。單頻正弦信號因容易實現、信噪比高以及對硬件部分要求較低，在 FOT 的技術早期得到了廣泛使用。但單頻正弦信號若要實現多個頻率點的呼吸阻抗測量，需要依次掃頻，存在檢測時間偏長和時間分辨率較低的問題。后來，相繼出現了偽隨機多正弦組合和周期性脈沖等激勵信號，它們能夠在較短的時間內產生呼吸系統辨識所需的持續激勵信號，使呼吸系統在多個頻率點被充分激勵，因此能夠有效縮短檢測時間，并具有較高的時間分辨率，但該方法對硬件電路和算法提出了更高要求[3, 10-12]。本文著重介紹偽隨機多正弦組合的相關研究，后續文章會介紹周期性脈沖信號的研究進展。
FOT 的偽隨機多正弦組合激勵信號本質上是多個正弦信號的疊加，存在相位組合問題。多正弦信號組合在頻域特性滿足要求的情況下，不同的相位組合會產生時域特性不同的信號，若相位均相同，將出現最為極端的情況，時域信號周期性疊加，信號幅值太大，極易導致系統飽和失真[3, 12]。另一方面，呼吸系統對于低頻激勵信號響應更為敏感，可以反映更多正常生理過程和病理結構改變的內容，因而臨床更為關注[7, 13-14]。而 FOT 振蕩單元產生的激勵信號的幅頻特性曲線存在低頻響應幅度過低的問題，因此需要對輸入激勵信號的低頻段的幅值進行補償。本文針對強迫振蕩技術肺功能儀研制過程中遇到的偽隨機多正弦組合信號產生、信號低頻段幅度補償以及相位優化問題進行了深入研究。
1 原理
1.1 振蕩肺功能基本原理
FOT 的經典裝置示意圖如圖 1 所示，壓力振蕩單元通常為揚聲器，產生單頻或者復合頻振蕩的壓力信號，加載在呼吸道上，受試者在測試過程中進行自主呼吸，同步采集自主呼吸過程中的壓力信號 Pao 和流量信號 V′[1-2]。運用系統辨識原理，將呼吸系統等效為一個“黑箱”，輸入信號為壓力，輸出信號為流量，去除自主呼吸產生的干擾，壓力和流量信號分別做傅里葉變換后相比得到呼吸系統的總阻抗 Zrs。
圖1
強迫振蕩肺功能系統示意圖
Figure1.
Schematic arrangement of the forced oscillatory pulmonary function system
1.2 振蕩單元的幅頻特性及補償
基于上述經典裝置示意圖，實驗室設計了振蕩肺功能儀，其振蕩單元包含揚聲器和功率放大器。其中揚聲器采用 JBL 的 145 型超低音喇叭，額定功率 200 W，阻抗 8 Ω，靈敏度為 92 dB/W。本研究首先測試了振蕩單元整體的幅頻特性，利用上位機軟件（LabVIEW2013）分別產生幅值為 300 mV、2～50 Hz（2 Hz 等間隔）正弦掃頻的數字信號，經過 NI 公司的 DAQ 設備 USB-6366 的 D/A 轉換成模擬信號，經功率放大器后驅動揚聲器振蕩[15]。輸入幅值為 300 mV，是經過前期振蕩單元性能測試以產生開路測量壓力為 2 cm H2O 的最佳輸入幅值。裝置輸出端空載時測量得到系統響應的壓力信號，如圖 2 所示。從圖中可以明顯得出，由于揚聲器本身的頻響特性以及低頻聲波傳播時在管路開口處能量泄漏要高于高頻聲波[16]，因此在輸入信號幅值相同的情況下，系統輸出壓力信號幅度在 2～18 Hz 頻率范圍內明顯衰減，會對我們后續研究呼吸系統在低頻段的激勵響應帶來很大影響，故需要對 2～18 Hz 的壓力信號輸入的幅值進行補償以滿足后續實驗要求。
圖2
振蕩單元壓力的頻率響應特性包絡曲線
Figure2.
Envelope curves of pressure response of oscillation unit with different frequencies
1.3 偽隨機信號的產生
偽隨機信號不僅具有近似于白噪聲的性質，而且工程上容易實現，能夠保證較好的系統辨識精度，因此是一種普遍采用的輸入信號。常見的是偽隨機二進制序列，它是一種廣泛應用的隨機變量，只有 0 或 1 兩種邏輯狀態的偽隨機序列。但在 FOT 領域，采用的是偽隨機多正弦組合信號[12]。
一個頻帶有限的偽隨機信號可以通過不同頻率成分的正弦信號疊加得到。數學方面可以寫成傅里葉級數或者三角多項式的形式：
 |
其中 T 是 U 的周期，決定了頻率分辨率為
，
是分量編號矢量（mk 是正數），
表示向量的轉置，
是分量的振幅矢量 [
]，
是分量的相位矢量。對于實際信號，T 和 M 是由要求的帶寬和頻率分辨率決定的。D 的組成是根據預期設置激勵信號的頻譜，通過確定 U 的總能量
得到[12]。
我們基于 FOT 研制的振蕩肺功能儀，通過上位機產生數字信號即偽隨機多正弦組合，通過 D/A 轉換器將數字信號轉變為模擬信號輸出，經功率放大器驅動揚聲器產生高頻振蕩的偽隨機多正弦組合信號。
多正弦疊加的偽隨機信號存在相位的組合問題，可能會導致信號的周期性幅值疊加[17]。本文對 4～40 Hz 頻段的信號進行了研究。若不提前對相位的組合進行選擇，會出現較為極端的現象。為體現這一現象，選用兩個典型的多正弦信號組合，兩者每個頻點幅值相同，其中一個信號的相位完全相同，另一個信號的相位均勻變化。如圖 3 所示，為 4～40 Hz（間隔 2 Hz）范圍內，有 19 個頻點、幅值均為 300 mV 的偽隨機多正弦組合信號的頻域、時域圖。其中上一行是每個頻點幅值均相同、相位均為 0 的偽隨機多正弦信號組合的頻域和時域圖；下一行是每個頻點幅值均相同、相位為 0～2π 等間隔取值組合而成的偽隨機多正弦信號組合的頻域和時域圖。對比兩行圖片，具有相同頻率成分但不同相位特性的多正弦信號在時域表現出巨大的差異性，上圖中幅度周期性疊加最為明顯，極易超出揚聲器的響應范圍，振蕩單元的非線性失真最為嚴重。另外，FOT 要求偽隨機多正弦組合激勵信號的時域波形應盡量平坦，以減小受試者正常呼吸感到的不適。因此，采用相位優化使輸出的偽隨機信號在時域能量（幅度）分布均勻變得尤為重要。
圖3
不同相位組合對應的偽隨機信號頻域、時域圖
Figure3.
Frequency and time domain of pseudorandom signals with different phase combination
2 方法
本文嘗試了兩種方法：基于計算機仿真產生隨機數的方法和時-頻域交換迭代優化算法。為衡量優化效果的好壞，本研究采用波峰因數來量化信號在時域的能量集中程度。定義 f（t），是一個周期性多正弦諧波的組合，它的離散傅里葉頻譜為
 |
波峰因數 Kr 的定義為
 |
，它與信號包含的總能量有關。M+ 和 M– 分別表示 f（t）最大的正值和最小的負值。
Kr 越小，表明信號在時域的能量分布越均勻，信號優化算法效果越好。
2.1 基于隨機相位組合的優化算法
借鑒 Monte Carlo 方法的思想利用隨機數進行數值模擬，利用 MATLAB 軟件 rand 函數乘以 2π 隨機產生 19 個 0～2π 的相位組合（分別對應于 4～40 Hz 間隔 2 Hz 的相位），多次循環選擇輸出 Kr 值最小的相位組合，作為優化結果。本研究考慮到程序運算時間、與另一種優化算法效果的對比以及最終波形的優化效果是否滿足實驗要求，將隨機數產生循環次數設置為 150 000 次。
2.2 基于時-頻域交換的優化算法[14 , 17 -18 ]
時域-頻域交換算法具體實現過程如圖 4 所示。首先輸入信號初始的相位、頻率和幅值，將其進行組合即可得到時域波形即進行了離散傅里葉逆變換；接著，在時域上設定一個削減幅度（閾值），將真實的信號中幅值超出這一閾值的削減到該閾值；然后對削減之后的時域信號進行離散傅里葉變換得到頻率對應的相位組合；最后將相位組合重新代入到最初的幅值和頻率重復上述過程，直至 Kr 不再減少后退出循環得到相應的偽隨機多正弦組合信號的相位組合。其中，削減的幅度應設置在 5%～25% 之間，低于 5% 將會導致算法收斂速度很慢，高于 25% 將會出現嚴重的削減導致頻譜的振幅發生改變，算法發散。本文利用 MATLAB 實現“時-頻域交換算法”，削減幅度為 15%，設置循環次數上限亦為 150 000 次。
圖4
時域-頻域交換算法流程圖
Figure4.
Flow chart of time-frequency domain swapping algorithm
3 結果
3.1 振蕩單元幅頻特性及信號補償
根據歐洲呼吸協會于 2003 年對 FOT 做出的相應技術規范，在振蕩肺功能儀開路時測量，振蕩信號在管路開口處產生的振蕩壓力信號的峰峰值應在 1～3 cm H2O 之間[2]，本實驗室自行研制的振蕩肺功能儀振蕩源經測試可以滿足上述要求，以開路測量壓力為 2 cm H2O 作為后續研究標準。圖 5a、b 中藍色線分別為不考慮揚聲器低頻響應特性的上位機輸入信號和管路開口處系統壓力響應，可以看出圖 5b 中 4～18 Hz 的壓力幅值低于 2 cm H2O 的要求，20～40 Hz 的壓力響應幅值平坦并均達到了 2 cm H2O，因此需要針對 4～18 Hz 的輸入信號幅值進行補償。根據前期對振蕩源頻率響應特性的測試，對 4～18 Hz 的輸入信號幅值進行修正，得到了修正后的上位機輸入信號，如圖 5a 中紅色線所示；對應的單頻測試系統壓力響應如圖 5b 中紅色線所示，經幅值補償 4～40 Hz 各個測試頻率的壓力響應均滿足后續測試的 2 cm H2O 要求。
圖5
振蕩肺功能儀系統輸入信號和響應
a. 上位機輸入未補償和補償后 4～40 Hz（間隔 2 Hz）單頻信號；b. 上位機對輸入信號的系統壓力響應
Figure5. Input signals and response of oscillatory lung function systema. input signals of 4–40 Hz (2 Hz step) without and with compensation; b. system pressure response to the input signals
3.2 偽隨機信號優化結果
根據 3.1 得到的上位機輸入信號，若偽隨機多正弦組合的相位不進行優化，會導致偽隨機多正弦組合信號的時域波形由于周期性疊加幅值超限使得硬件設施（振蕩單元）能量過充而失真。本文基于隨機相位結合的優化算法、時-頻域交換優化算法對偽隨機多正弦組合信號進行了波形優化，波峰因數 Kr 作為最后評價優化效果的指標。兩種算法的優化對象均為零相位組合的多正弦信號（Kr=3.95），設置循環次數均為 150 000 次，得到最優的 Kr 值，其中隨機相位結合的優化算法 Kr=1.69，時-頻域交換優化算法 Kr=1.53。兩種優化算法都能夠有效降低信號的 Kr 值，從而使信號時域能量分布更均勻。相比而言，時-頻域交換算法的性能略優于隨機相位組合算法。圖 6a 為幅值疊加效應最為明顯的相位為零的偽隨機信號，可以看到幅值周期性疊加極為明顯，最大幅值接近 5 V，已達到了振蕩源的工作飽和區域；圖 6b 為隨機相位結合優化后的結果，可以看出相對于圖 6a，周期性的幅值疊加減小，幅值在時域分配較均勻，幅值降低到 2 V；圖 6c 為時-頻域交換算法優化后的波形，對比于圖 6b，時-頻域交換優化算法優化后的波形比隨機相位結合的優化更均勻且 Kr 更小。
圖6
未優化的偽隨機信號及其兩種算法優化后結果
a. 4～40 Hz 相位均為零多正弦序列偽隨機信號；b. 隨機相位組合優化后的偽隨機信號；c. 時-頻域交換算法優化后的偽隨機信號
Figure6. Pseudorandom signals without or with optimization using two algorithmsa. the sinusoidal sequence of pseudorandom signals with zero phase of 4–40 Hz; b. optimized signals with random phase combination; c. optimized signals with time-frequency domain swapping algorithm
利用上述兩種方法優化相位后的偽隨機信號，經過上位機輸入可以得到經過振蕩單元的壓力流量信號。理想情況下，輸出的壓力信號幅頻曲線應該為 4～40 Hz 范圍內幅值為 1 cm H2O（峰峰值 2 cm H2O）的平坦信號。圖 7a 為低頻信號幅值未修正的上位機輸入信號和經揚聲器的輸出壓力信號的幅值，可以看到在低頻段 4～8 Hz 的幅值小于 1 cm H2O，不能滿足激勵信號峰峰值約為 2 cm H2O 的性能要求；圖 7b 為低頻信號幅值修正后經隨機相位組合優化算法得到的輸出信號，可以看到 4～8 Hz 頻段幅值相比于圖 7a 得到了提升，但仍有部分頻率幅值不能滿足性能要求；圖 7c 為低頻信號幅值修正后經時-頻域交換算法優化得到的輸出信號，產生了滿足性能要求的幅值接近或者大于 1 cm H2O 的激勵信號。因此，對比圖 7a、7b 和 7c，低頻信號的幅值修正后，時-頻域交換算法優化得到的激勵信號最為理想且滿足了性能要求。
圖7
不同輸入信號以及經不同算法優化后的輸出信號頻譜圖
a. 未修正幅值的輸入信號和輸出壓力信號；b. 幅值修正后輸入信號和經隨機相位組合算法優化的輸出信號；c. 幅值修正后輸入信號和經時-頻域交換算法優化的輸出信號
Figure7. The spectrums of different input signals and their output signals without or with different optimization algorithmsa. input signals without amending amplitude and output signals; b. input signals with amending amplitude and output signals optimized with random phase combination algorithm; c. input signals with amending amplitude and output signals optimized with time-frequency domain swapping algorithm
4 討論
FOT 是基于系統辨識原理對肺功能進行診斷的一種新技術，本文對該技術目前常用的激勵信號之一——偽隨機多正弦組合信號的產生和優化進行了研究，采用了隨機相位組合的優化算法和時-頻域交換的優化算法。零相位組合的信號經兩種算法優化后波峰因數 Kr 由 3.95 分別變為 1.69 和 1.53，提示兩種算法都能夠優化相位組合；相對而言，時-頻域交換優化算法效果更好。最終利用時-頻域交換優化算法得到的激勵信號，在 4～40 Hz 頻率范圍內的峰峰值均滿足接近 2 cm H2O 的要求。
強迫振蕩肺功能實現的關鍵技術有兩大部分：一是對振蕩單元的特性辨識產生滿足測量要求的激勵信號；二是在不同頻率信號激勵下呼吸系統不同部位的力學特性辨識。本文解決了關鍵技術實現的第一部分。呼吸系統的正常生理過程體現和病理結構改變在檢測頻率范圍小于 10 Hz 時最為敏感，因而對呼吸系統低頻段的響應測試具有重要的意義，而振蕩單元揚聲器部分的低頻段特性以及聲波在管路傳播時低頻段容易在開口處泄漏更多能量，使得振蕩單元產生的激勵信號在低頻段幅值不能滿足要求，因而需要對低頻段幅值進行補償。
本文采用隨機相位組合的優化算法和時-頻域交換的優化算法對偽隨機信號相位進行了優化。由于時-頻域交換優化算法的原理是對波形的迭代修正，與隨機相位組合的優化算法相比，是一種直接的優化算法，運算時間短（產生相同結果）并且優化效果理想。隨機相位組合的優化算法在性能上略低于時-頻域交換的優化算法，可能原因是偽隨機信號設置以及循環次數限制，沒有產生并篩選出足夠精細的相位組合。從理論上講，計算機可以產生一種相位組合，性能逼近甚至優于時-頻域交換的優化算法，但運算時間會比較長。因此，時-頻域交換的優化算法是一種相對優選的算法。
本文圖 7b、c 輸出的激勵信號部分頻率對應的幅值小于 1 cm H2O，部分高于 1 cm H2O，導致這一現象的原因可能有兩方面。一是本文的激勵信號輸入幅值是根據單頻正弦（4～40 Hz）信號通過振蕩單元進行測試獲得，而實際輸入為偽隨機多正弦組合信號。由于揚聲器本身的特性，對兩種不同類型輸入信號的響應存在一定差異，可能使最終振蕩單元產生的激勵信號的幅值并不完全相同。二是輸出信號中涵蓋了本研究頻域范圍之外的 40～100 Hz 的諧波成分，考慮可能是由于產生 4～40 Hz 激勵信號時同步產生了諧波干擾，導致在 4～40 Hz 的部分頻段幅值的提高和降低。
目前 FOT 肺功能檢測的另一種重要應用類型是脈沖振蕩法（impulse oscillometry system，IOS），其激勵信號是周期三角脈沖信號。本文通過對偽隨機多正弦組合激勵信號的研究，將為后續周期脈沖激勵信號的研究提供非常有利的幫助[2, 19]，同時為我們后續研究的重點——呼吸系統力學特性的辨識奠定扎實基礎。
引言
強迫振蕩技術（forced oscillation technique，FOT）是一種主動測量呼吸系統力學特性的技術，其核心思想是將人體的呼吸系統視為一個“黑箱”，利用系統辨識的原理，通過外界施加不同類型的壓力振蕩，測量呼吸系統的響應，即壓力、流量，進而辨識出呼吸系統阻抗[1-3]。FOT 是一種無創技術，無需受試者主動配合，測量過程中只需要受試者平靜地進行自主呼吸，即可非常靈敏地測試出呼吸系統力學特性參數，特別適合于老年人、兒童使用[4]；它可以分辨吸氣相和呼氣相阻力變化，對于辨別慢性阻塞性肺疾病和哮喘具有重要意義[5-8]；另外，它還有助于對吸煙患者呼吸系統病變的辨識，可以幫助臨床早期診斷由吸煙導致的呼吸功能改變[9]。因而，FOT 在呼吸系統疾病診斷和指導用藥方面具有廣闊的應用前景[3]。
FOT 通過壓力振蕩源施加連續的激勵信號以檢測系統的響應[3]。FOT 常采用的激勵信號有三種，分別為單頻正弦信號、偽隨機信號（多正弦組合）和周期脈沖信號。單頻正弦信號因容易實現、信噪比高以及對硬件部分要求較低，在 FOT 的技術早期得到了廣泛使用。但單頻正弦信號若要實現多個頻率點的呼吸阻抗測量，需要依次掃頻，存在檢測時間偏長和時間分辨率較低的問題。后來，相繼出現了偽隨機多正弦組合和周期性脈沖等激勵信號，它們能夠在較短的時間內產生呼吸系統辨識所需的持續激勵信號，使呼吸系統在多個頻率點被充分激勵，因此能夠有效縮短檢測時間，并具有較高的時間分辨率，但該方法對硬件電路和算法提出了更高要求[3, 10-12]。本文著重介紹偽隨機多正弦組合的相關研究，后續文章會介紹周期性脈沖信號的研究進展。
FOT 的偽隨機多正弦組合激勵信號本質上是多個正弦信號的疊加，存在相位組合問題。多正弦信號組合在頻域特性滿足要求的情況下，不同的相位組合會產生時域特性不同的信號，若相位均相同，將出現最為極端的情況，時域信號周期性疊加，信號幅值太大，極易導致系統飽和失真[3, 12]。另一方面，呼吸系統對于低頻激勵信號響應更為敏感，可以反映更多正常生理過程和病理結構改變的內容，因而臨床更為關注[7, 13-14]。而 FOT 振蕩單元產生的激勵信號的幅頻特性曲線存在低頻響應幅度過低的問題，因此需要對輸入激勵信號的低頻段的幅值進行補償。本文針對強迫振蕩技術肺功能儀研制過程中遇到的偽隨機多正弦組合信號產生、信號低頻段幅度補償以及相位優化問題進行了深入研究。
1 原理
1.1 振蕩肺功能基本原理
FOT 的經典裝置示意圖如圖 1 所示，壓力振蕩單元通常為揚聲器，產生單頻或者復合頻振蕩的壓力信號，加載在呼吸道上，受試者在測試過程中進行自主呼吸，同步采集自主呼吸過程中的壓力信號 Pao 和流量信號 V′[1-2]。運用系統辨識原理，將呼吸系統等效為一個“黑箱”，輸入信號為壓力，輸出信號為流量，去除自主呼吸產生的干擾，壓力和流量信號分別做傅里葉變換后相比得到呼吸系統的總阻抗 Zrs。
圖1
強迫振蕩肺功能系統示意圖
Figure1.
Schematic arrangement of the forced oscillatory pulmonary function system
1.2 振蕩單元的幅頻特性及補償
基于上述經典裝置示意圖，實驗室設計了振蕩肺功能儀，其振蕩單元包含揚聲器和功率放大器。其中揚聲器采用 JBL 的 145 型超低音喇叭，額定功率 200 W，阻抗 8 Ω，靈敏度為 92 dB/W。本研究首先測試了振蕩單元整體的幅頻特性，利用上位機軟件（LabVIEW2013）分別產生幅值為 300 mV、2～50 Hz（2 Hz 等間隔）正弦掃頻的數字信號，經過 NI 公司的 DAQ 設備 USB-6366 的 D/A 轉換成模擬信號，經功率放大器后驅動揚聲器振蕩[15]。輸入幅值為 300 mV，是經過前期振蕩單元性能測試以產生開路測量壓力為 2 cm H2O 的最佳輸入幅值。裝置輸出端空載時測量得到系統響應的壓力信號，如圖 2 所示。從圖中可以明顯得出，由于揚聲器本身的頻響特性以及低頻聲波傳播時在管路開口處能量泄漏要高于高頻聲波[16]，因此在輸入信號幅值相同的情況下，系統輸出壓力信號幅度在 2～18 Hz 頻率范圍內明顯衰減，會對我們后續研究呼吸系統在低頻段的激勵響應帶來很大影響，故需要對 2～18 Hz 的壓力信號輸入的幅值進行補償以滿足后續實驗要求。
圖2
振蕩單元壓力的頻率響應特性包絡曲線
Figure2.
Envelope curves of pressure response of oscillation unit with different frequencies
1.3 偽隨機信號的產生
偽隨機信號不僅具有近似于白噪聲的性質，而且工程上容易實現，能夠保證較好的系統辨識精度，因此是一種普遍采用的輸入信號。常見的是偽隨機二進制序列，它是一種廣泛應用的隨機變量，只有 0 或 1 兩種邏輯狀態的偽隨機序列。但在 FOT 領域，采用的是偽隨機多正弦組合信號[12]。
一個頻帶有限的偽隨機信號可以通過不同頻率成分的正弦信號疊加得到。數學方面可以寫成傅里葉級數或者三角多項式的形式：
|
其中 T 是 U 的周期，決定了頻率分辨率為
，
是分量編號矢量（mk 是正數），
表示向量的轉置，
是分量的振幅矢量 [
]，
是分量的相位矢量。對于實際信號，T 和 M 是由要求的帶寬和頻率分辨率決定的。D 的組成是根據預期設置激勵信號的頻譜，通過確定 U 的總能量
得到[12]。
我們基于 FOT 研制的振蕩肺功能儀，通過上位機產生數字信號即偽隨機多正弦組合，通過 D/A 轉換器將數字信號轉變為模擬信號輸出，經功率放大器驅動揚聲器產生高頻振蕩的偽隨機多正弦組合信號。
多正弦疊加的偽隨機信號存在相位的組合問題，可能會導致信號的周期性幅值疊加[17]。本文對 4～40 Hz 頻段的信號進行了研究。若不提前對相位的組合進行選擇，會出現較為極端的現象。為體現這一現象，選用兩個典型的多正弦信號組合，兩者每個頻點幅值相同，其中一個信號的相位完全相同，另一個信號的相位均勻變化。如圖 3 所示，為 4～40 Hz（間隔 2 Hz）范圍內，有 19 個頻點、幅值均為 300 mV 的偽隨機多正弦組合信號的頻域、時域圖。其中上一行是每個頻點幅值均相同、相位均為 0 的偽隨機多正弦信號組合的頻域和時域圖；下一行是每個頻點幅值均相同、相位為 0～2π 等間隔取值組合而成的偽隨機多正弦信號組合的頻域和時域圖。對比兩行圖片，具有相同頻率成分但不同相位特性的多正弦信號在時域表現出巨大的差異性，上圖中幅度周期性疊加最為明顯，極易超出揚聲器的響應范圍，振蕩單元的非線性失真最為嚴重。另外，FOT 要求偽隨機多正弦組合激勵信號的時域波形應盡量平坦，以減小受試者正常呼吸感到的不適。因此，采用相位優化使輸出的偽隨機信號在時域能量（幅度）分布均勻變得尤為重要。
圖3
不同相位組合對應的偽隨機信號頻域、時域圖
Figure3.
Frequency and time domain of pseudorandom signals with different phase combination
2 方法
本文嘗試了兩種方法：基于計算機仿真產生隨機數的方法和時-頻域交換迭代優化算法。為衡量優化效果的好壞，本研究采用波峰因數來量化信號在時域的能量集中程度。定義 f（t），是一個周期性多正弦諧波的組合，它的離散傅里葉頻譜為
|
波峰因數 Kr 的定義為
|
，它與信號包含的總能量有關。M+ 和 M– 分別表示 f（t）最大的正值和最小的負值。
Kr 越小，表明信號在時域的能量分布越均勻，信號優化算法效果越好。
2.1 基于隨機相位組合的優化算法
借鑒 Monte Carlo 方法的思想利用隨機數進行數值模擬，利用 MATLAB 軟件 rand 函數乘以 2π 隨機產生 19 個 0～2π 的相位組合（分別對應于 4～40 Hz 間隔 2 Hz 的相位），多次循環選擇輸出 Kr 值最小的相位組合，作為優化結果。本研究考慮到程序運算時間、與另一種優化算法效果的對比以及最終波形的優化效果是否滿足實驗要求，將隨機數產生循環次數設置為 150 000 次。
2.2 基于時-頻域交換的優化算法[14 , 17 -18 ]
時域-頻域交換算法具體實現過程如圖 4 所示。首先輸入信號初始的相位、頻率和幅值，將其進行組合即可得到時域波形即進行了離散傅里葉逆變換；接著，在時域上設定一個削減幅度（閾值），將真實的信號中幅值超出這一閾值的削減到該閾值；然后對削減之后的時域信號進行離散傅里葉變換得到頻率對應的相位組合；最后將相位組合重新代入到最初的幅值和頻率重復上述過程，直至 Kr 不再減少后退出循環得到相應的偽隨機多正弦組合信號的相位組合。其中，削減的幅度應設置在 5%～25% 之間，低于 5% 將會導致算法收斂速度很慢，高于 25% 將會出現嚴重的削減導致頻譜的振幅發生改變，算法發散。本文利用 MATLAB 實現“時-頻域交換算法”，削減幅度為 15%，設置循環次數上限亦為 150 000 次。
圖4
時域-頻域交換算法流程圖
Figure4.
Flow chart of time-frequency domain swapping algorithm
3 結果
3.1 振蕩單元幅頻特性及信號補償
根據歐洲呼吸協會于 2003 年對 FOT 做出的相應技術規范，在振蕩肺功能儀開路時測量，振蕩信號在管路開口處產生的振蕩壓力信號的峰峰值應在 1～3 cm H2O 之間[2]，本實驗室自行研制的振蕩肺功能儀振蕩源經測試可以滿足上述要求，以開路測量壓力為 2 cm H2O 作為后續研究標準。圖 5a、b 中藍色線分別為不考慮揚聲器低頻響應特性的上位機輸入信號和管路開口處系統壓力響應，可以看出圖 5b 中 4～18 Hz 的壓力幅值低于 2 cm H2O 的要求，20～40 Hz 的壓力響應幅值平坦并均達到了 2 cm H2O，因此需要針對 4～18 Hz 的輸入信號幅值進行補償。根據前期對振蕩源頻率響應特性的測試，對 4～18 Hz 的輸入信號幅值進行修正，得到了修正后的上位機輸入信號，如圖 5a 中紅色線所示；對應的單頻測試系統壓力響應如圖 5b 中紅色線所示，經幅值補償 4～40 Hz 各個測試頻率的壓力響應均滿足后續測試的 2 cm H2O 要求。
圖5
振蕩肺功能儀系統輸入信號和響應
a. 上位機輸入未補償和補償后 4～40 Hz（間隔 2 Hz）單頻信號；b. 上位機對輸入信號的系統壓力響應
Figure5. Input signals and response of oscillatory lung function systema. input signals of 4–40 Hz (2 Hz step) without and with compensation; b. system pressure response to the input signals
3.2 偽隨機信號優化結果
根據 3.1 得到的上位機輸入信號，若偽隨機多正弦組合的相位不進行優化，會導致偽隨機多正弦組合信號的時域波形由于周期性疊加幅值超限使得硬件設施（振蕩單元）能量過充而失真。本文基于隨機相位結合的優化算法、時-頻域交換優化算法對偽隨機多正弦組合信號進行了波形優化，波峰因數 Kr 作為最后評價優化效果的指標。兩種算法的優化對象均為零相位組合的多正弦信號（Kr=3.95），設置循環次數均為 150 000 次，得到最優的 Kr 值，其中隨機相位結合的優化算法 Kr=1.69，時-頻域交換優化算法 Kr=1.53。兩種優化算法都能夠有效降低信號的 Kr 值，從而使信號時域能量分布更均勻。相比而言，時-頻域交換算法的性能略優于隨機相位組合算法。圖 6a 為幅值疊加效應最為明顯的相位為零的偽隨機信號，可以看到幅值周期性疊加極為明顯，最大幅值接近 5 V，已達到了振蕩源的工作飽和區域；圖 6b 為隨機相位結合優化后的結果，可以看出相對于圖 6a，周期性的幅值疊加減小，幅值在時域分配較均勻，幅值降低到 2 V；圖 6c 為時-頻域交換算法優化后的波形，對比于圖 6b，時-頻域交換優化算法優化后的波形比隨機相位結合的優化更均勻且 Kr 更小。
圖6
未優化的偽隨機信號及其兩種算法優化后結果
a. 4～40 Hz 相位均為零多正弦序列偽隨機信號；b. 隨機相位組合優化后的偽隨機信號；c. 時-頻域交換算法優化后的偽隨機信號
Figure6. Pseudorandom signals without or with optimization using two algorithmsa. the sinusoidal sequence of pseudorandom signals with zero phase of 4–40 Hz; b. optimized signals with random phase combination; c. optimized signals with time-frequency domain swapping algorithm
利用上述兩種方法優化相位后的偽隨機信號，經過上位機輸入可以得到經過振蕩單元的壓力流量信號。理想情況下，輸出的壓力信號幅頻曲線應該為 4～40 Hz 范圍內幅值為 1 cm H2O（峰峰值 2 cm H2O）的平坦信號。圖 7a 為低頻信號幅值未修正的上位機輸入信號和經揚聲器的輸出壓力信號的幅值，可以看到在低頻段 4～8 Hz 的幅值小于 1 cm H2O，不能滿足激勵信號峰峰值約為 2 cm H2O 的性能要求；圖 7b 為低頻信號幅值修正后經隨機相位組合優化算法得到的輸出信號，可以看到 4～8 Hz 頻段幅值相比于圖 7a 得到了提升，但仍有部分頻率幅值不能滿足性能要求；圖 7c 為低頻信號幅值修正后經時-頻域交換算法優化得到的輸出信號，產生了滿足性能要求的幅值接近或者大于 1 cm H2O 的激勵信號。因此，對比圖 7a、7b 和 7c，低頻信號的幅值修正后，時-頻域交換算法優化得到的激勵信號最為理想且滿足了性能要求。
圖7
不同輸入信號以及經不同算法優化后的輸出信號頻譜圖
a. 未修正幅值的輸入信號和輸出壓力信號；b. 幅值修正后輸入信號和經隨機相位組合算法優化的輸出信號；c. 幅值修正后輸入信號和經時-頻域交換算法優化的輸出信號
Figure7. The spectrums of different input signals and their output signals without or with different optimization algorithmsa. input signals without amending amplitude and output signals; b. input signals with amending amplitude and output signals optimized with random phase combination algorithm; c. input signals with amending amplitude and output signals optimized with time-frequency domain swapping algorithm
4 討論
FOT 是基于系統辨識原理對肺功能進行診斷的一種新技術，本文對該技術目前常用的激勵信號之一——偽隨機多正弦組合信號的產生和優化進行了研究，采用了隨機相位組合的優化算法和時-頻域交換的優化算法。零相位組合的信號經兩種算法優化后波峰因數 Kr 由 3.95 分別變為 1.69 和 1.53，提示兩種算法都能夠優化相位組合；相對而言，時-頻域交換優化算法效果更好。最終利用時-頻域交換優化算法得到的激勵信號，在 4～40 Hz 頻率范圍內的峰峰值均滿足接近 2 cm H2O 的要求。
強迫振蕩肺功能實現的關鍵技術有兩大部分：一是對振蕩單元的特性辨識產生滿足測量要求的激勵信號；二是在不同頻率信號激勵下呼吸系統不同部位的力學特性辨識。本文解決了關鍵技術實現的第一部分。呼吸系統的正常生理過程體現和病理結構改變在檢測頻率范圍小于 10 Hz 時最為敏感，因而對呼吸系統低頻段的響應測試具有重要的意義，而振蕩單元揚聲器部分的低頻段特性以及聲波在管路傳播時低頻段容易在開口處泄漏更多能量，使得振蕩單元產生的激勵信號在低頻段幅值不能滿足要求，因而需要對低頻段幅值進行補償。
本文采用隨機相位組合的優化算法和時-頻域交換的優化算法對偽隨機信號相位進行了優化。由于時-頻域交換優化算法的原理是對波形的迭代修正，與隨機相位組合的優化算法相比，是一種直接的優化算法，運算時間短（產生相同結果）并且優化效果理想。隨機相位組合的優化算法在性能上略低于時-頻域交換的優化算法，可能原因是偽隨機信號設置以及循環次數限制，沒有產生并篩選出足夠精細的相位組合。從理論上講，計算機可以產生一種相位組合，性能逼近甚至優于時-頻域交換的優化算法，但運算時間會比較長。因此，時-頻域交換的優化算法是一種相對優選的算法。
本文圖 7b、c 輸出的激勵信號部分頻率對應的幅值小于 1 cm H2O，部分高于 1 cm H2O，導致這一現象的原因可能有兩方面。一是本文的激勵信號輸入幅值是根據單頻正弦（4～40 Hz）信號通過振蕩單元進行測試獲得，而實際輸入為偽隨機多正弦組合信號。由于揚聲器本身的特性，對兩種不同類型輸入信號的響應存在一定差異，可能使最終振蕩單元產生的激勵信號的幅值并不完全相同。二是輸出信號中涵蓋了本研究頻域范圍之外的 40～100 Hz 的諧波成分，考慮可能是由于產生 4～40 Hz 激勵信號時同步產生了諧波干擾，導致在 4～40 Hz 的部分頻段幅值的提高和降低。
目前 FOT 肺功能檢測的另一種重要應用類型是脈沖振蕩法（impulse oscillometry system，IOS），其激勵信號是周期三角脈沖信號。本文通過對偽隨機多正弦組合激勵信號的研究，將為后續周期脈沖激勵信號的研究提供非常有利的幫助[2, 19]，同時為我們后續研究的重點——呼吸系統力學特性的辨識奠定扎實基礎。
