核磁共振波谱仪字滤波器的设

核磁共振波谱仪字滤波器的设
摘要 介绍了最佳一致逼近理论的基本概念及Remez算法，描述了适用于棱磁共振波谱仪 的FIR—DF设计方法及实践经验一提供丁实验结果
近年来 数字处理技术逐渐被引人核磁共振(NMR)渡谱实验技术中，取得了可喜的成果． 图1是我们在NMR波谱仪中增加数字滤波器(DF)之后所傲试验得到的实际谱图，可见DF 在完成指定的功能(低通，带通等)时．其结果是令人满意的，即过渡带十分陡峭．实验与分析表 明、在NMR实验中增加DF具有如下优势：1)消除了谱图的折叠现象，从而抑制了伪峰的产 生．这是模拟滤波器所无法办到的；2)在整个频率范围内都保持线性相移．从而消减了基线波 动；3)平坦的幅频特性，既避免了谱图边缘峰的遗褥一又为精确的谱图积分提供了先决条件；4) 结合高速(例如，16倍奈奎斯特频率)采样，数据抽取技术，即能以高信噪比处理小谱峰．综上 所述．以数字滤波技术取代传统的模拟滤波技术，将使NMR实验技术的发展跃上一个新的台 阶．本文拟论述NMR DF的设计原理与应用．

数字滤波器可分为无限冲激响应(IIR)和有限冲激响应(FIR)两类．FIR DF具有如下特 点：1)满足指定的设计指标所要的阶数较低；2)具有严格的线性相位特性；3)宜采用非递归结 构，因而有限精度的运算稳定且误差较小：4)可直接采用快速付立叶变换(FFT)；5)设计灵活． 适应面广(除可做多通带滤波器之外，还可做差分器和希尔伯特变换器)．因此，尽管FIR—DF 的设计无法借助现有的模拟滤波器设计公式，笔者仍主张采用FIR—DF，因为其不足之处可以 由计算机辅助设计技术加以弥补．
众所周知，一个理想的滤波器应当具有的特性是： 即要求：1)幅频特性，在通带内的相对衰减为1，在阻带内的相对衰为0，且平坦．2)相频特性为 线性．3)通带边缘 与阻带边缘 相重合．例如一个理想的低通滤波器应具有图2实线所 示的特性．如上所述，FIR DF固有线性相位特性，但幅频特性则与其它滤波器一样·表现为一 定的波动(通带在士 内波动，阻带在士 内波动)．另外，不可能实现 = ．即存在 一定的过渡带 一I 一 I，如图2中虚 线所示． 因此，FIR—DF的设计目标就是获得一 个参数集，以保证波动值处于指定的范围 内(即波动≤ ± )，而过渡带满足指定的 和 值．
一个能实现这一目标的设计方法 是“一致逼近法”(也称等纹波逼近法)．其 涵义不难从图2得到说明，即在指定的通 带(或／和阻带)内其幅度的波动可能形成 多个极值点(极大或极小)．然而波动的幅 度 都是相等的(注意．一一般地 ≠ ：)，习 惯上称波动值为偏差． 一致逼近法的理论基础是“切比雪夫 最佳一致逼近定理”．这个定理的基本意思 是：对于给定的闭区间 ．6]上的连续函数 ，( )．在所有m 阶多项式的集合 ．寻找 一个多项式 ( )使它在 ，阳上对于，(z) 的最大偏差和其它属于 的多项式户(z) 对于，( )的最大偏差相比是最小的，即 I (z)一，(z)I — rain{maxI ( )一f(z)j) (2) 口≤ ^ 该理论指出，这样一个多项式户(z)是存在 的，且是唯一的．构造户(z)的方法就是有名 的“交错点组定理”．该定理是：设，(z)是 定义在_d．阳上的连续函数， (z)为多项式 集合 中一个阶次不超过 的多项式， 并令 一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)个交错点构成的 集合{z }．使得 一± ’ — 。’⋯ +。 l (3b) 3 (x )一一 (薯+1)， 一】，2，⋯ ， +2』 这个( +2)个点即是交错点组． ( )为偏差函数．式(3)的涵义是明确的．即： ． 。．·-- +。．都 是偏差函数 ( )的极值点(极大或极小)，它们对， )都具有最大的偏差，且其偏差值都等于 以，其它各处( 在{z ))的偏差都小于靠；而 )的各项系数可以由交错点组对应的极值组确 定．

数字信号处理理论表明：1)滤波器的系统特性． 由其单位冲激响应h(n)充分体现；2)h(n)在z平面 单位圆上的z变换H( )就是滤波器的频率响应特 性．显然，滤波器的设计是在频率域上进行的，关键 在于根据设定的指标，在频率域上求解 ( )．而 后，对H(d )进行反z变换，求得序列h( )． 对于长度为N 的序列^ )，其z变换为 H( )一Σ^(n)P一一 (4) 图2滤波器特性 简写为：以H( )代替通常的H( )形式．对于NMR Fig·2 Filter specifie[ty 实验，FIR～DF宜选取Ⅳ 为奇数；且选取h( )为偶对称+即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)联系起来：1)以H ( )代替， )；2)以H(∞)代替； )； 3)以 代替z，在数字信号处理领域中， 称为数字角频率，它的取值范围由0到 ( 到2 与 其对称)，且构成闭区间[0， ]．于是有 一n1ax IH ( )一H ( )ll j (8) ( )一H ( ) (∞) ．I 式中下标F-[o， ]．针对图2的具体情况，且注意到通常n≠ ：，命 ( )：j{，通带内 一蓑 (9) 1，阻带内 则式(8)可写成 8(oJ)一Ⅳ )IH ( )一H (∞)1 (10) ㈨ 聊 对 ” 一 将交错点组定理应用到式(10)，就是说：H( )在区间 ’上对H ( )唯一最佳一致逼近舶 充要条件是偏差函数 (m)在 、上至少呈现(M+2)个交错点．
构成频率点的集合{ ’．使得 l ( )l— r + )r一 如一max J f )J F 如果已知在F上的集合 }中这( 一2)个频率点，则由式(1 0)有 W ( )E-h’ ( )一H ( ) =(一】) 代人 ( )则是(注意，求和上限的M一竺 ) ( ) ( )～ 2J口(n)cosO )]= 【一1) 以 l1 式(11)是一个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)．然而这仅是数学上的认识，当M 稍大一点．解这个方程谈何容易．为此Mc． Clallan等人EzJ利用数值分析中的Remez算法，靠一次一次迭代求得一组交错点组，而且每一 次迭代过程中避免求解式(12)．Remez算法可归纳如下： 1)在频率集F上等间隔地取( +2)个频率点foot,． ， ． 测位置，然后计算 + L — 1 一0n译 t cOS％— C而OSm~ j + Σ a(k)tt ( ) = — — Σ(-1) ( )／ ( ) 叫 一 做为交错点组的初始猜 (1 3a) (10b) ^ 0 即是此次猜测的交错点组集合{ }形成的偏差，它本质上就是图2中的岛．得到8 后．根据 重心形式的拉格朗日插值公式求 一H ( )一(- 1) ． 一。，1，⋯ ，M (13c) 超’~ a(k )7 gz 0 s 一cos (1 3(1) 然后，依次检查 ㈤ ㈤ ∞ ” " 甜 ～ { ( ))中的每一个元素是否满足(Id( )t≤ t t)．若满足，此次猜测韵频率交错点组{ }就 是最佳一致逼近交错点组，计算可告结束．但是，一般需要进行多次迭代，即转入下．步． 圈3 Remez算法流程圈 Fig．3 Remez algorithm flowchart Ⅳ ， ，¨ ， ·(∞)， (m) 技木参嘲I人 工 在[O， 上分格点 ~l× I6 I 试探 +】十撮值氨 算法 ． 1 得到 (∞) 计算^(-) 上 jI出 ，赴，^(^) 圈4 最佳一致逼近计算流程囝 Fig 4 O U．A calculate flowchart 3)对此次猜测的{ )中每一个元素q，在其 甜近依式(13d)计算 ( )，搜索是否有f ( )f> 的频率点 ，．若有则以 ，迭代 ．如此构成⋯ 个新的频率点集合{ 。)取代{ ) 再返回计算 (】)、(2)以至(3)，直到满足Ia( )j≤ j以I才结柬 计算．即可得到 =d 和H )的最佳一致逼近 系数,if／(0)，口(1)，d(2)．⋯ ，n( )．然后．根据式(6 b)计算出^(0)，^(1)．⋯ ⋯，h(Ⅳ一1)序列则完成 了反z变换(注：在循环迭代中d 将增加)． 图3为Remez算法流程图．运行Remez程序时，为保证 )是对 ( )的最佳一致逼 近，误差函数8( )必须呈现肘+2个“交错”，即d( )至少有( +2)个极值，且交替改变符号． 式(105中由于Ⅳ( )和H ( )都是常数，所以 ( )的极值就是H( )的极值．由于式(7)是～ 余弦函数．对其微分则为“一sin( ”．可见，在 一0和OA= 处必为极大值或极小值，因此．在 区间 o， ]上，H( )至多有(肘+1j个极值．图2中虚线相应于N一1 3( 一6)的情况． 图4是计算最佳一致逼近的计算机程序流程图．程序的输人参数是：^(”)序列的长度Ⅳ ； 通带边缘频率 ；阻带边频 ；以及H (。)和W( )特性(均由专门的子程序输入)频带的(遥、 阻)数目可以选择(一般限制在10以内)．频率轴。的均匀分隔点数也是可选择的(可约定为 16)．由于是以归一 化为1的频率采样，所以 的数值域 O． ]对应于[O．0．5]．N 的最大值～ 般限定在128以内．程序的输出是a-、赴 ^( )． 由于 是由程序运行之后得出的，所以当 它们不能达到指定的偏差值时．可增加 的 输人值以提高度计算精度．由图2可以看出， 过渡带的宽度4 大致为半个板值交错点阃的 宽度，所以，当M (一 )较大时．边缘是十 分陡峭的．

图5是笔者利用一个8 Hz的点阻FIR DF滤去某样品中的水溶剂峰之后所得到的实 际NMR谱图．由图可见．基本上消除了基线 漂移，使谱图更为精细．笔者的经验表明，宜采 用高、低通组合的方式实现点阻滤波+否剐将 由于过渡带的存在而难 实现窄带的点阻滤 波器． 实验中运行的程序是用C++语言编写 的．顺便指出，文献[2]给出的程序实际上还可 以运行差分器和希尔伯特变换器的设计程序．