三維線性相位FIR數(shù)字濾波器對(duì)稱特性的設(shè)計(jì)

文/殷仕淑 武岳 常郝

1 引言

對(duì)多維信號(hào)處理的研究工作最早可以追溯到上世紀(jì)60年代。D.P. Petersen 與D. Middleton在文獻(xiàn)[1]中首次開展了對(duì)多維頻帶有限信號(hào)進(jìn)行抽樣和重建工作的理論研究。此后，在多維數(shù)字濾波器的研究領(lǐng)域里，致力于二維數(shù)字濾波器理論及設(shè)計(jì)方法的研究工作占據(jù)絕大部分分量，發(fā)展很快。相較而言，由于數(shù)據(jù)量大，系統(tǒng)復(fù)雜度增加，處理三維空間信號(hào)（例如三維醫(yī)學(xué)影像信號(hào)，地球物理空間信號(hào)等）的三維數(shù)字濾波器（Three dimensional digital filters: 3D DF）的相關(guān)理論及設(shè)計(jì)工作進(jìn)展不大。設(shè)計(jì)一個(gè)長度為10的FIR一維數(shù)字濾波器，最多有10個(gè)未知變量待求。設(shè)計(jì)一個(gè)三維長度都為10的FIR三維數(shù)字濾波器，待求的未知變量數(shù)為1000。雖然濾波器的設(shè)計(jì)問題描述相同，但龐大的未知變量數(shù)對(duì)于算法的有效收斂是個(gè)極大考驗(yàn)。雖然可以利用低維濾波器去處理三維信號(hào)，但考慮到數(shù)據(jù)之間的相關(guān)性，尤其是在移動(dòng)圖像處理、機(jī)器視覺、醫(yī)學(xué)圖像處理、地震信號(hào)處理和其他地球物理等三維空間信號(hào)處理領(lǐng)域中，發(fā)展三維數(shù)字濾波器的概念和研究三維數(shù)字濾波器的設(shè)計(jì)問題仍舊是不可避免的。

本文主要討論線性相位三維FIR數(shù)字濾波器的分類和設(shè)計(jì)工作。與IIR濾波器相比，F(xiàn)IR系統(tǒng)的單位脈沖長度有限，滿足絕對(duì)可和系統(tǒng)穩(wěn)定性判斷條件。其次，F(xiàn)IR數(shù)字濾波器可以實(shí)現(xiàn)線性相位響應(yīng)。最初FIR三維數(shù)字濾波器的設(shè)計(jì)是以一維數(shù)字濾波器為原型，采用McClellan變換(McClellan transform)進(jìn)行轉(zhuǎn)換得到的[2-7]。遺憾的是，這種方法無法精確控制三維數(shù)字濾波器的頻率特性。優(yōu)化迭代算法是另一種比較常見的用來設(shè)計(jì)FIR 三維數(shù)字濾波器方法。濾波器的設(shè)計(jì)問題寫成線性規(guī)劃問題求解，且不用約束其系統(tǒng)穩(wěn)定性[8,9]。在文獻(xiàn)[10,11]中，解析最小二乘法用來求解FIR三維數(shù)字濾波器。為了降低設(shè)計(jì)工作的計(jì)算復(fù)雜度，F(xiàn)IR三維數(shù)字濾波器的設(shè)計(jì)工作充分利用三維空間的對(duì)稱特性[4,5,8-13]。即便如此，在上述文獻(xiàn)中所獲得的FIR三維數(shù)字濾波器頻譜特性仍舊不太令人滿意。此外還有一些工作嘗試采用比較新穎的算法思想來解決FIR多維數(shù)字濾波器的設(shè)計(jì)問題，如模擬退火算法[14]，凸集投影法[15]，但遺憾的是三維數(shù)字濾波器的實(shí)驗(yàn)數(shù)據(jù)并沒有被成功獲取。

研究三維數(shù)字濾波器的對(duì)稱特性可以從頻域（幅度響應(yīng)，相位響應(yīng)）和空域（或稱時(shí)域，即系統(tǒng)單位脈沖響應(yīng)或?yàn)V波器系數(shù)）兩個(gè)方面入手，如：球面對(duì)稱和立方對(duì)稱等屬于頻域概念[16-18]。基于平面對(duì)稱概念所定義的多種不同幅度響應(yīng)對(duì)稱特性以列表的形式報(bào)告于[19]。在文獻(xiàn)[8]中，作者對(duì)濾波器的幅度平方響應(yīng)對(duì)稱性進(jìn)行研究，得出三維數(shù)字濾波器的幅度平方響應(yīng)共有47種不同類型對(duì)稱性。而在文獻(xiàn)[11]中，基于單位脈沖響應(yīng)三維數(shù)據(jù)的長度奇偶性和對(duì)稱奇偶性的不同組合，得出共有64種不同類型的線性相位FIR三維數(shù)字濾波器類型。類似地，在1998年，線性相位一維數(shù)字濾波器單位脈沖響應(yīng)對(duì)稱性的四種組合被移植至多維系統(tǒng)[20]。本文中，在文獻(xiàn)[11,20]工作的基礎(chǔ)上，擬對(duì)單位脈沖響應(yīng)三維數(shù)據(jù)對(duì)稱性的64種組合進(jìn)行分析，剔除冗余數(shù)，初步提出20種不同類型的線性相位FIR三維數(shù)字濾波器。結(jié)合零點(diǎn)位置分析，給出設(shè)計(jì)不同種類（低通、高通、帶通和帶阻）線性相位FIR三維數(shù)字濾波器的單位脈沖響應(yīng)對(duì)稱特性的選取。為了得到線性相位三維FIR數(shù)字濾波器，本文把此類型濾波器的設(shè)計(jì)問題寫出凸問題，采用約束凸規(guī)劃進(jìn)行求解[21]，獲得具有良好頻率特性的線性相位三維FIR數(shù)字濾波器。最后給出設(shè)計(jì)實(shí)例。

2 3D FIR線性相位數(shù)字濾波器

對(duì)于一個(gè)因果的FIR三維數(shù)字濾波器，其單位脈沖響應(yīng)可以用三維矩陣數(shù)據(jù)來表示：其系統(tǒng)函數(shù)表示為：

系統(tǒng)的頻率響應(yīng)定義為濾波器單位脈沖響應(yīng)的離散時(shí)間傅里葉變換（Discrete-time Fourier transform: DTFT），是式（1）中z變量在z平面單位圓上取值，即，

類似于一維數(shù)字濾波器單位脈沖響應(yīng)的對(duì)稱性分析， FIR三維數(shù)字濾波器在三維方向?yàn)V波器長度在每一維度方向上可能是奇數(shù)長度也可能是偶數(shù)長度，且在每一維度方向可能是奇對(duì)稱，也可能是偶對(duì)稱的。因此，F(xiàn)IR三維數(shù)字濾波器的單位脈沖響應(yīng)每一維度方向有四種可能組合：

表1：?jiǎn)挝幻}沖響應(yīng)三維方向?qū)ΨQ特性組合種類

（1）長度為奇數(shù)且偶對(duì)稱；

（2）長度為偶數(shù)且偶對(duì)稱；

（3）長度為奇數(shù)且奇對(duì)稱；

（4）長度為偶數(shù)且奇對(duì)稱。

也就是文獻(xiàn)[11]中所說的一共有64種組合。我們可以用數(shù)組（111），（112），（121），。。。，（444）來表示這64種組合。

利用線性相位濾波器單位脈沖響應(yīng)的數(shù)據(jù)對(duì)稱性，我們可以很方便地寫出表1中三維FIR線性相位濾波器的頻率響應(yīng)。例如，現(xiàn)有某FIR三維數(shù)字濾波器的單位脈沖響應(yīng)三維方向?qū)ΨQ特性為表1中組合“111”，則有：

其中，

表2：例1中三維線性相位FIR濾波器系數(shù)

發(fā)現(xiàn)：

圖1：三維FIR數(shù)字濾波器幅頻特性

圖2

再例如：對(duì)于組合（444），利用濾波器系數(shù)對(duì)稱特性，濾波器的頻率響應(yīng)可以寫成：

其中，

對(duì)于一維線性相位FIR濾波器的四種類型，類型1可用于設(shè)計(jì)低通（LP），帶通（BP），高通（HP），和帶阻（BS）濾波器；類型2可用于設(shè)計(jì)LP和BP濾波器；類型3僅能用于設(shè)計(jì)BP濾波器；類型4可用于設(shè)計(jì)HP和BP濾波器。因此，表1中20種類型中適用于設(shè)計(jì)在三個(gè)維度里具有相同通阻特性的線性相位FIR濾波器為：低通（111，112，122，222）、帶通（全部組合）、高通（111，114，144，444）和帶阻（111）。如果濾波器的三個(gè)維度具有不同的通阻特性，則根據(jù)具體要求可選擇不同組合種類。

在接下來的小節(jié)，我們把三維線性相位FIR濾波器組的設(shè)計(jì)問題寫出極小化極大值問題，帶入約束凸規(guī)劃求解。

3 線性相位3D FIR數(shù)字濾波器的設(shè)計(jì)

基于CVX工具箱的約束凸規(guī)劃法（DCP），用來求解各類凸問題[21]。與半定規(guī)劃（SDP）法和二階錐規(guī)劃（SOCP）法相比，DCP對(duì)問題描述簡(jiǎn)單直觀，求解復(fù)雜問題的能力更強(qiáng)。為了采用CVX求解三維線性相位FIR數(shù)字濾波器的設(shè)計(jì)問題，我們先把FIR三維數(shù)字濾波器的設(shè)計(jì)問題表達(dá)成極小化極大值問題：

采用DCP設(shè)計(jì)FIR一維或二維數(shù)字濾波器時(shí)，問題描述相同。唯一不同的是在（7）式中待求未知變量數(shù)是三維變量。在濾波器每維度系數(shù)長度一樣時(shí)FIR三維數(shù)字濾波器的待求變量是一維數(shù)字濾波器的倍。式（7）可以用來求任意相位響應(yīng)的FIR濾波器。對(duì)于線性相位FIR濾波器，可以利用其濾波器系數(shù)的對(duì)稱性，來降低待求變量數(shù)。只是對(duì)于表1中的不同組合形式，式5要做不同的改寫。例如組合111，式（7）改寫如下：

4 三維線性相位FIR數(shù)字濾波器的設(shè)計(jì)實(shí)例

設(shè)計(jì)三維線性相位FIR數(shù)字濾波器時(shí)，必須預(yù)先了解以下幾點(diǎn)：

（1）所期望的濾波器單位脈沖響應(yīng)長度；

（2）所期望的濾波器分段頻率特性（低通、高通等）及截止頻率點(diǎn)；

（3）所期望的濾波器通帶/阻帶形狀（球形、立方形等）；

（4）所期望的濾波器通帶或阻帶所能接受的最大波紋誤差。

一般來說，濾波器系數(shù)的長度、過渡帶寬度（通帶與阻帶之間的間隔）和通帶、阻帶上的波紋誤差之間是一種相互制衡的關(guān)系。濾波器系數(shù)越長，過渡帶越寬，則通帶和阻帶波紋誤差則越小。而通帶和阻帶波紋誤差之間也有相互制衡的關(guān)系，壓縮通帶的波紋誤差，也就是要獲得更好的通帶平坦性，則在其他參數(shù)不變的前提下，必然以增大阻帶波紋誤差作為代價(jià)。

例1：設(shè)計(jì)一個(gè)長度(5，5，5)，球面對(duì)稱，線性相位的低通FIR三維數(shù)字濾波器。通帶截止頻率和阻帶截止頻率分別為（0.3π，0.6π），濾波器的阻帶最大波紋誤差設(shè)為-30dB。把濾波器設(shè)計(jì)問題寫成式（8）進(jìn)行求解。待求系數(shù)c(n1,n2,n3)大小為（3×3×3）個(gè)。圖1中顯示其幅度頻率響應(yīng)。

由于所得濾波器在三個(gè)維度方向上具有對(duì)稱性，從另外兩個(gè)維度看進(jìn)去的濾波器幅頻特性與圖1中是一致的。其濾波器系數(shù)列于表2中，以供讀者參考。仔細(xì)觀察表2中數(shù)據(jù)，可發(fā)現(xiàn)每頁數(shù)據(jù)具有關(guān)于對(duì)角線對(duì)稱的特性；并且第一頁第二列數(shù)據(jù)與第二頁第一列數(shù)據(jù)相同，以此類推。表2中濾波器系數(shù)的對(duì)稱性緣于本例中所要求設(shè)計(jì)的濾波器在三個(gè)維度上通阻帶截止頻率點(diǎn)相同，濾波器長度一致，即濾波器在幅頻上具有對(duì)稱性。這種幅頻對(duì)稱特性也可以用來降低濾波器設(shè)計(jì)過程的未知變量數(shù)。但這種幅頻對(duì)稱特性在單位脈沖響應(yīng)上的影響較為復(fù)雜，尤其是濾波器比較長的時(shí)候。例2：設(shè)計(jì)一個(gè)長度(12，12，12)，通帶區(qū)域方形對(duì)稱，線性相位的低通FIR三維數(shù)字濾波器。通帶截止頻率和阻帶截止頻率分別為（0.3π，0.6π），濾波器的阻帶最大波紋誤差設(shè)為-30dB。待求系數(shù)c(n1,n2,n3)大小為63個(gè)。圖2顯示的是所得三維FIR數(shù)字濾波器幅頻響應(yīng)圖，其在通帶的最大波紋誤差為2.048e-3。由于篇幅所限，濾波器系數(shù)沒有列出來。

5 結(jié)束語

在處理三維數(shù)字信號(hào)時(shí)，考慮到數(shù)據(jù)之間的相關(guān)性，三維數(shù)字濾波器的設(shè)計(jì)不可避免。三維數(shù)字濾波器設(shè)計(jì)工作的難點(diǎn)在于龐大的待求未知變量數(shù)。本文從一維線性相位FIR濾波器系數(shù)的對(duì)稱性出發(fā)，首次提出三維線性相位FIR濾波器共有20種不同類型。利用濾波器系數(shù)的對(duì)稱性，大幅度減少濾波器設(shè)計(jì)過程中待求未知變量淑，把三維線性相位FIR數(shù)字濾波器的設(shè)計(jì)問題寫成極大值極小化問題，帶入凸規(guī)劃求解。實(shí)驗(yàn)結(jié)果顯示，這種方法可以獲得具有良好頻率特性的三維線性相位FIR濾波器。