核磁共振波譜儀數字濾波器的設計
摘要 介紹了最佳一致逼近理論的基本概念及Remez算法，描述了適用于棱磁共振波譜儀 的FIR—DF設計方法及實(shí)踐經(jīng)驗一提供丁實(shí)驗結果
I 近年來(lái) 數字處理技術(shù)逐漸被引人核磁共振(NMR)渡譜實(shí)驗技術(shù)中，取得了可喜的成果． 圖1是我們在NMR波譜儀中增加數字濾波器(DF)之后所傲試驗得到的實(shí)際譜圖，可見(jiàn)DF 在完成指定的功能(低通，帶通等)時(shí)．其結果是令人滿(mǎn)意的，即過(guò)渡帶十分陡峭．實(shí)驗與分析表 明、在NMR實(shí)驗中增加DF具有如下優(yōu)勢：1)消除了譜圖的折疊現象，從而抑制了偽峰的產(chǎn) 生．露點(diǎn)儀| 亮度計| 溫度記錄儀| 溫濕度記錄儀| 光功率計| 粒子計數器| 粉塵計|這是模擬濾波器所無(wú)法辦到的；2)在整個(gè)頻率范圍內都保持線(xiàn)性相移．從而消減了基線(xiàn)波 動(dòng)；3)平坦的幅頻特性，既避免了譜圖邊緣峰的遺褥一又為精確的譜圖積分提供了先決條件；4) 結合高速(例如，16倍奈奎斯特頻率)采樣，數據抽取技術(shù)，即能以高信噪比處理小譜峰．綜上 所述．以數字濾波技術(shù)取代傳統的模擬濾波技術(shù)，將使NMR實(shí)驗技術(shù)的發(fā)展躍上一個(gè)新的臺 階．本文擬論述NMR DF的設計原理與應用．

數字濾波器可分為無(wú)限沖激響應(IIR)和有限沖激響應(FIR)兩類(lèi)．FIR DF具有如下特 點(diǎn)：1)滿(mǎn)足指定的設計指標所要的階數較低；2)具有嚴格的線(xiàn)性相位特性；3)宜采用非遞歸結 構，因而有限精度的運算穩定且誤差較�。�4)可直接采用快速付立葉變換(FFT)；5)設計靈活． 適應面廣(除可做多通帶濾波器之外，還可做差分器和希爾伯特變換器)．因此，盡管FIR—DF 的設計無(wú)法借助現有的模擬濾波器設計公式，筆者仍主張采用FIR—DF，因為其不足之處可以 由計算機輔助設計技術(shù)加以彌補． 眾所周知，一個(gè)理想的濾波器應當具有的特性是： c 一 窩 ( )： K ·60 ’ (1b) 60 一 (Ic) 即要求：1)幅頻特性，在通帶內的相對衰減為1，在阻帶內的相對衰為0，且平坦．2)相頻特性為 線(xiàn)性．3)通帶邊緣 與阻帶邊緣 相重合．例如一個(gè)理想的低通濾波器應具有圖2實(shí)線(xiàn)所 示的特性．如上所述，FIR DF固有線(xiàn)性相位特性，但幅頻特性則與其它濾波器一樣·表現為一 定的波動(dòng)(通帶在士 內波動(dòng)，阻帶在士 內波動(dòng))．另外，不可能實(shí)現 = ．即存在 一定的過(guò)渡帶 一I 一 I，如圖2中虛 線(xiàn)所示． 因此，FIR—DF的設計目標就是獲得一 個(gè)參數集，以保證波動(dòng)值處于指定的范圍 內(即波動(dòng)≤ ± )，而過(guò)渡帶滿(mǎn)足指定的 和 值．

一個(gè)能實(shí)現這一目標的設計方法 是“一致逼近法”(也稱(chēng)等紋波逼近法)．其 涵義不難從圖2得到說(shuō)明，即在指定的通 帶(或／和阻帶)內其幅度的波動(dòng)可能形成 多個(gè)極值點(diǎn)(極大或極小)．然而波動(dòng)的幅 度 都是相等的(注意．一一般地 ≠ ：)，習 慣上稱(chēng)波動(dòng)值為偏差． 一致逼近法的理論基礎是“切比雪夫 最佳一致逼近定理”．這個(gè)定理的基本意思 是：對于給定的閉區間 ．6]上的連續函數 ，( )．在所有m 階多項式的集合 ．尋找 一個(gè)多項式 ( )使它在 ，陽(yáng)上對于，(z) 的最大偏差和其它屬于 的多項式戶(hù)(z) 對于，( )的最大偏差相比是最小的，即 I (z)一，(z)I — rain{maxI ( )一f(z)j) (2) 口≤ ^ 該理論指出，這樣一個(gè)多項式戶(hù)(z)是存在 的，且是唯一的．構造戶(hù)(z)的方法就是有名 的“交錯點(diǎn)組定理”．該定理是：設，(z)是 定義在_d．陽(yáng)上的連續函數， (z)為多項式 集合 中一個(gè)階次不超過(guò) 的多項式， 并令 一max I (z)-f(x)I 1 ) (3a) (z)一Ip(x)一f )f J a)沒(méi)有施加濾波b)帶通濾渡c)低通濾波 Fig．1 NM R spectrum 則戶(hù)( )是，(z)的最佳一致逼近的充要件是戶(hù)(z)在[ ．6]上至少存在(m-2)個(gè)交錯點(diǎn)構成的 集合{z }．使得 一± ’ — �！�⋯ +。 l (3b) 3 (x )一一 (薯+1)， 一】，2，⋯ ， +2』 這個(gè)( +2)個(gè)點(diǎn)即是交錯點(diǎn)組． ( )為偏差函數．式(3)的涵義是明確的．即： ． 。．·-- +。．都 是偏差函數 ( )的極值點(diǎn)(極大或極小)，它們對， )都具有最大的偏差，且其偏差值都等于 以，其它各處( 在{z ))的偏差都小于靠；而 )的各項系數可以由交錯點(diǎn)組對應的極值組確 定．

數字信號處理理論表明：1)濾波器的系統特性． 由其單位沖激響應h(n)充分體現；2)h(n)在z平面 單位圓上的z變換H( )就是濾波器的頻率響應特 性．顯然，濾波器的設計是在頻率域上進(jìn)行的，關(guān)鍵 在于根據設定的指標，在頻率域上求解 ( )．而 后，對H(d )進(jìn)行反z變換，求得序列h( )． 對于長(cháng)度為N 的序列^ )，其z變換為 H( )一Σ^(n)P一一 (4) 圖2濾波器特性 簡(jiǎn)寫(xiě)為：以H( )代替通常的H( )形式．對于NMR Fig·2 Filter specifie[ty 實(shí)驗，FIR～DF宜選�、� 為奇數；且選取h( )為偶對稱(chēng)+即h(一)一h(Ⅳ 一】 )，稍加演算，即 可將式(4)寫(xiě)成 ( )-C- Σ 口 )cos( ) (5) n一0 f^( )， 一。 舯 ⋯ ⋯ ⋯ ，丁N 1 (6a) 一 亦即^( ) 1 ( ～ ) 卜 _ _～ ，H”—一oO， 1， - t／_ 1)一 ， 式(5)中的相移 ( )：一( ) 一一K 就是線(xiàn)性相位特性；其幅頻特性則為 H( )：Σn )cos(noJ) (7) 通過(guò)式(3)，可以把式(1)和式(7)聯(lián)系起來(lái)：1)以H ( )代替， )；2)以H(∞)代替； )； 3)以 代替z，在數字信號處理領(lǐng)域中， 稱(chēng)為數字角頻率，它的取值范圍由0到 ( 到2 與 其對稱(chēng))，且構成閉區間[0， ]．于是有 一n1ax IH ( )一H ( )ll j (8) ( )一H ( ) (∞) ．I 式中下標F-[o， ]．針對圖2的具體情況，且注意到通常n≠ ：，命 ( )：j{，通帶內 一蓑 (9) 1，阻帶內 則式(8)可寫(xiě)成 8(oJ)一Ⅳ )IH ( )一H (∞)1 (10) ㈨ 聊 對 ” 一 將交錯點(diǎn)組定理應用到式(10)，就是說(shuō)：H( )在區間 ’上對H ( )唯一最佳一致逼近舶 充要條件是偏差函數 (m)在 、上至少呈現(M+2)個(gè)交錯點(diǎn)．構成頻率點(diǎn)的集合{ ’．使得 l ( )l— r + )r一 如一max J f )J F 如果已知在F上的集合 }中這( 一2)個(gè)頻率點(diǎn)，則由式(1 0)有 W ( )E-h’ ( )一H ( ) =(一】) 代人 ( )則是(注意，求和上限的M一竺 ) ( ) ( )～ 2J口(n)cosO )]= 【一1) 以 l1 式(11)是一個(gè)r 1。2)×f +2)的矩陣 l cos(~、) cos(~：) _ ● ● cos(~M) c0s(w~十 ， ) cos(2w ) - c∞(2 ee ) ● ● _ c0s(2 M) cos(2 M+】) COS(Mo,。) cos(M eo J) cos(Mw!) cos(M Ⅵ1 cos(朋 1) 1／W(o~ ) 1／w( ) 1／w ( 1) ／W (％D ) ( +J) 解此方程組，可以唯一地求得系數n(o)，a(1)，⋯ ， ( )及 ．因此得到最佳濾波器的頻率響應 (幅頻特性)H(m)．然而這僅是數學(xué)上的認識，當M 稍大一點(diǎn)．

解這個(gè)方程談何容易．為此Mc． Clallan等人EzJ利用數值分析中的Remez算法，靠一次一次迭代求得一組交錯點(diǎn)組，而且每一 次迭代過(guò)程中避免求解式(12)．Remez算法可歸納如下： 1)在頻率集F上等間隔地取( +2)個(gè)頻率點(diǎn)foot,． ， ． 測位置，然后計算 + L — 1 一0n譯 t cOS％— C而OSm~ j + Σ a(k)tt ( ) = — — Σ(-1) ( )／ ( ) 叫 一 做為交錯點(diǎn)組的初始猜 (1 3a) (10b) ^ 0 即是此次猜測的交錯點(diǎn)組集合{ }形成的偏差，它本質(zhì)上就是圖2中的島．得到8 后．根據 重心形式的拉格朗日插值公式求 一H ( )一(- 1) ． 一。，1，⋯ ，M (13c) 超’~ a(k )7 gz 0 s 一cos (1 3(1) 然后，依次檢查 ㈤ ㈤ ∞ ” " 甜 ～ { ( ))中的每一個(gè)元素是否滿(mǎn)足(Id( )t≤ t t)．若滿(mǎn)足，此次猜測韻頻率交錯點(diǎn)組{ }就 是最佳一致逼近交錯點(diǎn)組，計算可告結束．但是，一般需要進(jìn)行多次迭代，即轉入下．步． 圈3 Remez算法流程圈 Fig．3 Remez algorithm flowchart Ⅳ ， ，¨ ， ·(∞)， (m) 技木參嘲I人 工 在[O， 上分格點(diǎn) ~l× I6 I 試探 +】十撮值氨 算法 ． 1 得到 (∞) 計算^(-) 上 jI出 ，赴，^(^) 圈4 最佳一致逼近計算流程囝 Fig 4 O U．A calculate flowchart 3)對此次猜測的{ )中每一個(gè)元素q，在其 甜近依式(13d)計算 ( )，搜索是否有f ( )f> 的頻率點(diǎn) ，．若有則以 ，迭代 ．如此構成⋯ 個(gè)新的頻率點(diǎn)集合{ 。)取代{ ) 再返回計算 (】)、(2)以至(3)，直到滿(mǎn)足Ia( )j≤ j以I才結柬 計算．即可得到 =d 和H )的最佳一致逼近 系數,if／(0)，口(1)，d(2)．⋯ ，n( )．然后．根據式(6 b)計算出^(0)，^(1)．⋯ ⋯，h(Ⅳ一1)序列則完成 了反z變換(注：在循環(huán)迭代中d 將增加)． 圖3為Remez算法流程圖．運行Remez程序時(shí)，為保證 )是對 ( )的最佳一致逼 近，誤差函數8( )必須呈現肘+2個(gè)“交錯”，即d( )至少有( +2)個(gè)極值，且交替改變符號． 式(105中由于Ⅳ( )和H ( )都是常數，所以 ( )的極值就是H( )的極值．由于式(7)是～ 余弦函數．對其微分則為“一sin( ”．可見(jiàn)，在 一0和OA= 處必為極大值或極小值，因此．

在 區間 o， ]上，H( )至多有(肘+1j個(gè)極值．圖2中虛線(xiàn)相應于N一1 3( 一6)的情況． 圖4是計算最佳一致逼近的計算機程序流程圖．程序的輸人參數是：^(”)序列的長(cháng)度Ⅳ ； 通帶邊緣頻率 ；阻帶邊頻 ；以及H (。)和W( )特性(均由專(zhuān)門(mén)的子程序輸入)頻帶的(遙、 阻)數目可以選擇(一般限制在10以?xún)?．頻率軸。的均勻分隔點(diǎn)數也是可選擇的(可約定為 16)．由于是以歸一 化為1的頻率采樣，所以 的數值域 O． ]對應于[O．0．5]．N 的最大值～ 般限定在128以?xún)龋�程序的輸出是a-、赴 ^( )． 由于 是由程序運行之后得出的，所以當 它們不能達到指定的偏差值時(shí)．可增加 的 輸人值以提高度計算精度．由圖2可以看出， 過(guò)渡帶的寬度4 大致為半個(gè)板值交錯點(diǎn)閫的 寬度，所以，當M (一 )較大時(shí)．邊緣是十 分陡峭的．

圖5是筆者利用一個(gè)8 Hz的點(diǎn)阻FIR DF濾去某樣品中的水溶劑峰之后所得到的實(shí) 際NMR譜圖．由圖可見(jiàn)．基本上消除了基線(xiàn) 漂移，使譜圖更為精細．筆者的經(jīng)驗表明，宜采 用高、低通組合的方式實(shí)現點(diǎn)阻濾波+否剮將 由于過(guò)渡帶的存在而難 實(shí)現窄帶的點(diǎn)阻濾 波器． 實(shí)驗中運行的程序是用C++語(yǔ)言編寫(xiě) 的．順便指出，文獻[2]給出的程序實(shí)際上還可 以運行差分器和希爾伯特變換器的設計程序．