用Matlab 實現隨機信號的實時處理

魯曉東，紀玉川

（1.浙江海洋大學東海科學技術學院，浙江舟山 316000；2.舟山市海大科學技術研究院，浙江舟山 316000）

Matlab 由于其強大的計算和繪圖功能，在工程領域中常被用于各種設計方案、處理算法的仿真[1-2]，以簡化計算過程，提高工作效率。同時Matlab 對硬件的支持[3]也是可圈可點的，例如Matlab 的數據采集工具箱(DAQ)就包含了市面上的各種數據采集卡驅動。由于它是通過自身的MEX 工具包調用操作系統的底層API 去驅動計算機的外設，對用戶來說可以不必深究硬件如何工作，僅調用工具包的相關函數就能獲取外部傳感器的數據。所以Matlab 也非常適合于信號的實時處理[4-5]，例如在嵌入式開發前期，硬件系統并未建立完善，當需知某種處理方法在設計的系統中是否有效，可以直接借助PC 機上的Matlab 對該算法進行驗證，因為系統處理的是真實信號，只是算法在不同的平臺進行實現而已，所以仿真的結果幾近于實際效果。

1 Matlab信號處理實現原理

1.1 DAQ工作原理

DAQ(Data AcQuisition)是Matlab的數據采集工具箱，該工具箱由3個部分組成:功能函數.m文件、DAQ引擎和硬件驅動程序。為了方便用戶的使用，Matlab對其內部的DAQ 底層驅動作了函數的封裝，通過調用這些函數便可完成設備數據的采集。由于Matlab包含當前市面上多數的硬件驅動，因此，其DAQ 有較好的兼容性，而用戶使用時感覺不到驅動的存在，只須調用函數就可直接控制設備[6]。例如要完成的是對音頻實時信號的處理，使用系統集成的聲卡作為DAQ 設備，Matlab 對聲卡的功能封裝如圖1 所示，Matlab 通過實例化一個audio Device Reader 對象，設置對象的屬性并使用對象提供的方法驅動聲卡工作以獲得數據。Matlab應用的基本方法可以描述為：

圖1 音頻數據采集卡功能結構

1）創建音頻采集卡DAQ 設備對象

MyDaq=audio Device Reader，

其中，audio Device Reader 是Matlab 實例化系統音頻卡的功能模塊，MyDaq 即為設備對象的句柄。

2）設置DAQ 對象的屬性參數

屬性參數主要包括采集卡工作時數據幀的幀緩沖大小、采樣率、通道數、采樣點數據類型等，屬性值可以根據系統工作的需要進行設置。例如：

MyDaq.Samples PerFrame=44100//采樣率MyDaq.Num Channels=2//通道數為2

MyDaq.BitDepth=′16-bit integer′//16 位整數

3）數據的采集控制

使用audio Device Reader對象的record 或step 方法，就可以讀取音頻采集卡中緩沖區里的數據。緩沖區一般存放一幀的數據，軟件讀取的間隔一般都遠慢于硬件的采樣。所以會產生溢出，讀取的是緩沖區最后一幀數據。

[x，numOverrun]=record(MyDaq)//record 在MyDaq讀取數據并存放在x，返回幀溢出數numOverrun。

4）釋放內存

release(MyDaq)；//釋放該采集卡MyDaq 所占用的內存資源。

1.2 Matlab信號處理原理

以一個在實際環境中常見的反射信號的處理為例。通過Matlab 提供的信號處理工具計算其反射混響的時間參數，對其濾波并消除該反射信號。

不妨設x1(n)為信號源信號，經過某終端產生一系列衰減系數為αk延遲nk的反射信號，它與原信號疊加后的波形為：

x(n)可看成x1(n)與沖激串∑δ(n-nk)的卷積，令h(n)=δ(n)+α1δ(n-n1)+…+αkδ(n-nk)，則式（1）可以改寫為:

h(n)實際上就是信道的沖激響應。若想從x(n)中分離出x1(n)，就必須知道h(n) 的相關信息，才可以對式（1）作反卷積或解差分方程。由于混響信號的特點是自身信號的延時疊加，所以可以借助信號的倒譜來觀測h的延時參數nk。先把式（2）變換到z域：

兩邊取對數得到：

所以對于延遲nk個序列的反射信號一定在倒譜上的第nk個序列值有一個峰值，因為此時的卷積值取得極值。通過尋找數值的峰值便可以確定延遲的參數nk，從而確定了式（1）所表達的差分方程。設僅存在一個反射波混合信號，利用Matlab 提供的濾波器函數filter（b，a，x）求解該差分方程，可以解算出原始信號x(n) 。基本步驟如下：

1）用倒譜函數rceps 計算出混合信號的倒譜序列。

cepsig=rceps(xecho)；//xecho為混合時域信號，cepsig 為倒譜序列。

2）用findpeaks()檢索倒譜的峰值，確定反射信號的延遲參數delay。

[px，nk]=findpeaks(cepsig，′ Threshold′，0.01，′ Min PeakDistance′，10)；

delay=nk(2)-1；

函數findpeaks()將在倒譜cepsig 中尋找峰值，并返回峰值的序列位置nk以及相應峰值px。參數MinPeakDistance 是設置檢峰的最小時間尺度，用于過濾噪聲干擾。Threshold 則是設置各峰值必須與相鄰數據點的差值。這兩個參數會影響檢測的結果，操作中可以根據實際進行調整。倒譜序列的峰值在第二個即為反射序列出現的位置，因第一個位置nk(1)一定為1，延時為0，所以第二個位置nk(2)的延時要進行減1 調整。

3）用filter()函數解算差分方程并求得原信號。根據filter(b，a，x)函數的時域表達式[8]:

可確定存在一個反射信號的混響信號x(n)=x1(n)+αkx1(n-nk)的a，b參數。用Matlab描述為：

2 信號實時處理實現示例

2.1 信號處理的界面

為提高編程效率，采用Matlab GUI 方式[9]建立信號處理的控制界面如圖3 所示。在面板上布置4 個AXES（圖形坐標系）對象即程序中的signalraw、fsignalraw、cepstrum、firsignal，分別用于顯示信道信號、信號的頻譜、倒譜、以及濾波后的信號。右邊是濾波器控制面板、信號采集控制的按鈕等，用于實現對各控件函數的回調。

2.2 處理控制的實現

處理控制主要是完成信號的實時處理和數據的可視化處理。在建立采集卡對象后，便可以操作該對象提供的DAQ函數，完成數據的采集、處理與顯示。

1）初始化DAQ 對象和圖形對象

以下是數據可視化圖形句柄初始化：

2）實時處理

由于Matlab 是一種單線程的應用程序，因此控制腳本使用了循環結構，用循環結構中的各條件變量決定信號處理的流程和方式，如圖2 所示。

圖2 信號處理流程結構

這是人機對話的控制處理，主要完成數據采集的啟動、停止以及濾波器的設置。控制的內容都是實時處理腳本中的一些參數，如啟動、停止的flag 變量，flag 為1 則while 中不斷采集數據并處理；選擇是否需要手動改變濾波系數[b，a]，設置的變量是default，當default為0則系統自動估算濾波參數，否則在面板上手動調節截止頻率，并定義處理的函數filter_set()計算調節后相應的參數。

3 實時處理結果

用一個拾音器采集一個音頻信號，并給于一定的回聲效果，時延約為0.05 s，通過聲卡的A/D 轉換，在Matlab 上得到波形如圖3（a）所示，同時可以看到其頻譜如圖3（b）所示，在頻譜中看不到回聲延遲的特征，但可以把它轉為倒譜如圖3（c）所示，可以看出在原點和0.043 處各有一個峰值，原點是信號本身，而第二個就是延時的信號。經過濾波得到信號的恢復與測試的原始信號作對比，如圖3（d）所示，效果是比較理想的。這個過程也可以通過界面上的濾波器手動設置來實現，信號處理結果也可實時顯示[10-17]。

圖3 實時處理結果

4 結束語

從上述處理過程可見，利用Matlab 進行信號的實時處理非常方便，主要表現在其對硬件良好的兼容性，使用戶不必去關注硬件的驅動，只需關心軟件的算法效果。而且它封裝了大量的工具函數，便于計算和調試，可以提高實現處理的效率。值得注意的是Matlab 是一種解釋性的語言，當有多路信號同時要求處理時，Matlab 的這種循環結構的編程方式會降低效率，實時性就不那么理想，但是其作為一種對實時處理算法的有效性檢驗手段是非常值得應用的。