无限精度下滤波器传递函数有多种等价实现结构，但在有限精度下，由于有限字长(finite word length，FWL)效应，滤波器系数以截断或舍入方式进行量化，其频响特性与理论值有偏差，因此滤波器实际性能达不到要求．然而滤波器实际性能与其实现结构有关系，不同结构具有不同的抗FWL效应能力[1-2]，灵敏度为衡量实现结构FWL效应指标之一．基于状态空间实现结构，文献[3]从系数误差灵敏度的角度优化滤波器结构，提出L2范数优化准则．文献[4]指出L2范数准则下的灵敏度最小化问题是非线性问题．针对非线性问题，文献[5-6]改变L2约束条件以优化灵敏度函数，得到了最佳实现结构，此结构的复杂度(非平凡参数个数)为N²+2N+1(N为滤波器阶数)，其复杂度随N增加而增加．针对复杂度问题，文献[7]提出复杂度为5N-1的滤波器结构，但灵敏度不如最佳实现结构．以上文献均基于实现结构优化无限冲激响应(infinite impulse response，IIR)滤波器的FWL效应．相比IIR滤波器，有限脉冲响应(FIR)滤波器具有稳定、易实现线性相位和适合多采样率转换等特点，被广泛应用于宽带无线通信系统领域．研究人员多采用算法优化FIR滤波器的设计，文献[8-11]利用启发式算法直接优化滤波器系数，改善滤波器在有限字长约束下的幅频和相频响应．针对FWL效应，文献[12]基于快速卷积算法提出一种低复杂度的并行FIR滤波器结构，节省了硬件资源但有限字长性能下降．文献[13]指出模拟退火法可优化FIR滤波器级联实现结构中的FWL效应，但优化效果不明显．为进一步降低FWL效应对滤波器性能的影响，文献[14]提出一种改进FIR滤波器结构，但其结构系数灵敏度不理想．以上有关FIR滤波器的文献对于复杂度的讨论没有给出具体表达式或列表，相应结论不直观，不便于比较和讨论．本研究根据模拟域与数字域变换理论，结合增量差分方法与双线性变换导出一种低灵敏度FIR数字滤波器的优化结构，以提高FIR滤波器的抗FWL效应能力．1 优化结构实现N阶线性时不变FIR数字滤波器传递函数为H(z)=f0zN+f1zN-1+⋯+fN-1z+fNa0zN+a1zN-1+⋯+aN-1z+aN, (1)式中：fi和ai（i=1，2，…，N）为H(z)分子和分母系数；a0=1,ai=0(∀i≠0)．增量差分Delta算子定义为δ=(z-1)/T，其中T为采样周期．δ的传递函数表示为H(δ)=β0δN+β1δN-1+⋯+βN-1δ+βNα0δN+α1δN-1+⋯+αN-1δ+αN，(2)由式(1)和(2)得到H(z)与H(δ)系数的对应关系为[β0,β1,⋯,βN]T=S[f0,f1,⋯,fN]T;[α0,α1,⋯,αN]T=S[a0,a1,⋯,aN]T,式中：βi和αi为H(δ)分子分母系数；S为一个下三角矩阵．将式(2)对应的状态空间实现(Aδ,bδ,cδ,dδ)表述为δx(n)=Aδx(n)+bδu(n);y(n)=cδx(n)+dδu(n),式中：x(n)为系统状态矢量；u(n)与y(n)为系统输入和输出．状态空间与传递函数的关系为H(δ)=dδ+cδ(δI-Aδ)-1bδ．(3)结合Delta算子的定义，将连续s域与离散δ域的双线性变换表示为s=δ/(Tδ+2), (4)则H(δ)转为H(s)的关系等式为H(s)=H(δ)δ=2s/(1-Ts)．(5)当采用状态空间模型描述H(s)时，记状态空间实现为(Φ,K,L,D)．根据式(4)与(5)可得：Φ=(2I+TAδ)-1Aδ;K=2(2I+TAδ)-1bδ;L=2cδ(2I+TAδ)-1;D=dδ-cδ(2I+TAδ)-1Tbδ. (6)与离散系统一样，连续时间状态空间实现满足H(s)=D+L(sI-Φ)-1K．(7)引入文献［15］中N阶模拟滤波器结构，其状态空间实现记为(Φin,Κin,Lin,Din)，满足式(7)，其中：Φin=0ε100⋯00-ε10ε20⋯000-ε20ε3⋯00⋯0000⋯0εN-10000⋯-εN-1-εN; (8)Kin=[0,0,0,⋯,2εN]T,式中εi∈R且大于零．下面计算参数εi．将式(7)中H(s)表示为H(s)=ψ(s)/Γ(s)，将H(s)的分母按幂次分解为Γ(s)=E(s)+O(s)，其中：E(s)表示偶次幂项和；O(s)表示奇次幂项和．定义Z(s)=O(s)/E(s)   (N为奇数);E(s)/O(s)   (N为偶数), 并将Z(s)作连分式展开，得到Z(s)=rNs+1rN-1s+1rN-2s+1⋮1r1s，式中ri为参数．则参数εi[15]为εi=1/(riri+1)   (i=1,2,⋯,N-1);1/ri   (i=N). (9)根据上式即可确定系数矩阵Φin与Kin．假定(Φc,Kc,Lc,Dc)为H(s)的任意可控状态空间实现，则(Φc,Kc,Lc,Dc)到(Φin,Κin,Lin,Din)的相似变换矩阵为Jin=ΛcΛin-1，其中：Λc=[Kc,ΦcKc,⋯,ΦcN-1Kc];Λin=[Kin,ΦinKin,⋯ ,ΦinN-1Kin].根据Lin=LcJin求出系数矩阵Lin．再由式(6)求得模拟滤波器结构离散化后对应的δ域状态空间实现结构，将其称为δin结构，其状态空间实现记为(Aδin,bδin,cδin,dδin)，表示如下：Aδin=2Φin(1-TΦin)-1;bδin=(2I+TAδin)Kin/2;cδin=2Lin(2I+TAδin)/2;dδin=Din+cδin(2I+TAδ)-1Tbδin. (10)由式(8)求得Aδin中(1-TΦin)-1是全参数化的，因此对矩阵Aδin中(1-TΦin)-1进行稀疏化分解．首先定义N0=2(N-1),φi+1=-Tεi+1χi，χi=1/(1-Tεiφi)   (i=1,2,⋯,N-2)．其中：φ1=-Tε1；χN-1=1/(1+TεN-TεN-1φN-1)．由于I-TΦin=1-Tε100⋯0Tε11-Tε20⋯00Tε2 1-Tε3⋯0⋯00000⋯-TεN-1000⋯TεN-11+TεN，若定义XN0=μ(1,2,Tε1)μ(2,2,χ1)，μ(i,j,x)为除了第(i,j)个元素为x的N×N维单位矩阵，则有    (I-TΦin)XN0=1000⋯0-φ11-Tε20⋯00Tε2 1-Tε3⋯0⋯0000⋯0-TεN-1000⋯TεN-11+TεN，若再定义XN0+1-i=μ(i,i+1,Tεi)μ(i+1,i+1,χi), (11)式中i=1,2,⋯N-1．XN0+1-i为单位矩阵，除了XN0+1-i(i:i+1,i:i+1)=1υi0χi，(12)其中υi=Tεiχi．经过计算得到    (I-TΦin)XN0XN0-1⋯XN0+1-i⋯XN=1000⋯0-φ1100⋯00-φ210⋯0⋯00000⋯0000⋯-φN-11=g, (13)记矩阵g-1=XN-1XN-2⋯X1，其中Xi定义为Xi=μ(i+1,i,φi)    (i=1,2,⋯,N-1)．(14)结合式(13)有(I-TΦin)-1=XN0XN0-1⋯X1=∏k=1N0Xk，当1≤k≤N-1时，Xk为式(14)；当N≤k≤N0，Xk为式(11)．综上得到FIR滤波器的稀疏优化实现结构，记为δnew结构，将其状态空间方程表示为x(n+1)=2Φinx(n)+Yu(n);y(n)=cδinx(n)+dδinu(n), (15)式中Y=2Kin为含一个非平凡参数的稀疏矩阵．2 灵敏度表达式推导将滤波器结构中的参数用集合Ωτ={τk}表示，当τk用有限位二进制表示时，实际τ˜k与理想τk有量化误差Δτk，即τ˜k=τk+Δτk．Δτk使实际系统H˜(δ)与理想系统H(δ)存在偏差ΔH(δ)，    ΔH(δ)=H(δ)-H˜(δ)≈∑τk∈Ωτ,∀k(∂H(δ)/∂τk)Δτk, (16)式中∂H(δ)/∂τk为滤波器传递函数对第k个τk的灵敏度函数．记∂H(δ)/∂τk=Mτk(δ)，M为滤波器的传递函数H(δ)关于τk的灵敏度，M=∑τk∈Ωτ,∀k||Mτk(δ)||22，(17)式中Mτk(δ)的状态空间实现(Aδ,bδ,cδ,dδ)满足式(7)，由留数定理得到||Mτ(δ)||22=tr(dδdδT+cδWcδcδT)，(18)式中W0δ和Wcδ为可观和可控格莱姆矩阵．在δnew结构中，参数包括{υi,χi,φi}，εN及cδin,dδin中的参数．若σ表示δnew结构中的任意一个结构参数，则传递函数H(δ)对σ的灵敏度函数为    Mσ(δ)=cδin(δI-Aδin)-1∂Aδin∂σ(δI-Aδin)-1∙bδin+∂cδin∂σ(δI-Aδin)-1bδin+∂dδin∂σ+cδin(δI-Aδin)-1∂bδin∂σ. (19)首先计算Mcδin(δ)与Mdδin(δ)．矩阵cδin与dδin相互独立，则由式(19)可得MCδin(δ)=γiT(δI-Aδin)-1bδin,MDδin(δ)=1, (20)式中γi为1×N维列向量，除第i个元素为1，其余元素均为0．根据式(18)得到||Mcδin(δ)||22=tr(γiTWcinγi)；||Mdδin(δ)||22=1，(21)式中Wcin=I是δin结构的可控格莱姆矩阵．然后，计算Mφi(δ)．参数φi与矩阵cδin,dδin无关，则由式(19)可得    Mφi(δ)=cδin(δI-Aδin)-1∂Aδin∂φi(δI-Aδin)-1bδin+cδin(δI-Aδin)-1∂bδin∂φi. (22)由式(14)知，只有Xi=μ(i+1,i,φi)中含φi，则∂Aδin∂φi=2ΦinKiγi+1γiTFi；∂bδin∂φi=Kiγi+1γiTFiY，(23)式中：Ki=∏k=i+1N0Xk    (i=1,2,⋯,N-1)；Fi=I    (i=1),∏k=1i-1Xk    (i=2,⋯,N-1).因此可以得到    Mφi(δ)=cδin(δI-Aδin)-1[2ΦinKiγi+1γiTFi*(δI-Aδin)-1bδin+Kiγi+1γiTFiY]. (24)假定ξ(δ)=[c2(δI2-A2)-1b2+d2][c1(δI1-A1)-1b1+d1]， 其中I1和I2为适当维度单位矩阵．若上式中A1=Aδin，b1=bδin,c1=2ΦinEiγi+1γiTFi，d1=Eiγi+1γiTFiY，A2=Aδin，c2=cδin，d2=0，则式(24)又可以表示为    Mφi(δ)=(d2c1, c2 )δI100I2-A10b2c1A2-1∙(b1,b2d1)T+d2d1 =c¯(δI-A¯)-1b¯+d¯. (25)利用式(18)即可求得||Mφi(δ)||22，接着计算Mε¯i(δ)．由式(12)知υi只存在于XN0+1-i中，则可得∂Aδin/∂υi=2ΦinQiγiγi+1TPi;∂bδin/∂υi=Qiγiγi+1TPiJ，(26)式中：Pi=∏k=1N0-iXk    (i=1,2,⋯,N-1);Qi=I    (i=1),∏k=N0+2-iN0Xk    (i=2,⋯,N-1),因此    Mυi(δ)=cδin(δI-Aδin)-1Qiγiγi+1TPiJ+2ΦinQiγiγi+1TPi*(δI-Aδin)-1bδin.取A1=Aδin，b1=bδin，C1=2ΦinQiγiγi+1TPi，D1=Qiγiγi+1TPiY，A2=Aδin,b2=I,c2=cδin,d2=0，由式(18)得||Mε¯i(δ)||22．接着，计算Mχi(δ)．由式(11)可知χi也只存在于XN0+1-i(i=1,2,⋯N-1)中，可求得    Mχi(δ)=cδin(δI-Aδin)-1[Qiγi+1γi+1TPiY+2ΦinQiγi+1γi+1TPi(δI-Aδin)-1bδin],取A1=Aδin，b1=bδin，c1=2ΦinQiγi+1γi+1TPi,d1=Qiγi+1γi+1TPiY，A2=Aδin,b2=I,c2=cδin,d2=0，同理得||Mχi(δ)||22．最后计算MεN(δ)．因为bδin=(I-TAδin)-1Y，Y=2εNγN=υNγN，所以有MεN(δ)=cδin(δI-Aδin)-1(I-TΦin)-1γN，根据式(18)得到||MεN(δ)||22=tr[γNT(I-TΦin)-TWoin(I-TΦin)-1γN]，式中Woin为δin结构的可观格莱姆矩阵．综上，δnew结构的总灵敏度为    Mnew=∑i=1N||Mcδin(δ)||22+||Mdδin(δ)||+||MεN(δ)||22+∑i=1N-1(||Mυi(δ)||22+||Mχi(δ)||22+||Mφi(δ)||22). (27)3 舍入噪声增益理论分析在无限精度下，δnew结构输入输出方程表示为：x(0)(n)=2Φinx(n)+Yu(n);x(k)(n)=Xkx(k-1)(n);x(n+1)=x(N0)(n);y(n)=cδinx(n)+dδinu(n). (28)当定点制实现时，信号s(n)与非平凡参数τk相乘后经过量化器引入舍入噪声es(n)(舍入噪声增益G衡量)，将其建模成均值为0、方差为σc2的高斯白噪声，其定义为es(n)=Q[τks(n)]-τks(n)．设Δy(n)表示实际输出值与理想输出值偏差，G定义为G=E[(Δy(n))2]E[(es(n))2]=tr(dδinTdδin+bδinTWoinbδin)．(29)在δnew结构中，非平凡参数包括{υi,χi,φi}，εN及cδin,dδin中的参数．首先，考虑非平凡参数εN引入的舍入噪声eεN(n)，式(28)变为：x(0)(n)=2Φinx(n)+Yu(n);x(k)(n)=Xkx(k-1)(n);x(n+1)=x(N0)(n);y(n)=cδinx(n)+dδinu(n). (30)将式(30)与式(28)相减，并定义：Δx(n)=x¯(n)-x(n);Δx(k)(n)=x¯(k)(n)-x(k)(n)Δy(n)=y¯(n)-y(n).；(31)结合式(28)、(30)及(31)可得：Δx(0)(n)=2ΦinΔx(n)+eεN(n);Δx(k)(n)=XkΔx(k-1)(n);Δx(n+1)=Δx(N0)(n);Δy(n)=cδinΔx(n), (32)式中eεN(n)=Q[Yu(n)]-Yu(n)．则状态空间方程为Δx(n+1)=AδinΔx(n)+γNeεN(n);Δy(n)=cδinΔx(n). (33)对比式(28)，eεN(n)为输入，与Δy(n)的传递函数为HεN(δ)，则HεN(δ)=cδin(δI-Aδin)-1(I-TΦin)-1γN，由式(29)求得舍入噪声增益为GεN=tr(γNT(I-TΦin)-TWoin(I-TΦin)-1γN). (34)然后，考虑非平凡参数{υi,χi,φi}．由以上分析可知{υi,χi,φi}均存在于Xi中，则式(30)变为：x¯(0)(n)=2Φinx¯(n)+Yu(n) ;x¯(k)(n)=Xkx¯(k-1)(n) ;x¯(n+1)=x¯(N0)(n) ;y(n)=cδinx¯(n)+dδinu(n). (35)将式(35)与式(28)相减，并结合式(31)，定义ei(n)=Q[Xix(i-1)(n)]-Xix(i-1)(n)，可得：Δx(0)(n)=2ΦinΔx(n);Δx(k)(n)=XkΔx(k-1)(n)   (∀k≠i);Δx(i)(n)=XiΔx(i-1)(n)+ei(n);Δx(n+1)=Δx(N0)(n);Δy(n)=cδinΔx(n).则状态空间方程表示为Δx(n+1)=AδinΔx(n)+Viei(n);Δy(n)=cδinΔx(n),式中：ei(n)=     ei+1eφi(n)    (1≤i≤N-1),    ekeυk(n)+ek+1eχk(n)    (i=N0+1-k,1≤k≤N-1);Vi=I    (i=N0),∏k=i+1N0Xk    (i≠N0).由结构参数{υi,χi,φi}引入的舍入噪声与Δy(n)的传递函数Hi(δ)=cδin(δI-Aδin)-1Vi，由式（29）求得舍入噪声增益为Gi=tr(ViTWoinVi)．最后考虑cδin和dδin的舍入噪声增益，得到Gcδin=ς1;Gdδin=ς2，式中ς1与ς2分别为cδin和dδin中非平凡参数的个数．综上，优化结构的总舍入噪声增益为Gδnew=∑i=1N0Gi+ς1+ς2．4 仿真实例及分析4.1　复杂度对比以16阶低通FIR数字滤波器为例，仿真参数：通带截止频率0.2π，阻带截止频率0.4π，通带最大衰减1 dB，阻带最小衰减40 dB．为了更简便对比实现复杂度，采用非平凡参数个数来衡量．对于N阶FIR滤波器，z算子结构的复杂度为2N+1，δDFIIt结构[13]为3N+1，δnew结构为4N-1．显然，所提δnew结构复杂度高于前两种结构．4.2　灵敏度与舍入噪声增益分析表1为不同采样周期下，Mz=17和Gz=25.948 6时各实现结构系数灵敏度和舍入噪声增益对比，其中z算子结构与δDFIIt结构参数灵敏度与舍入噪声增益分别由文献[16]与文献[3]定义．δnew结构灵敏度由式(27)求得，舍入噪声增益由式(39)求得．10.13245/j.hust.210110.T001表1各实现结构系数灵敏度和舍入噪声增益对比TMz=17Gz=25.948 6MδDFIItMδnewGδDFIItGδnew1×1002.093 0×10854.486 06.278 9×1081.067 3×1021×10-11.555 2×1061.045 34.665 6×10618.955 31×10-21.551 3×1041.008 94.653 8×10410.951 81×10-3156.120 01.001 0467.350 010.921 61×10-42.551 21.000 06.653 510.926 01×10-51.015 51.000 02.640 510.926 0由表1可知：随着T减小，MδDFIIt，Mδnew，GδDFIIt及Gδnew逐渐减小．当T=1时，所提δnew结构较z算子结构在灵敏度与舍入噪声性能上没有优势，但其数值远小于δDFIIt结构；当T≤1×10-1时，MδnewMz，GδnewGz，而δDFIIt结构只有当T≤1×10-4时，MδDFIItMz，GδDFIItGz．原因是高速采样下，Delta算子比z算子有更好数值特性．另外，当T≤1×10-4时，GδDFIItGδnew，这表明δnew结构在高速采样下系数灵敏度的减小是以增加舍入噪声增益为代价．舍入噪声增益的增加并不影响滤波器的频响特性，下面通过仿真验证该结论．4.3　频响特性分析为验证上述结论并分析各结构抵抗FWL效应的能力，将系数小数部分的字长量化为Bc．当T=1×10-4，Bc=6，8与12 bit时，z算子结构、δDFIIt结构及δnew结构得到的幅频响应A(w)与相频响应Φ(w)如图1~3所示．10.13245/j.hust.210110.F001图1频率响应特性比较(T=1×10-4，Bc=6 bit)10.13245/j.hust.210110.F002图2频率响应特性比较(T=1×10-4，Bc=8 bit)10.13245/j.hust.210110.F003图3频率响应特性比较(T=1×10-4，Bc=12 bit)从图1~3中观察到：当T=1×10-4时，有限字长量化对z算子结构与δDFIIt结构的频响特性均有影响，而δnew结构的幅频和相频响应更稳定，即与理想频响特性曲线趋于重合．随着Bc从6 bit增加到12 bit，z算子结构与δDFIIt结构频响特性曲线偏离理想频响特性曲线的程度减小．此外，z算子结构与δDFIIt结构的幅频和相频响应接近理想频响特性所需的字长分别为8与10 bit，而δnew结构只需要z算子结构的1/2．因为所提δnew结构的字长被充分利用，降低了量化对频响特性影响．这表明所提δnew结构更节约系统资源．当T=1×10-3，1×10-5，Bc=6，8 bit时，z算子结构、δDFIIt结构与δnew结构的幅频响应A(w)与相频响应Φ(w)如图4~7所示．10.13245/j.hust.210110.F004图4频率响应特性比较(T=1×10-3，Bc=6 bit)10.13245/j.hust.210110.F005图5频率响应特性比较(T=1×10-3，Bc=8 bit)10.13245/j.hust.210110.F006图6频率响应特性比较(T=10-5，Bc=6 bit)10.13245/j.hust.210110.F007图7频率响应特性比较(T=10-5，Bc=8bit)由图4和5可知：当T=1×10-3时，量化对z算子结构和δDFIIt结构的频响特性均有影响，但δDFIIt结构曲线偏离更严重，这是因为Delta算子在高速采样下有较好数值特性．在图6与图7中，当T=1×10-5时，δDFIIt结构的频响特性更接近理想频响特性，而δnew结构始终与理想曲线基本符合．这是因为δnew结构具有L2缩放特性，使得结构参数在量化后的误差降低，从而使得频响特性偏差减小．这表明δnew结构比δDFIIt结构的频响特性更优，而且在高速采样时Delta算子比z算子有更好数值特性．同时，结合图1与图2可以发现：T变化不会影z算子结构频响特性，仅有量化字长对其影响更明显．表2列出了z算子、δDFIIt结构及δnew结构所得到的幅频响应与理想幅频响应在各频率点的数值差(取平均值)ϕz，ϕδDFIIt与ϕδnew．10.13245/j.hust.210110.T002表2各结构幅频响应与理想幅频响应的差值Bc/bitTϕz/dBϕδDFIIt/dBϕδnew/dB61×10-36.468 88.301 90.396 21×10-46.468 81.157 10.178 01×10-56.468 80.545 80.172 181×10-32.414 54.544 80.109 91×10-42.414 51.157 00.163 81×10-52.414 50.194 70.172 6由表2可知：当Bc=6 bit时，随着T减小，ϕz值不变，而ϕδDFIIt与ϕδnew减小，表明δDFIIt结构与δnew结构幅频响应与理想幅频响应的程度减小．当Bc=8 bit时，ϕz不变，与z算子结构和δDFIIt结构相比，所提δnew结构与理想幅频响应的差值最小，表明所提结构的幅频响应更趋近于理想幅频响应．综上分析结果可知：采样周期的变化对z算子的灵敏度、舍入噪声增益和频响特性无影响．当考虑有限字长约束时，δnew结构的FIR滤波器抗FWL效应能力优于z算子与δDFIIt结构，且δnew结构的频响特性接近理想频响特性时所需的量化字长更少，更节省系统资源．5 结语从实现结构入手改善FIR滤波器有限精度实现时FWL效应影响，提出一种低灵敏度FIR滤波器优化结构．对比分析了基于传统移位z算子结构、增量差分Delta算子结构和所提优化结构的灵敏度、舍入噪声性能及频响特性．实验结果表明：与z算子和Delta算子结构相比，优化结构的系数灵敏度随采样频率的增加而减小，且在高速采样下系数灵敏度的减少是以增加舍入噪声增益为代价；在系数有限精度表示下所提优化结构的频响特性更优，具有更好抗FWL效应的能力．