基于幅频谱估计系统时域响应的改进方法研究

陈 翔 汪连栋 陈永光 曾勇虎 韩 慧 董 俊

（1.电子信息系统复杂电磁环境效应国家重点实验室，河南 洛阳471003；2.北京跟踪与通信技术研究所，北京100094）

引 言

现行标准及人们习惯中，一些电磁兼容测量结果常以幅频曲线表征，如材料的屏蔽效能、滤波器的插入损耗、接收天线的频域特性以及传感器探头的校正系数等等.在实验条件受限，而人们又希望仅通过频域测量结果估计系统的时域特性时，就需要采取一定的方法对忽略的相频信息进行重构.可行的方法是，假设系统为最小（或最大）相位系统，根据Hilbert变换重构其相位信息［1-2］.虽然不是所有系统都满足最小相位假设，但任何系统都可表示为一个最小相位系统和一个全通系统的级联，而全通系统并不会改变传输信号的幅度特性，因此可在最小相位假设条件下对系统进行初步考察.

最小相位法陆续被引入到相关研究领域解决了各自面临的实际问题，2000年，石立华等［3］提出使用最小相位系统模型，利用幅-频曲线重构系统的相位信息，进而估计系统的时域响应，找到了一条系统的频域指标到时域指标的转换途径.2004年，谢彦召等［4-5］根据最小相位原理，基于频域幅度谱数据成功反演了高空核爆电磁脉冲和余（正）弦阻尼振荡的时域波形，对于非最小相位信号，重建波形与原始波形有一定差异，但还是可以提供一些波形、累计能量、幅值量级的信息.2012年，曹景阳等［6］将传感器视为满足最小相位假设的信号传输系统，通过传感器的修正系数计算了相频特性，得到了系统的冲激响应，进而恢复了待测电磁脉冲的时域波形，解决了脉冲敏感度试验中门限电平不易确定的难题.

以上常规的最小相位估计时域响应方法，在实现过程中由于反复使用快速傅里叶变换（Fast Fourier Transform，FFT）在时-频域间转换信号，会产生额外的数值误差，此时必须仔细选择FFT的相关参数，减小引入误差.一些改进的处理方法是使用prony法［3］、最佳逼近元法［7］或者是滤波器建模［8］等对频域传递函数进行参数化估计，得到系统的离散传递函数模型，由频域参数模型估计系统的冲激响应.但是，prony法等参数建模方法也存在着一些的局限性［9］，如拟合精度差、所得极点不稳定、算法的收敛性特别依赖于初始值的选择等.文献［6］和文献［7］的讨论中还提到常规最小相位法对于高频噪声较为敏感，噪声对重建效果影响较大.

为克服prony法等参数估计方法的不足，减小最小相位法对噪声的敏感性，提出了一种将向量拟合法（Vector Fitting，VF）与最小相位法（Minimum-Phase，MP）相结合的由幅频谱估计时域响应的改进方法（简称VF-MP方法）.首先给出VF-MP方法的基本思路，然后在介绍最小相位法与向量拟合基本原理的基础上，以材料屏蔽效能测试为研究对象，使用材料屏蔽效能的仿真数据验证VF-MP方法的可行性.

1 基本原理

1.1 最小相位法

在假定系统满足最小相位的条件下，可通过Hilbert变换由系统的幅频信息恢复相频信息.利用FFT可以很方便地将数据在时域与频域间转换，也就可以使用Kolmogoroff方法（简称柯氏法）［10］求取幅度谱对应的最小相位序列，其基本流程如图1所示.其中，Xc（n）为实倒谱）复倒谱，x（n）为最小相位序列.首先，对幅度谱取自然对数，然后用傅里叶逆变换（Inverse Fast Fourier Transform，IFFT）和FFT实现Hilbert变换，再取指数就可得x（n）的傅里叶变换X（ω），经IFFT后得到x（n）.

图1 柯氏法求解最小相位序列

计算时需注意以下几点：1）取对数后的幅度谱在进行IFFT之前，需对数据进行对称延拓处理.2）使用FFT得到的实倒谱（K）是倒谱混叠后的偶数部分，需经频率不变线性滤波从Xc（n）恢复.3）计算的在N很大时才是复倒谱的有效近似.

式中，lmin（n）＝2u（n）-δ（n）为加窗滤波函数.

柯氏法采用FFT算法，实现简单、计算快速.但是当最小相位系统的零点比较靠近单位圆时，幅度谱的值接近于0，采用对数运算会产生大的误差，可使用改进的柯氏法［10］，通过迭代估计来减小误差.

1.2 向量拟合

如果直接使用图1中得到的最小相位序列作为系统的时域模型，最后估计的结果可能受到高频噪声的影响.较好的处理方法是根据最小相位法重构的复频域传递函数，选取一种模型对其进行参数估计.与一些常用的参数建模方法相比，B.Gustavsen提出的向量拟合法［11-12］，能够对频域测试数据进行稳定的有理函数逼近，已成功应用于电源系统、电路系统、天线等方面的建模中［13-15］.

复频域的离散数据，可以由给定阶次的两个多项式比值的有理函数式来拟合为

式中：a0，a1，…，aM和b0，b1，…，bN为待辨识的实系数；s＝jω.

假设有一组复频域的采样数据ωi、f（si），i＝1，…，N，将数据带入式（4）可以得到形如Ax＝b的线性方程组.当阶数较高、频率较宽时将产生一个“病态”的方程组.这对于频谱范围宽、频率上限可高达GHz以上的电磁脉冲估计来说，会产生严重的偏差，向量拟合法［11］通过将式（2）表示成式（3）的部分分式展开形式克服了这一缺点.

式中：cn是留数；an是实数或者共轭复数对的极点；d和h是实数.使用最小二乘法估计式（3）中的参数，分成极点辨识和留数辨识两步实现，

引入辅助函数σ（s）

或者写为

问题的求解就被转化为线性问题，其中cn，d，h，~cn未知.根据给定的一组频率点的采样值，可将式（6）整理成线性方程组的形式

x未知，使用最小二乘法求解式（7）.

由式（8）可得

一次求解得到的极点往往不够精确，将所得极点作为假定极点，重复式（4）到式（7）的步骤，进行迭代，得到更精确的极点.为了获得更为精确的结果，根据式（3）的展开式计算留数.以上关于极点和留数的辨识，都只涉及到线性方程组的运算，方程系数矩阵的条件数一般也不会随阶数的提高增大，与其它方法相比向量拟合法更适用于高阶或宽频带的有理函数拟合.

当曲线较为平滑时，仅使用实数形式的初始极点即可.当曲线有谐振峰时，初始极点应使用式（10）形式的复共轭极点形式，极点的频率应该覆盖测试的范围.

式中，α＝β／100.

1.3 VF-MP方法的基本流程

不考虑材料的磁饱和等非线性特性，将屏蔽体对电磁波的屏蔽近似为线性时不变系统.以电场屏蔽效能为研究对象，屏蔽体的电场屏蔽效能S（ω）与其传递函数H（ω）的关系可表示为式（11）.

式中：ω为角频率；Eu（ω）为未屏蔽的电场强度；Es（ω）为屏蔽后的电场强度.

传递函数的幅频特性可写为式（12）.

图2为VF-MP方法的流程图，即在使用最小相位法得到频域传递函数之后，采用向量拟合法得到频域传递函数的有理函数式参数模型，然后通过拉普拉斯逆变换得到屏蔽体的指数叠加形式的冲激响应函数，再与待估电磁脉冲直接在时域卷积，最终得到屏蔽体的电磁脉冲响应.

图2 采用VF-MP方法的参数模型估计流程图

使用向量拟合法对屏蔽体的电场传递函数H（jω）进行拟合后，可以将其表示为式（13）.

对式（13）进行拉普拉斯反变换得到屏蔽体冲激响应的双指数叠加形式.

由式（14）和时域卷积的定义，采用梯形积分法可得到x（t）和h（t）的递归卷积式（15）.

式中：Δt为采样间隔；tk时刻的响应由tk-1时刻的响应和两时刻激励共同确定.

递归卷积算法可以缩短计算时间，与向量拟合法相结合还可以减少使用FFT和IFFT在时频域间的转换次数，克服了一些常规参数模型估计的局限性，减小了数值误差.

2 平面材料仿真实验验证

以无限大平面材料的仿真数据为例，研究利用材料屏蔽效能的幅频曲线估计其电磁脉冲响应.设置理想平面材料的厚度为1mm，相对介电常数和相对磁导率为1，电导率为1×103S／m，材料的理论屏蔽效能曲线如图3（a）所示.图3（b）为满足IEC61000-4-4标准的入射电磁脉冲波形，使用双指数函数拟合得到，峰值场强为50kV／m，修正系数为1.03，脉冲前后沿参数分别为4×107s-1和6×108s-1，使用CST（Computer Simulation Technology）软件计算经材料衰减后的电磁脉冲波形（见图4）.

图3 建模用仿真数据

使用常规最小相位法计算图3（b）中透过材料后的脉冲波形，结果也列于图4中.可见最小相位法估计的电磁脉冲响应波形与理论计算的波形非常接近.理想无限大平面材料的屏蔽行为可看作一个良好的线性系统，在没有噪声的情况下，仅依靠常规最小相位法即可很好地估计材料的电磁脉冲响应.

图4 最小相位法估计的材料电磁脉冲响应

然而在实际测试中，数据会不可避免地遭受噪声污染，屏蔽效能曲线不会像图3（a）中那样平滑.通过给屏蔽效能曲线叠加上不同程度的高斯白噪声，来模拟实测频域曲线，研究噪声对最小相位法使用效果的影响.给图3（a）中的屏蔽效能曲线叠加上平均值为0、标准差分别为其最大值2%和4%的噪声.使用受噪声污染的屏蔽效能曲线，采用常规最小相位法估计材料的电磁脉冲响应，在同一噪声水平下一共计算了三次，结果如图5所示.可看出当噪声较大时，估计的时域波形与没有噪声污染时的波形相比，峰值误差增大，脉冲的上升沿部分受噪声的影响较小，但脉冲下降沿畸变严重；即使对同一水平的噪声，由于每次叠加噪声为随机生成，所估计波形也有很大不同.

图5 不同噪声水平下常规最小相位法估计的波形

仍以叠加噪声的屏蔽效能曲线为例，尝试使用VF-MP方法估计材料的电磁脉冲响应.由于理想无限大平面材料的屏蔽是一个良好的线性系统，所需模型的阶次较低，极点设置为实数即可.经试验发现，当拟合阶次为4时，进行5次迭代即可很好地拟合屏蔽效能传递函数曲线.使用4阶模型、迭代5次，分别对几组2%和4%噪声水平下的数据进行拟合，结果均保持了很好的稳定性，并未因为每次叠加噪声的随机性而发生大的变化，图6给出了2%噪声水平下传递函数幅度谱的拟合结果.可预见，估计的材料电磁脉冲响应波形也会保持很好的重复性，不会如图5中那样有剧烈变化.

图6 向量拟合法对传递函数幅频曲线的拟合

使用式（15）给出的递归卷积公式直接与图3（b）的时域波形卷积，即可得到材料的电磁脉冲响应，结果如图7所示.可见采用VF-MP法估计的材料电磁脉冲响应波形，与理论计算的波形非常接近，受噪声影响较小.

图7 不同噪声水平下VF-MP法估计的时域波形

3 暗箱窗口法仿真实验验证

继续使用暗箱窗口法的仿真数据来检验VFMP方法，测试材料的电磁参数与上一节保持不变.仿真模型使用文献［16］中的设置，箱体为1m×1m×1m的正方体，窗口尺寸为60cm×60cm，测试点A0位于窗口中心轴线上距窗口10cm处.图8（a）为测得的材料对于垂直入射平面波的插入损耗曲线；图8（b）为幅值1V／m、上升时间1ns、平顶保持47ns、下降时间4ns的方波脉冲激励下，窗口未加载材料时A0点的时域波形.

根据图8中的数据，使用VF-MP方法估计加载材料后A0点的时域波形，由于幅频曲线不平滑，应使用复共轭极点和较高阶次进行向量拟合.以均方误差来衡量向量拟合法对传递函数的拟合效果，分别尝试使用10阶、20阶、30阶、40阶模型，每一阶次迭代20次，计算相应阶次和迭代次数下的均方误差.经试验发现在阶次为30、迭代10次时，模型的均方误差较小，拟合精度已足够好.由VF-MP方法估计的时域波形如图9所示，可见VF-MP方法大致能够估计出加载材料之后的时域波形，但与实际响应波形相比，峰值偏小，拟合度不够高.使用VF-MP方法估计的材料峰值屏蔽效能为47.5dB，而根据CST软件计算的结果为44.6dB，相差了大约3dB.这主要是因为暗箱窗口法测试比无限大平面材料情况要复杂得多，不满足最小相位系统的假设.由于测试中窗口的存在，测得的脉冲不只受材料衰减的影响，窗口的截止波导作用同样会改变电磁脉冲各频率分量的相位，而且在箱体内还存在着些许谐振的影响，最终导致实际情况与最小相位假设不符.

图8 建模用仿真数据

图9 VF-MP法估计波形与仿真波形

4 结 论

采用向量拟合法与最小相位法相结合的VFMP方法是对常规最小相位法的有益补充，减小了因反复进行时频转换而产生的数值误差，估计结果受噪声影响小，重复性好.受最小相位假设的限制，对于不满足最小相位近似的系统，估计结果会存在一定的偏差，但可用作前期的粗略估计.

［1］HAYES M H，LIM J S，OPPENHEIM A V.Signal reconstruction from phase or magnitude［J］.IEEE Transactions on Acoustics，Speech，Signal Processing，1980，28：672-680.

［2］TESCHE F M.On the use of Hibert transform for processing measure CW data［J］.IEEE Transactions on Electromagnetic Compatibility，1992，34（3）：259-266.

［3］石立华，周璧华，陈 彬，等.基于幅-频曲线的系统时域响应特性评价方法［J］.电波科学学报，2000，15（4）：467-471.SHI Lihua，ZHOU Bihua，CHEN Bin，et al.Timedoamin characterization of a system based on the magnitude of its frequency response［J］.Chinese Journal of Radio Science，2000，15（4）：467-471.（in Chinese）

［4］谢彦召，王赞基，王群书，等.基于频域幅度谱数据重建电磁脉冲时域波形［J］.强激光与粒子束，2004，16（3）：320-324.XIE Yanzhao，WANG Zanji，WANG Qunshu，et al.Reconstruction of electromagnetic pulse waveform based on the amplitude spectrum data［J］.High Power Laser and Particle Beams，2004，16（3）：320-324.（in Chinese）

［5］谢彦召，刘顺坤，孙蓓云，等.电磁脉冲传感器的时域和频域标定方法及其等效性［J］.核电子学与探测技术，2004，24（4）：395-399.XIE Yanzhao，LIU Shunkun，SUN Beiyun，et al.Methodology of time-domain ＆frequency-domain calibration and equivalence for EMP sensor［J］.Nuclear Electronics ＆Detection Technology，2004，24（4）：395-399.（in Chinese）

［6］曹景阳，谢树果，苏东林.基于最小相位法重建电磁脉冲时域波形［J］.电波科学学报，2011，26（6）：1102-1106.CAO Jingyang，XIE Shuguo，SU Donglin.Application of minimum phase method in a pulse measurement［J］.Chinese Journal of Radio Science，2011，26（6）：1102-1106.（in Chinese）

［7］石立华，司 卿，张 琦，等.屏蔽室脉冲波屏蔽效能模型估计方法［J］.电波科学学报，2013，28（1）：125-129.SHI Lihua，SI Qing，ZHANG Qi，et al.Model estimation metohd for shielding effectiveness of shielding enclosure in pulse field［J］.Chinese Journal of Radio Science，2013，28（1）：125-129.（in Chinese）

［8］郭东义，石立华，郭曜华，等.基于频响函数测量脉冲磁场屏蔽效能的新方法［J］.信息与电子工程，2010，8（1）：49-52.GUO Dongyi，SHI Lihua，GUO Yaohua，et al.Measurement for shielding effectiveness of pulsed magnetic field based on frequency response function［J］.Information and Electronic Engineering，2010，8（1）：49-52.（in Chinese）

［9］张 黎，李庆民，王 伟，等.用于传递函数参数辨识的改进复向量拟合法［J］.高电压技术，2008，34（8）：1542-1546.ZHANG Li，LI Qingmin，WANG Wei，et al.Improved vector method for parameter identification of transfer functions［J］.High Voltage Engineering，2008，34（8）：1542-1546.（in Chinese）

［10］李衍达，常 迥.信号重构理论及其应用［M］.北京：清华大学出版社，1991：78-80.

［11］GUSTAVSEN B，SEMLYEN A.Rational approximation of frequency of frequency domain responses by vector fitting［J］.IEEE Transactions on Power Delivery，1999，14（3）：1052-1061.

［12］GUSTAVSEN B.Improving the pole relocating properties of vector fitting［J］.IEEE Transaction on Power Delivery，2006，21（3）：1587-1592.

［13］DESCHRIJVER D，MROZOWSKI M，DHAENE T，et al.Macromodeling of multiport systems using a fast implementation of the vector fitting method［J］.IEEE Microwave and Wireless Components Letters，2008，18（6）：383-385.

［14］郭 剑，邹 军，何金良，等.基于Vector Fitting的金属薄壳电磁脉冲屏蔽效能的计算［J］.清华大学学报：自然科学版，2004，44（10）：1302-1305.GUO Jian，ZOU Jun，HE Jinliang，et al.Evaluation of EMP shielding effectiveness for metallic sheels based on the vector fitting method［J］.J Tsinghua University：Science and Technology，2004，44（10）：1302-1305.（in Chinese）

［15］魏 琰，郭裕顺.基于向量拟合法的多端口网络函数有理逼近及其瞬态分析［J］.电路与系统学报，2008，13（1）：67-72.WEI Yan，GUO Yushun.Rational approximation of multiport network functions by the vector fitting algorithm［J］.Journal of Circuits and Systems，2008，13（1）：67-72.（in Chinese）

［16］陈 翔，陈永光，魏 明，等.基于屏蔽暗箱窗口法的材料电磁脉冲屏蔽效能时域测试［J］.高电压技术，2013，39（3）：668-674.CHEN Xiang，CHEN Yongguang，WEI Ming，et al.Time-domain test for electromagnetic pulse shielding effectiveness of materials based on the method of windowed semi-anechonic box［J］.High Voltage Engineering，2013，39（3）：668-674.（in Chinese）