核磁共振波譜儀字濾波器的設(shè)

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

數(shù)字濾波器可分為無限沖激響應(yīng)(IIR)和有限沖激響應(yīng)(FIR)兩類．FIR DF具有如下特 點(diǎn)：1)滿足指定的設(shè)計(jì)指標(biāo)所要的階數(shù)較低；2)具有嚴(yán)格的線性相位特性；3)宜采用非遞歸結(jié) 構(gòu)，因而有限精度的運(yùn)算穩(wěn)定且誤差較�。�4)可直接采用快速付立葉變換(FFT)；5)設(shè)計(jì)靈活． 適應(yīng)面廣(除可做多通帶濾波器之外，還可做差分器和希爾伯特變換器)．因此，盡管FIR—DF 的設(shè)計(jì)無法借助現(xiàn)有的模擬濾波器設(shè)計(jì)公式，筆者仍主張采用FIR—DF，因?yàn)槠洳蛔阒�處可�?由計(jì)算機(jī)輔助設(shè)計(jì)技術(shù)加以彌補(bǔ)．
眾所周知，一個(gè)理想的濾波器應(yīng)當(dāng)具有的特性是： 即要求：1)幅頻特性，在通帶內(nèi)的相對(duì)衰減為1，在阻帶內(nèi)的相對(duì)衰為0，且平坦．2)相頻特性為 線性．3)通帶邊緣 與阻帶邊緣 相重合．例如一個(gè)理想的低通濾波器應(yīng)具有圖2實(shí)線所 示的特性．如上所述，F(xiàn)IR DF固有線性相位特性，但幅頻特性則與其它濾波器一樣·表現(xiàn)為一 定的波動(dòng)(通帶在士 內(nèi)波動(dòng)，阻帶在士 內(nèi)波動(dòng))．另外，不可能實(shí)現(xiàn) = ．即存在 一定的過渡帶 一I 一 I，如圖2中虛 線所示． 因此，F(xiàn)IR—DF的設(shè)計(jì)目標(biāo)就是獲得一 個(gè)參數(shù)集，以保證波動(dòng)值處于指定的范圍 內(nèi)(即波動(dòng)≤ ± )，而過渡帶滿足指定的 和 值．
一個(gè)能實(shí)現(xiàn)這一目標(biāo)的設(shè)計(jì)方法 是“一致逼近法”(也稱等紋波逼近法)．其 涵義不難從圖2得到說明，即在指定的通 帶(或／和阻帶)內(nèi)其幅度的波動(dòng)可能形成 多個(gè)極值點(diǎn)(極大或極小)．然而波動(dòng)的幅 度 都是相等的(注意．一一般地 ≠ ：)，習(xí) 慣上稱波動(dòng)值為偏差． 一致逼近法的理論基礎(chǔ)是“切比雪夫 最佳一致逼近定理”．這個(gè)定理的基本意思 是：對(duì)于給定的閉區(qū)間 ．6]上的連續(xù)函數(shù) ，( )．在所有m 階多項(xiàng)式的集合 ．尋找 一個(gè)多項(xiàng)式 ( )使它在 ，陽上對(duì)于，(z) 的最大偏差和其它屬于 的多項(xiàng)式戶(z) 對(duì)于，( )的最大偏差相比是最小的，即 I (z)一，(z)I — rain{maxI ( )一f(z)j) (2) 口≤ ^ 該理論指出，這樣一個(gè)多項(xiàng)式戶(z)是存在 的，且是唯一的．構(gòu)造戶(z)的方法就是有名 的“交錯(cuò)點(diǎn)組定理”．該定理是：設(shè)，(z)是 定義在_d．陽上的連續(xù)函數(shù)， (z)為多項(xiàng)式 集合 中一個(gè)階次不超過 的多項(xiàng)式， 并令 一max I (z)-f(x)I 1 ) (3a) (z)一Ip(x)一f )f J a)沒有施加濾波b)帶通濾渡c)低通濾波 Fig．1 NM R spectrum 則戶( )是，(z)的最佳一致逼近的充要件是戶(z)在[ ．6]上至少存在(m-2)個(gè)交錯(cuò)點(diǎn)構(gòu)成的 集合{z }．使得 一± ’ — �！�⋯ +。 l (3b) 3 (x )一一 (薯+1)， 一】，2，⋯ ， +2』 這個(gè)( +2)個(gè)點(diǎn)即是交錯(cuò)點(diǎn)組． ( )為偏差函數(shù)．式(3)的涵義是明確的．即： ． 。．·-- +。．都 是偏差函數(shù) ( )的極值點(diǎn)(極大或極小)，它們對(duì)， )都具有最大的偏差，且其偏差值都等于 以，其它各處( 在{z ))的偏差都小于靠；而 )的各項(xiàng)系數(shù)可以由交錯(cuò)點(diǎn)組對(duì)應(yīng)的極值組確 定．

數(shù)字信號(hào)處理理論表明：1)濾波器的系統(tǒng)特性． 由其單位沖激響應(yīng)h(n)充分體現(xiàn)；2)h(n)在z平面 單位圓上的z變換H( )就是濾波器的頻率響應(yīng)特 性．顯然，濾波器的設(shè)計(jì)是在頻率域上進(jìn)行的，關(guān)鍵 在于根據(jù)設(shè)定的指標(biāo)，在頻率域上求解 ( )．而 后，對(duì)H(d )進(jìn)行反z變換，求得序列h( )． 對(duì)于長度為N 的序列^ )，其z變換為 H( )一Σ^(n)P一一 (4) 圖2濾波器特性 簡寫為：以H( )代替通常的H( )形式．對(duì)于NMR Fig·2 Filter specifie[ty 實(shí)驗(yàn)，F(xiàn)IR～DF宜選�、� 為奇數(shù)；且選取h( )為偶對(duì)稱+即h(一)一h(Ⅳ 一】 )，稍加演算，即 可將式(4)寫成 ( )-C- Σ 口 )cos( ) (5) n一0 f^( )， 一。 舯 ⋯ ⋯ ⋯ ，丁N 1 (6a) 一 亦即^( ) 1 ( ～ ) 卜 _ _～ ，H”—一oO， 1， - t／_ 1)一 ， 式(5)中的相移 ( )：一( ) 一一K 就是線性相位特性；其幅頻特性則為 H( )：Σn )cos(noJ) (7) 通過式(3)，可以把式(1)和式(7)聯(lián)系起來：1)以H ( )代替， )；2)以H(∞)代替； )； 3)以 代替z，在數(shù)字信號(hào)處理領(lǐng)域中， 稱為數(shù)字角頻率，它的取值范圍由0到 ( 到2 與 其對(duì)稱)，且構(gòu)成閉區(qū)間[0， ]．于是有 一n1ax IH ( )一H ( )ll j (8) ( )一H ( ) (∞) ．I 式中下標(biāo)F-[o， ]．針對(duì)圖2的具體情況，且注意到通常n≠ ：，命 ( )：j{，通帶內(nèi) 一蓑 (9) 1，阻帶內(nèi) 則式(8)可寫成 8(oJ)一Ⅳ )IH ( )一H (∞)1 (10) ㈨ 聊 對(duì) ” 一 將交錯(cuò)點(diǎn)組定理應(yīng)用到式(10)，就是說：H( )在區(qū)間 ’上對(duì)H ( )唯一最佳一致逼近舶 充要條件是偏差函數(shù) (m)在 、上至少呈現(xiàn)(M+2)個(gè)交錯(cuò)點(diǎn)．
構(gòu)成頻率點(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) 解此方程組，可以唯一地求得系數(shù)n(o)，a(1)，⋯ ， ( )及 ．因此得到最佳濾波器的頻率響應(yīng) (幅頻特性)H(m)．然而這僅是數(shù)學(xué)上的認(rèn)識(shí)，當(dāng)M 稍大一點(diǎn)．解這個(gè)方程談何容易．為此Mc． Clallan等人EzJ利用數(shù)值分析中的Remez算法，靠一次一次迭代求得一組交錯(cuò)點(diǎn)組，而且每一 次迭代過程中避免求解式(12)．Remez算法可歸納如下： 1)在頻率集F上等間隔地取( +2)個(gè)頻率點(diǎn)foot,． ， ． 測(cè)位置，然后計(jì)算 + L — 1 一0n譯 t cOS％— C而OSm~ j + Σ a(k)tt ( ) = — — Σ(-1) ( )／ ( ) 叫 一 做為交錯(cuò)點(diǎn)組的初始猜 (1 3a) (10b) ^ 0 即是此次猜測(cè)的交錯(cuò)點(diǎn)組集合{ }形成的偏差，它本質(zhì)上就是圖2中的島．得到8 后．根據(jù) 重心形式的拉格朗日插值公式求 一H ( )一(- 1) ． 一。，1，⋯ ，M (13c) 超’~ a(k )7 gz 0 s 一cos (1 3(1) 然后，依次檢查 ㈤ ㈤ ∞ ” " 甜 ～ { ( ))中的每一個(gè)元素是否滿足(Id( )t≤ t t)．若滿足，此次猜測(cè)韻頻率交錯(cuò)點(diǎn)組{ }就 是最佳一致逼近交錯(cuò)點(diǎn)組，計(jì)算可告結(jié)束．但是，一般需要進(jìn)行多次迭代，即轉(zhuǎn)入下．步． 圈3 Remez算法流程圈 Fig．3 Remez algorithm flowchart Ⅳ ， ，¨ ， ·(∞)， (m) 技木參嘲I人 工 在[O， 上分格點(diǎn) ~l× I6 I 試探 +】十撮值氨 算法 ． 1 得到 (∞) 計(jì)算^(-) 上 jI出 ，赴，^(^) 圈4 最佳一致逼近計(jì)算流程囝 Fig 4 O U．A calculate flowchart 3)對(duì)此次猜測(cè)的{ )中每一個(gè)元素q，在其 甜近依式(13d)計(jì)算 ( )，搜索是否有f ( )f> 的頻率點(diǎn) ，．若有則以 ，迭代 ．如此構(gòu)成⋯ 個(gè)新的頻率點(diǎn)集合{ 。)取代{ ) 再返回計(jì)算 (】)、(2)以至(3)，直到滿足Ia( )j≤ j以I才結(jié)柬 計(jì)算．即可得到 =d 和H )的最佳一致逼近 系數(shù),if／(0)，口(1)，d(2)．⋯ ，n( )．然后．根據(jù)式(6 b)計(jì)算出^(0)，^(1)．⋯ ⋯，h(Ⅳ一1)序列則完成 了反z變換(注：在循環(huán)迭代中d 將增加)． 圖3為Remez算法流程圖．運(yùn)行Remez程序時(shí)，為保證 )是對(duì) ( )的最佳一致逼 近，誤差函數(shù)8( )必須呈現(xiàn)肘+2個(gè)“交錯(cuò)”，即d( )至少有( +2)個(gè)極值，且交替改變符號(hào)． 式(105中由于Ⅳ( )和H ( )都是常數(shù)，所以 ( )的極值就是H( )的極值．由于式(7)是～ 余弦函數(shù)．對(duì)其微分則為“一sin( ”．可見，在 一0和OA= 處必為極大值或極小值，因此．在 區(qū)間 o， ]上，H( )至多有(肘+1j個(gè)極值．圖2中虛線相應(yīng)于N一1 3( 一6)的情況． 圖4是計(jì)算最佳一致逼近的計(jì)算機(jī)程序流程圖．程序的輸人參數(shù)是：^(”)序列的長度Ⅳ ； 通帶邊緣頻率 ；阻帶邊頻 ；以及H (。)和W( )特性(均由專門的子程序輸入)頻帶的(遙、 阻)數(shù)目可以選擇(一般限制在10以內(nèi))．頻率軸。的均勻分隔點(diǎn)數(shù)也是可選擇的(可約定為 16)．由于是以歸一 化為1的頻率采樣，所以 的數(shù)值域 O． ]對(duì)應(yīng)于[O．0．5]．N 的最大值～ 般限定在128以內(nèi)．程序的輸出是a-、赴 ^( )． 由于 是由程序運(yùn)行之后得出的，所以當(dāng) 它們不能達(dá)到指定的偏差值時(shí)．可增加 的 輸人值以提高度計(jì)算精度．由圖2可以看出， 過渡帶的寬度4 大致為半個(gè)板值交錯(cuò)點(diǎn)閫的 寬度，所以，當(dāng)M (一 )較大時(shí)．邊緣是十 分陡峭的．

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