林偉毅
(華南理工大學自動化科學與工程學院,廣東廣州510640)
近年來,Donoho和Candes等人提出了一種采樣與壓縮同步進行的理論——壓縮感知理論(Compressive Sensing,CS)[1],與傳統(tǒng)的“先采樣,后壓縮”不同,CS理論是“邊采樣,邊壓縮”的方法。更重要的是,CS理論指出,信號的采樣頻率可以遠遠低于信號頻譜中最高頻率的2倍,而不影響信號的精確重建,只要信號是稀疏的,以低于奈奎斯特頻率對信號進行采樣的方式,即次奈奎斯特采樣。將次奈奎斯特采樣應用于超聲波成像中,可以有效的減少采樣點和采樣頻率,而這對于縮小超聲波成像系統(tǒng)的儀器體積以及減少工作過程中所消耗的電能,都有著積極的意義。
FRI[2]模型是傳統(tǒng)的采樣理論與壓縮感知相結合的采集信號的新方案,它最早由Vetterli等人提出來。FRI模型的提出是基于以下的設想的:在自然界中,很多信號都可以分解為一系列已知形狀的短脈沖(假設脈沖數(shù)為L)。這樣的信號,可以由各個短脈沖信號的延遲時間和幅度進行完備的表示。因此在特定的時間T內,只有有限的稀疏度(2L),文獻[2]將稀疏度2L與時間T的比,定義為ROI(the Rate of Invovation),這樣的信號叫做FRI(Finite Rate of Innovation)信號。ROI是理論上,F(xiàn)RI信號能夠完美重建的最低采樣頻率,它遠遠低于奈奎斯特采樣頻率。超聲波成像系統(tǒng)的反射信號可以看成由一系列不同延遲時間和幅度的脈沖信號的疊加。因此,可以利用FRI模型中的采樣方法進行信號采樣和重建。不幸的是,在之前的相關文獻中,提出的一些對FRI信號進行采樣的方案,比較著名的像B-Spline濾波器,E-Spline濾波器[3],Intergrators濾波器[4]以及Exponential濾波器[5],對于帶噪聲的高ROI信號(L≥10)的重建效果并不理想。所以,尋找一種適合高噪聲、高ROI信號的次奈奎斯特采樣方案,對于次奈奎斯特采樣在超聲波成像中的應用至關重要,本文的目的正在于此。
假設我們得到了如下的一個持續(xù)時間為(0,T)的有限長FRI信號:
其中,h(t)是形狀已知的脈沖并且滿足關系式:
{tl,為未知的延遲時間和幅度,L為短脈沖的個數(shù)。假設我們得到了如下所示的M個傅里葉系數(shù):
其中H(ω)為h(t)的連續(xù)傅里葉變換,κ是長度為M連續(xù)整數(shù)取值區(qū)間。令H為M×M的對角矩陣,其中第k個對角元素為令V(t)為的矩陣,它的第(k,j)元素為其中向量t={t1,t2,…,tL},令向量a={a1,a2,…,aL},令向量x表示傅里葉系數(shù)向量。那么(2)式可以重新寫成矩陣形式:
我們可以選擇區(qū)間,使得矩陣H是可逆的,定義y=H-1x,則有,
矩陣V是Vandermonde矩陣,因此只要滿足條件2L≤M,V則是列滿秩的。將方程組(5)進行重寫,則有
由式(3),(4),(5)可知,一旦得到了傅里葉系數(shù)向量x,并且保證2L≤M,以及ti≠tj(i≠j),就可以根據(jù)方程組(6)解出,本文采用Annihilating Filter Method[3]進行信號重建,當然也可以使用其它的譜分析[6]。
下面將討論一下如何在時間域上對FRI信號x(t)進行采樣,以得到傅里葉系數(shù)向量x。
我們采用的是多通道混合采樣的方式,使得每個通道的采樣輸出都是傅里葉系數(shù)向量x的線性組合,以防止通道失效,傅里葉系數(shù)丟失的問題。有限長FRI信號采樣框架如圖1所示。
圖1 有限長FRI信號采樣框架Fig.1 Sampling scheme of Finite stream of pulses
假設有P條通道,由圖1可知,第i條通道的載波信號為:
其中,每一條通道中的k∈κ不盡相同,設采樣向量為c={c1,c2,…,cP},則有
令S為P×M的矩陣,它的第(i,k)元素為sik,那么(8)式可以改寫為矩陣形式:
當且僅當S列滿秩的時候,即P≥M時,S左可逆,則有:
S由我們所選擇的載波信號si(t)所決定。
我們采用方波pi(t)產生載波si(t),因為在實際電路中,方波的發(fā)生電路比較簡單。方波的數(shù)學模型如下所示:
其中,ai[n]為長度是N的系數(shù)向量(N≥M),它的元素只能取±1,p(t)為:
它的連續(xù)傅里葉變換為:
下面,我們將會介紹:以濾波的方式,從方波信號pi(t)得到實際載波信號si(t),以及推導出這種情況下,混合矩陣S的計算公式。
1)由方波信號pi(t)產生si(t)
由式(14)可以看出,pi(t)是一個周期函數(shù)。我們知道,周期函數(shù)可以由傅里葉系數(shù)如下表示:
其中,傅里葉系數(shù)由di[k]由以下公式計算得出:
可以看出,與式(7)相比,式(14)是無限相加的形式,因此,為了產生si(t),必須加一個濾波器g(t),對pi(t)進行濾波
g(t)的連續(xù)傅里葉變換必須有如下限制:
那么,可得,
由式(7)和式(17)可知(18),下面推導S時將會用到:
2)混合矩陣S的計算公式
由式(11)、(15)計算可得:
結合(19)和(20),可得:
其中k′=-(k-■M/2」)。那么混合矩陣S可以表示為:
A為P×N的矩陣,第(i,n)元素Ain=ai[n];W為N×M的矩
結合式(13)、(17),(23)可以改寫為:
1)利用超聲波工業(yè)成像系統(tǒng)采集一組一維原始信號,其中探頭的中心頻率為fc=3.423 5 MHz,系統(tǒng)的采樣頻率fs=40 MHz,成像深度為Rmax=0.16 m,在被測物體中,超聲波速C=1 540,那么可以計算出信號的持續(xù)時間為T=2×Rmax/C=2.08×10-4sec,因此可以計算出,若用傳統(tǒng)采樣方法,每次成像的數(shù)據(jù)點數(shù)目為SamplingNo_classic=T×fs=8 320,我們的仿真,只采用了一半數(shù)據(jù),因此T=T/2=1.04×10-4sec,SamplingNo_classic=SamplingNo_classic/2=4 160;
2)根據(jù)超聲波工業(yè)探頭的特性,可知h(t)為高斯脈沖是比較合適的,其中高斯參數(shù),σ=3×10-7
3)利用圖1所示的框架對原始信號進行混合積分采樣,設L=4,P=N=M=2×L×oversampling+1,由于原始信號噪聲很大,因此設oversampling=5,則:P=N=M=4,κ={-20,-19,…,0,…,19,20};
4)我們選擇方波信號作為載波,其中方波信號符合式(11)、(12)、方波的個數(shù)為41,周期為T,寬度為T/N,濾波器g(t)符合式(17),該濾波器為低通濾波器,帶寬為2π/T;
5)積分采樣后,得到向量c,根據(jù)方波信號寫出A矩陣,根據(jù)式(22)、(24)得到S矩陣;
6)利用公式(10)計算出傅里葉系數(shù)向量x,利用公式,y=H-1x,求出向量y,從而得到方程組(6);
7)利用Annihilating Filter Method求解方程(6),得到參數(shù),重建效果,如圖2所示,其中信號大小和時間都歸1了。
圖2 信號重建效果圖(虛線部分)Fig.2 Applying our method on real ultrasound imaging data.Results are shown vs.full demodulated signal
由圖2可知,延遲時間的重建十分準確(其它的波峰波谷為噪聲),而幅值的重建則有偏差,尤其是第二個高斯脈沖。但這并不會影響到超聲波工業(yè)成像的效果,因為超聲波工業(yè)成像是為了找尋缺陷的位置,而延遲時間包含了位置的全部信息。比較傳統(tǒng)的香農定理采樣與我們的采樣方法可知,我們方法的采樣頻率約是實際采樣頻率的0.024%(1/T/fs),采樣數(shù)據(jù)點數(shù)目約是傳統(tǒng)采樣方法的0.8%(41/4 160)。
文中提出了一種新型的多通道FRI信號的采集方法,詳細地介紹了它的信號采樣過程以及信號重建過程,豐富了次奈奎斯特采樣方式;為次奈奎斯特采樣在超聲波成像系統(tǒng)中的應用邁出了堅實的一步,并為以后的研究工作指明方向。
[1] Candes E,Wakin M.An introduction to compressive sampling[J].IEEE Signal Process Magazine,2008,25(2):21-30.
[2] Vetterli M,Marziliano P,Blu T.Sampling signals with?nite rate of innovation[J].IEEE Transactions on Signal Processing,2002,50(6):1417-1428.
[3] Dragotti P L,Vetterli M,Blu T.Sampling moments and reconstructing signals of fnite rate of innovation:Shannon meets strang-fix[J].IEEE Transactions on Signal Processing,2007,55(5):1741-1757.
[4] Kusuma J,Goyal V.Multichannel sampling of parametric signals with a successive approximation property[J].IEEE Int.Conf.Image Process,2006,50(6):1265-1268.
[5] Olkkonen H,Olkkonen J.Measurement and reconstruction of impulse train by parallel exponential filters[J].IEEE Signal Process.Lett.,2008,15(1):241-244.
[6] Stoica P,Moses R.Introduction to Spectral Analysis[M].3rd ed.Englewood Cliffs,NJ:Prentice-Hall,1997.