CN102707143A - 采用灰色多周期模型提取电磁频谱中谐波分量的方法 - Google Patents

Abstract

本发明公开了一种采用灰色多周期模型提取电磁频谱中谐波分量的方法，首先采集干扰设备的频谱数据；然后利用灰色系统模型拟合出频谱数据的总体趋势分量；然后利用周期模型对频谱数据去除总体趋势分量后的数据提取谐波分量，获得谐波的量化表达式；最后将谐波量化表达式中的谐波频率与干扰设备中电路元器件的工作频率进行比对，若与某个元器件的工作频率相等则可确定对应谐波分量产生的根源。本发明方法将频谱数据中谐波干扰的特性用数学函数式表达出来，物理特征明确，便于电磁干扰频谱数据部分特性的认知，起到指导干扰设备电磁兼容性分析和设计的作用。

Description

采用灰色多周期模型提取电磁频谱中谐波分量的方法

技术领域

本发明涉及一种从电磁干扰发射频谱中提取出谐波分量的方法，属于电磁兼容技术领域。

背景技术

随着科学技术的发展，电气电子设备数量越来越多，这些设备在工作时通常会产生一些有用或无用的电磁发射，这些发射可能会对其他设备造成潜在电磁干扰。这些电磁干扰会造成设备的性能降级，甚至可能造成设备的毁坏。因此要对设备的电磁干扰进行控制。首先需要查找电磁干扰源，然后提取其电磁特征进行分析，而后采取相应的电磁兼容性设计。目前查找干扰源的方法是首先对发射设备采取电磁兼容检测试验，根据国军标(GJB151A)《军用设备和分系统电磁发射和敏感度要求》，电磁发射和敏感度要求与极限中5.3.15.2极限规定：“电场辐射发射不应超过图RE102-1～RE102-3所示的值。在30MHz以上，水平极化场和垂直极化场均应满足这些极限要求。”来判断设备是否满足电磁发射要求，一旦发射超过限值，则采用经验法或者近场探头检测法来判断产生干扰的原因。经验法对检测人员的专业技能要求较高，而近场探头检测法是一种非常耗时的方法，通常还不能追溯电磁干扰产生的根源。

在电磁干扰测试中发现，干扰频谱虽然复杂，一些元器件产生的具有稳定重复频率的窄带干扰信号在被测频谱中占据主导地位，这些信号一般是元器件的非线性产生的谐波分量，分析该类谐波分量特性和来源具有重要意义。

发明内容

本发明针对电磁干扰发射频谱中的主要贡献者——谐波分量，提出了一种采取灰色多周期模型对干扰设备发射频谱数据分析的方法。首先采集干扰设备的频谱数据；然后利用灰色系统模型拟合出频谱数据的总体趋势分量；然后利用周期模型对频谱数据去除总体趋势分量后的数据提取谐波分量，获得谐波的量化表达式；最后将谐波量化表达式中的谐波频率与干扰设备中电路元器件的工作频率进行比对，若与某个元器件的工作频率相等则可确定对应谐波分量产生的根源，为干扰的分析提供理论依据，从而指导干扰设备的电磁兼容性设计与整改。

本发明采用灰色多周期模型提取电磁频谱中谐波分量的方法包括有下列步骤：

步骤一：对频谱进行分段

采用斜率一致性原则对频率AF为横坐标，场强值AM为纵坐标绘制出的频谱幅度图中的频谱曲线进行分段，得到互不交叠的频谱频段F＝{f₁，f₂，…，f_w}；

步骤二：选取超标频段

将所述的频谱频段F＝{f₁，f₂，…，f_w}与国军标GJB151A中规定的发射极限进行对比，提取出超标频段FX＝{Xf₁，Xf₂，…，Xf_u}；

步骤三：构建总体趋势序列

任意一个超标频段Xf_u中包括有多个频率、以及每个频率对应的场强值，将所述的任意一个超标频段Xf_u用频谱场强序列表示为 $X_{\mathrm{test}}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(n\right)\};$

步骤301：对测试得到的任意一个超标频段Xf_u上的频谱场强序列

进行一次累加，生成一次累加生成序列 $X_{\mathrm{AGO}}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{AGO}}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{AGO}}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{AGO}}^{{\mathrm{Xf}}_u}\left(n\right)\};$

步骤302：对所述的一次累加生成序列X_AGO(Xf_u)进行紧邻均值生成，得到紧邻均值生成序列 $Z_{\mathrm{avg}}\left({\mathrm{Xf}}_u\right)=\{z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(1\right),z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(2\right),...,z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(n\right)\};$

步骤303：令

为GM(1，1)模型的原始形式，其中a表示发展系数的负数，b为灰色作用量；将a和b用参数列表示为AA＝[a，b]^T，T表示坐标转置。采取最小二乘法对参数列AA＝[a，b]^T进行求解得到AA＝(Z_add ^T(Xf_u)×Z_add(Xf_u))^-1×Z_add ^T(Xf_u)×X_extract(Xf_u)，因此可求得灰色GM(1，1)模型的时间响应序列为

经过一次累减还原得到总体趋势序列 $X_{\mathrm{trend}}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{trend}}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{trend}}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{trend}}^{{\mathrm{Xf}}_u}\left(n\right)\};$

步骤四：提取超标频谱中的谐波分量

步骤401：计算频率上的场强值与总体趋势序列在频率上的模拟值之差，得到第一次残差序列 $X_{\mathrm{diff}1}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(n\right)\};$

步骤402：对第一次残差 $X_{\mathrm{diff}1}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(n\right)\}$ 采用方差分析法进行优势周期的提取，并将提取出来的优势周期记为m；然后采用均值生成公式

计算序列X_diff1(Xf_u)中每个元素的均值生成函数；对均值生成函数

按优势周期m作周期性延拓，得到周期延拓函数

最后将GM(1，1)模型和周期延拓函数

叠加，得到灰色单周期模型序列 $X_{\mathrm{gp}1}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{gp}1}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{gp}1}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{gp}1}^{{\mathrm{Xf}}_u}\left(n\right)\};$

步骤403：计算曲线中的其他谐波分量，求灰色单周期模型序列X_gp1(Xf_u)与频谱场强序列X_test(Xf_u)的第二次残差序列 $X_{\mathrm{diff}2}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{diff}2}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{diff}2}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{diff}2}^{{\mathrm{Xf}}_u}\left(n\right)\};$

重复步骤三中对残差进行周期性建模的方法，再次提取残差序列X_diff2的优势周期m₂；如果m2_与m相等，则判定频谱场强序列X_test无其他谐波分量，建模工作结束，得到灰色双周期序列X_gp2；

采用步骤402中对残差进行周期性建模的方法，对第二次残差序列X_diff2(Xf_u)进行周期建模，分别得到第二次残差序列X_diff2(Xf_u)的优势周期记为m₂、周期延拓函数

灰色双周期模型序列 $X_{\mathrm{gp}2}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{gp}2}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{gp}2}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{gp}2}^{{\mathrm{Xf}}_u}\left(n\right)\};$

步骤404：若m₂与m相同，则谐波分量的提取结束；

若m₂与m不相同，则继续进行第三次残差序列的周期性建模，直至与上一次提取的优势周期相同，则结束谐波分量的提取。

最后将具有总体趋势的GM(1，1)模型序列与所有具有周期趋势的谐波分量叠加，得到任意一个超标频段Xf_u的灰色多周期拟合模型序列

且任意一个频率点上的灰色多周期模型拟合值 $x_{\mathrm{gps}}^{{\mathrm{Xf}}_u}\left(n\right)=x_{\mathrm{trend}}^{{\mathrm{Xf}}_u}\left(n\right)+E_m^{{\mathrm{Xf}}_u}\left(n\right)+E_{m2}^{{\mathrm{Xf}}_u}\left(n\right)+...+E_{\mathrm{ms}}^{{\mathrm{Xf}}_u}\left(n\right);$

步骤五：谐波分量与器件频率的匹配确定干扰源

将通过步骤四分析得出的频谱曲线表达式 $x_{\mathrm{gps}}^{{\mathrm{Xf}}_u}\left(n\right)=x_{\mathrm{trend}}^{{\mathrm{Xf}}_u}\left(n\right)+E_m^{{\mathrm{Xf}}_u}\left(n\right)+E_{m2}^{{\mathrm{Xf}}_u}\left(n\right)+...+E_{\mathrm{ms}}^{{\mathrm{Xf}}_u}\left(n\right)$ 中的谐波分量

的谐波的周期m、m₂、…m_s与干扰设备中所含的元器件的频率Freq进行比对；

(A)首先选取谐波周期m进行比对，若存在某个元器件的频率与其相等(即m＝Freq)，则认为该谐波分量

是由此元器件产生的干扰，确定该元器件为一个干扰源；

(B)若不存在元器件的频率与m相等(即m≠Freq)，则需进一步分析干扰设备中可能产生此频率间隔的元器件。

按(A)和(B)的对比法找出每个周期分量对应的的各个元器件，并按其对应的谐波提取顺序，得到具有不同间隔周期的谐波分量对整个干扰发射频谱曲线的贡献大小，按照所述谐波提取顺序依次采取干扰抑制措施，直至干扰设备的电磁干扰测试满足军标要求。

本发明方法的优点在于：

①通过电磁干扰测试手段得到干扰设备的频谱数据作为仿真数据来源，并从频谱数据中利用灰色系统模型拟合出频谱数据的总体趋势分量。该总体趋势分量反映了设备的宽带干扰与环境电平的叠加，作为提取谐波分量的基础。

②针对干扰设备的频谱数据去除总体趋势分量后的数据，利用周期模型提取出在干扰发射频谱中占主导地位的谐波分量，有针对性的简化频谱分析的过程。

③使用Matlab仿真软件，读入测试频谱数据，进行仿真，直接输出频谱数据中包含的所有谐波分量的频率，与干扰设备中所含元器件的频率吻合即可找出干扰源，使干扰源查找更加方便快捷。

④对频谱数据中的谐波分量进行量化表征，物理意义明确，对电磁干扰频谱部分特性进行认知，起到指导设备电磁兼容性的分析和设计的作用。

附图说明

图1是传统电磁干扰测试设备配置图。

图2是采用传统电磁干扰测试方法得到的频谱发射曲线与GJB151A国军标极限值对比图。

图3是采用本发明方法得到的频谱场强序列的一次累加生成曲线。

图4是采用本发明方法得到的频谱场强序列与GM(1，1)模型曲线对比图。

图5是采用本发明方法得到的频谱场强序列与灰色单周期模型曲线对比图。

图6是采用本发明方法得到的频谱场强序列与两次周期建模曲线对比图。

具体实施方式

下面将结合附图和实施例对本发明做进一步的详细说明。

参见图1所示，电磁干扰测试所需设备包括有接收天线、测量接收机和数据记录装置(计算机)，测量接收机分别与接收天线和计算机采用有线方式连接；接收天线通过无线方式检测干扰设备的辐射发射量。

在本发明中，不同的频段采用不同的接收天线。在频段10kHz～30MHz范围内使用具有阻抗匹配网络的拉杆天线，在频段30～200MHz范围内使用双锥天线，在频段200MHz～18GHz范围内采用双脊喇叭天线。在本发明中，使测量接收机在对应天线的接收频段范围内扫描，得到干扰设备产生辐射的频率(单位Hz)，以及所述频率上对应的场强值(单位dBμV/m)。对测试接收机接收到的信息采用集合形式表示为A_in＝{AF，AM}，AF表示干扰设备辐射干扰的频率，AM表示在该频率AF上对应的场强值(简称为场强值)。

本发明中的数据记录装置包括有一台计算机以及运行于计算机内的关联软件，所述关联软件包括有电磁兼容性检测软件EMC2003、Matlab 7.5.0软件。计算机是一种能够按照事先存储的程序，自动、高速地进行大量数值计算和各种信息处理的现代化智能电子设备。最低配置为CPU 2GHz，内存2GB，硬盘180GB；操作系统为windows 2000/2003/XP。在本发明中使用的计算机中运行的电磁兼容性检测软件EMC2003，将计算机与测试接收机相连，并启动软件，则可以记录测试中获取的数据信息A_in。软件会自动以AF为横坐标(频率，单位Hz)，AM为纵坐标(场强值，单位dBμV/m)，绘制对应的频谱幅度图，如图2所示。

在图2中，实线代表干扰设备测试得到的频谱场强曲线，虚线代表国军标GJB151A中的RE102的极限值。实线超过虚线的部分即为频谱超标范围。

本发明采用灰色多周期模型提取电磁频谱中谐波分量的方法包括有下列步骤：

步骤一：对频谱进行分段

采用斜率一致性原则对频率AF为横坐标，场强值AM为纵坐标绘制出的频谱幅度图中的频谱曲线进行分段，得到互不交叠的多个频段F＝{f₁，f₂，…，f_w}，简称为频谱频段F＝{f₁，f₂，…，f_w}；

所述的频谱频段F＝{f₁，f₂，…，f_w}中w表示频段标识号，f₁表示第1个频段、f₂表示第2个频段、f_w表示第w个频段，也称为任意一个频段。

步骤二：选取超标频段

将所述的频谱频段F＝{f₁，f₂，…，f_w}与国军标GJB151A中规定的发射极限进行对比，提取出超出军标极限的频段FX＝{Xf₁，Xf₂，…，Xf_u}，简称为超标频段FX＝{Xf₁，Xf₂，…，Xf_u}；

所述的超标频段FX＝{Xf₁，Xf₂，…，Xf_u}中u表示超标频段标识号，Xf₁表示第1个超标频段、Xf₂表示第2个超标频段、Xf_u表示第u个超标频段，也称为任意一个超标频段。

在本发明中，任意一个超标频段Xf_u中包括有多个频率、以及每个频率对应的场强值，将所述的任意一个超标频段Xf_u用频谱场强序列表示为 $X_{\mathrm{test}}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(n\right)\};$

所述的频谱场强序列 $X_{\mathrm{test}}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(n\right)\}$ 中表示在任意一个超标频段Xf_u下的第1个频率上的场强值、

表示任意一个超标频段Xf_u下的第2个频率上的场强值、

表示任意一个超标频段Xf_u下的第n个频率上的场强值，n表示任意一个超标频段Xf_u下的频率标识号。

同理可得，第1个超标频段Xf₁的频谱场强序列

第2个超标频段Xf₂的频谱场强序列 $X_{\mathrm{test}}\left({\mathrm{Xf}}_2\right)=\{x_{\mathrm{test}}^{{\mathrm{Xf}}_2}\left(1\right),x_{\mathrm{test}}^{{\mathrm{Xf}}_2}\left(2\right),...,x_{\mathrm{test}}^{{\mathrm{Xf}}_2}\left(n\right)\}.$

为了便于理解，下面以任意一个超标频段Xf_u进行列举说明。

步骤三：构建总体趋势序列

在本发明中，从任意一个超标频段Xf_u中提取出总体趋势序列X_trend(Xf_u)的步骤如下：

步骤301：对测试得到的任意一个超标频段Xf_u上的频谱场强序列进行一次累加，生成一次累加生成序列 $X_{\mathrm{AGO}}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{AGO}}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{AGO}}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{AGO}}^{{\mathrm{Xf}}_u}\left(n\right)\};$

所述的一次累加生成序列 $X_{\mathrm{AGO}}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{AGO}}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{AGO}}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{AGO}}^{{\mathrm{Xf}}_u}\left(n\right)\};$ 中

表示任意一个超标频段Xf_u下的第1个频率点上的累加生成值、表示任意一个超标频段Xf_u下的第2个频率点上的累加生成值、

表示任意一个超标频段Xf_u下的第n个频率点上的累加生成值，也称为任意一个频率点上的累加生成值，n表示任意一个超标频段Xf_u下的频率标识号。

步骤302：对所述的一次累加生成序列X_AGO(Xf_u)进行紧邻均值生成，得到紧邻均值生成序列 $Z_{\mathrm{avg}}\left({\mathrm{Xf}}_u\right)=\{z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(1\right),z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(2\right),...,z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(n\right)\};$

所述的紧邻均值生成序列 $Z_{\mathrm{avg}}\left({\mathrm{Xf}}_u\right)=\{z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(1\right),z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(2\right),...,z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(n\right)\}$ 中

表示任意一个超标频段Xf_u下的第1个频率点上的紧邻均值生成值、

表示任意一个超标频段Xf_u下的第2个频率点上的紧邻均值生成值、

表示任意一个超标频段Xf_u下的第n个频率点上的紧邻均值生成值，也称为在任意一个超标频段Xf_u下的任意一个频率点上的紧邻均值生成值，n表示任意一个超标频段Xf_u下的频率标识号。

步骤303：令

为GM(1，1)模型的原始形式，其中a表示发展系数的负数，b为灰色作用量；将a和b用参数列表示为AA＝[a，b]^T，T表示坐标转置。采取最小二乘法对参数列AA＝[a，b]^T进行求解得到AA＝(Z_add ^T(Xf_u)×Z_add(Xf_u))^-1×Z_add ^T(Xf_u)×X_extract(Xf_u)，所述X_extract(Xf_u)表示列向量，即 $X_{\mathrm{extract}}\left({\mathrm{Xf}}_u\right)=\left[\begin{array}{c}{x}_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(2\right)\\ x_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(3\right)\\ .\\ .\\ .\\ x_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(n\right)\end{array}\right]$ 为从频谱场强序列X_test(Xf_u)中提取的第2个到第n个频谱场强值组成的列向量，矩阵 $Z_{\mathrm{add}}\left({\mathrm{Xf}}_u\right)=\left[\begin{array}{cc}-z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(2\right)& 1\\ -z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(3\right)& 1\\ .& .\\ .& .\\ .& .\\ -z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(n\right)& 1\end{array}\right]$ 的第一列为从紧邻均值生成序列Z_avg(Xf_u)中提取的第2个到第n个元素取负组成的列向量，第二列为n-1个1组成的单位列向量，因此可求得灰色GM(1，1)模型的时间响应序列为e^-a(n-1)表示第n-1个频率上的指数分量，经过一次累减还原得到总体趋势序列 $X_{\mathrm{trend}}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{trend}}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{trend}}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{trend}}^{{\mathrm{Xf}}_u}\left(n\right)\},$ 其中

表示任意一个超标频段Xf_u下的第1个频率点上的总体趋势值、

表示任意一个超标频段Xf_u下的第2个频率点上的总体趋势值、

表示任意一个超标频段Xf_u下的第n个频率点上的总体趋势值，也称为任意一个频率点上的总体趋势值，n表示任意一个超标频段Xf_u下的频率标识号。

步骤四：提取超标频谱中的谐波分量

步骤401：计算频率上的场强值与总体趋势序列在频率上的模拟值之差，得到第一次残差序列 $X_{\mathrm{diff}1}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(n\right)\};$

所述第一次残差序列 $X_{\mathrm{diff}1}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(n\right)\}$ 中

表示任意一个超标频段Xf_u下的第1个频率点上的第一次残差值，

表示任意一个超标频段Xf_u下的第2个频率点上的第一次残差值，

表示任意一个超标频段Xf_u下的第n个频率点上的第一次残差值，也称为任意一个频率点上的第一次残差值；

步骤402：对第一次残差序列 $X_{\mathrm{diff}1}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(n\right)\}$ 采用方差分析法进行优势周期的提取，并将提取出来的优势周期记为m；然后采用均值生成公式

计算序列X_diff1(Xf_u)中每个元素的均值生成函数；对均值生成函数

按优势周期m作周期性延拓，得到周期延拓函数最后将GM(1，1)模型和周期延拓函数

叠加，得到灰色单周期模型序列 $X_{\mathrm{gp}1}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{gp}1}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{gp}1}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{gp}1}^{{\mathrm{Xf}}_u}\left(n\right)\};$

其中，n为频率标识号，l为小于n/m的最大整数，j为从0到l-1的整数，m为周期长度，t为均值生成函数的序号，其取值范围是[1，m]；周期延拓函数中k为延拓函数序列的序列号，取值范围是[1，n]，ta为对应的k对周期m取模后得到的值；

表示任意一个超标频段Xf_u下的第1个频率点上的灰色单周期模型拟合值；

表示任意一个超标频段Xf_u下的第2个频率点上的灰色单周期模型拟合值；

表示任意一个超标频段Xf_u下的第n个频率点上的灰色单周期模型拟合值，也称为任意一个频率点上的灰色单周期模型拟合值。

步骤403：计算曲线中的其他谐波分量，求灰色单周期模型序列X_gp1(Xf_u)与频谱场强序列X_test(Xf_u)的第二次残差序列 $X_{\mathrm{diff}2}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{diff}2}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{diff}2}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{diff}2}^{{\mathrm{Xf}}_u}\left(n\right)\};$

表示任意一个超标频段Xf_u下的第1个频率点上的第二次残差值；

表示任意一个超标频段Xf_u下的第2个频率点上的第二次残差值；

表示任意一个超标频段Xf_u下的第n个频率点上的第二次残差值，也称为任意一个频率点上的第二次残差值；

采用步骤402中对残差进行周期性建模的方法，对第二次残差序列X_diff2(Xf_u)进行周期建模，分别得到第二次残差序列X_diff2(Xf_u)的优势周期记为m₂、周期延拓函数

灰色双周期模型序列 $X_{\mathrm{gp}2}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{gp}2}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{gp}2}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{gp}2}^{{\mathrm{Xf}}_u}\left(n\right)\};$

步骤404：若m₂与m相同，则谐波分量的提取结束；

若m₂与m不相同，则继续进行第三次残差序列的周期性建模，直至与上一次提取的优势周期相同，则结束谐波分量的提取。

在本发明中，最后将具有总体趋势的GM(1，1)模型序列与所有具有周期趋势的谐波分量叠加，得到任意一个超标频段Xf_u的灰色多周期拟合模型序列 $X_{\mathrm{gps}}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{gps}}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{gps}}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{gps}}^{{\mathrm{Xf}}_u}\left(n\right)\},x_{\mathrm{gps}}^{{\mathrm{Xf}}_u}\left(1\right)$ 表示任意一个超标频段Xf_u下的第1个频率点上的灰色多周期模型拟合值，

表示任意一个超标频段Xf_u下的第2个频率点上的灰色多周期模型拟合值，

表示任意一个超标频段Xf_u下的第n个频率点上的灰色多周期模型拟合值，也称为任意一个频率点上的灰色多周期模型拟合值，其表达式为 $x_{\mathrm{gps}}^{{\mathrm{Xf}}_u}\left(n\right)=x_{\mathrm{trend}}^{{\mathrm{Xf}}_u}\left(n\right)+E_m^{{\mathrm{Xf}}_u}\left(n\right)+E_{m2}^{{\mathrm{Xf}}_u}\left(n\right)+...+E_{\mathrm{ms}}^{{\mathrm{Xf}}_u}\left(n\right);$ 其中为任意一个超标频段Xf_u下的具有周期值m的周期延拓函数，

为任意一个超标频段Xf_u下的具有周期值m₂的周期延拓函数，

为任意一个超标频段Xf_u下的具有周期值m_s的周期延拓函数。不同周期的谐波分量提取出来的顺序即定义为谐波序，它表示具有不同间隔周期的谐波按从大到小的顺序对整个干扰发射频谱曲线的贡献。

步骤五：谐波分量与器件频率的匹配确定干扰源

将通过步骤四分析得出的频谱曲线表达式 $x_{\mathrm{gps}}^{{\mathrm{Xf}}_u}\left(n\right)=x_{\mathrm{trend}}^{{\mathrm{Xf}}_u}\left(n\right)+E_m^{{\mathrm{Xf}}_u}\left(n\right)+E_{m2}^{{\mathrm{Xf}}_u}\left(n\right)+...+E_{\mathrm{ms}}^{{\mathrm{Xf}}_u}\left(n\right)$ 中的谐波分量的谐波的周期m、m₂、…m_s与干扰设备中所含的元器件的频率Freq进行比对；

(A)首先选取谐波周期m进行比对，若存在某个元器件的频率与其相等(即m＝Freq)，则认为该谐波分量

是由此元器件产生的干扰，确定该元器件为一个干扰源；

(B)若不存在元器件的频率与m相等(即m≠Freq)，则需进一步分析干扰设备中可能产生此频率间隔的元器件。

按(A)和(B)的对比法找出每个周期分量对应的的各个元器件，并按其对应的谐波提取顺序，得到具有不同间隔周期的谐波分量对整个干扰发射频谱曲线的贡献大小，按照所述谐波提取顺序依次采取干扰抑制措施，直至干扰设备的电磁干扰测试满足军标要求。

在本发明中，重复步骤三、步骤四和步骤五遍历超标频段FX＝{Xf₁，Xf₂，…，Xf_u}中所有频段，得到所有超标频段的灰色多周期拟合模型序列FX_gps＝{X_gps(Xf₁)，X_gps(Xf₂)，…，X_gps(Xf_u)}，X_gps(Xf₁)表示第1个超标频段Xf₁的灰色多周期拟合模型序列，X_gps(Xf₁)表示第2个超标频段Xf₂的灰色多周期拟合模型序列，X_gps(Xf_u)表示第u个超标频段Xf_u的灰色多周期拟合模型序列，也称任意一个超标频段Xf_u的灰色多周期拟合模型序列。

本发明提出的一种采用灰色多周期模型提取电磁频谱中谐波分量的方法，首先采集干扰设备的频谱数据；然后利用灰色系统模型拟合出频谱数据的总体趋势分量；然后利用周期模型对频谱数据去除总体趋势分量后的数据提取谐波分量，获得谐波的量化表达式；最后将谐波量化表达式中的谐波频率与干扰设备中电路元器件的工作频率进行比对，若与某个元器件的工作频率相等则可确定对应谐波分量产生的根源。本发明方法将频谱数据中谐波干扰的特性用数学函数式表达出来，物理特征明确，便于电磁干扰频谱数据部分特性的认知，起到指导干扰设备电磁兼容性分析和设计的作用。

实施例

下面通过一个具体的实施例来对本发明提供的谐波分析方法做出进一步的阐释。具体处理步骤如下：

步骤一：将电磁干扰检测试验中得到的某干扰设备的频谱发射曲线的频率、场强值读入Matlab 7.5.0，通过图形绘制命令将其显示出来如图2所示。将曲线按照斜率一致性进行分段，得到六个频率分段序号，如下表所示。

序号

f1

f2

f3

f4

f5

f6

频段起始点(Hz)

0.3×108

0.55×108

0.8×108

1.1×108

1.5×108

1.8×108

频段终止点(Hz)

0.55×108

0.8×108

1.1×108

1.5×108

1.8×108

2×108

步骤二：根据与军标极限比对从f1～f6提取需要分析的超标频段。从图2可以看出f4频段曲线超标严重，选取其作为研究对象。取频率间隔为47kHz，得到f4频段的1000个频率序列作为分析的频谱场强序列，将1～1000作为该频段上频率的序列号。对频谱场强数据序列进行一次累加生成，得到曲线如图3。图中横坐标为频率序列号，纵坐标为场强值(单位dBμV/m)。

步骤三：由图3曲线可以看出斜率没有发生突变，总体趋势明确，可以直接对该段频谱数据进行GM(1，1)建模。频谱场强数据的时间响应序列为：

用最小二乘法计算估计参数列得到： $\mathrm{AA}=\left[\begin{array}{c}a\\ b\end{array}\right]=\left[\begin{array}{c}-0.0001\\ 34.9886\end{array}\right],$ 累减得到还原值序列的表达式为： $x_{\mathrm{trend}}\left(n\right)={\hat{x}}_{\mathrm{GM}}\left(n\right)-{\hat{x}}_{\mathrm{GM}}(n-1)$

$=(1-e^a)\×(x_{\mathrm{test}}\left(1\right)-\frac{b}{a})\×e^{-a(n-1)}.$

$=34.9907\×e^{0.0001(n-1)}$

将测试得到的频谱场强值序列与GM(1，1)模型拟合序列示于图4，图中横坐标为频率序列号，纵坐标为场强值(单位dBμV/m)，可以看到GM(1，1)模型描述出频谱场强数据总体趋势。

步骤四：提取超标频段f4中的谐波分量。首先计算残差数列，并对残差序列进行周期性建模。通过方差分析计算得到的优势周期为117，于是可得灰色单周期模型的表达式为：x_gp1(n)＝x_trend(n)+E₁₁₇(n)。灰色单周期模型曲线与频谱场强序列的对比如图5，图中横坐标为频率序列号，纵坐标为场强值(单位dBμV/m)。从图中可以看出有周期趋势的分量已经很明显，除了个别频点的幅度差异，总体令人满意，且误差超过10％的数据量只有5％。

计算曲线中的其他谐波分量。计算出二次残差序列，对二次残差序列再一次进行周期建模，提取的优势周期仍为117，一方面说明第一次周期模型的幅值拟合存在一些偏差，二次周期提取可以减弱幅度误差；另一方面说明序列已经不存在其他的周期分量，周期拟合过程可以结束，即原序列只存在一个谐波频率。图6为经过二次周期拟合的序列与原序列的对比图形，可以看出拟合比较精确，二次周期建模的平均误差，下降为3％；图中横坐标为频率序列号，纵坐标为场强值(单位dBμV/m)。由此得出初步结论：干扰设备在f4频段的超标是由周期间隔为117的谐波分量引起的。

步骤五：将分析得出的谐波频率m与干扰设备中的元器件的频率Freq比对。只有一个谐波周期为117，而每个数据间的频谱间隔为47kHz，因此可知是某器件产生了频率间隔为5.5MHz的谐波干扰。用过比对发现电路中存在66MHz的晶振时钟，每12个周期产生一个CPU机器周期，所以将电磁干扰源定位到CPU板上。该发明提供了一种快速诊断的手段，为后续的电磁兼容性整改给出依据。

Claims (3)

1.一种采用灰色多周期模型提取电磁频谱中谐波分量的方法，其特征在于：所述电磁频谱中谐波分量的提取有下列步骤：

步骤一：对频谱进行分段

采用斜率一致性原则对频率AF为横坐标，场强值AM为纵坐标绘制出的频谱幅度图中的频谱曲线进行分段，得到互不交叠的频谱频段F＝{f₁，f₂，…，f_w}；

所述的频谱频段F＝{f₁，f₂，…，f_w}中w表示频段标识号，f₁表示第1个频段、f₂表示第2个频段、f_w，表示第w个频段，也称为任意一个频段；

步骤二：选取超标频段

将所述的频谱频段F＝{f₁，f₂，…，f_w}与国军标GJB151A中规定的发射极限进行对比，提取出超标频段FX＝{Xf₁，Xf₂，…，Xf_u}；

所述的超标频段FX＝{Xf₁，Xf₂，…，Xf_u}中u表示超标频段标识号，Xf₁表示第1个超标频段、Xf₂表示第2个超标频段、Xf_u表示第u个超标频段，也称为任意一个超标频段；

任意一个超标频段Xf_u中包括有多个频率、以及每个频率对应的场强值，将所述的任意一个超标频段Xf_u用频谱场强序列表示为 $X_{\mathrm{test}}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(n\right)\};$

所述的频谱场强序列 $X_{\mathrm{test}}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(n\right)\}$ 中表示在任意一个超标频段Xf_u下的第1个频率上的场强值、

表示任意一个超标频段Xf_u下的第2个频率上的场强值、

表示任意一个超标频段Xf_u下的第n个频率上的场强值，n表示任意一个超标频段Xf_u下的频率标识号；

步骤三：构建总体趋势序列

从任意一个超标频段Xf_u中提取出总体趋势序列X_trend(Xf_u)的步骤如下：

步骤301：对测试得到的任意一个超标频段Xf_u上的频谱场强序列

进行一次累加，生成一次累加生成序列 $X_{\mathrm{AGO}}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{AGO}}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{AGO}}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{AGO}}^{{\mathrm{Xf}}_u}\left(n\right)\};$

所述的一次累加生成序列 $X_{\mathrm{AGO}}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{AGO}}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{AGO}}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{AGO}}^{{\mathrm{Xf}}_u}\left(n\right)\}$ 中表示任意一个超标频段Xf_u下的第1个频率点上的累加生成值、表示任意一个超标频段Xf_u下的第2个频率点上的累加生成值、

表示任意一个超标频段Xf_u下的第n个频率点上的累加生成值，也称为任意一个频率点上的累加生成值，n表示任意一个超标频段Xf_u下的频率标识号；

步骤302：对所述的一次累加生成序列X_AGO(Xf_u)进行紧邻均值生成，得到紧邻均值生成序列 $Z_{\mathrm{avg}}\left({\mathrm{Xf}}_u\right)=\{z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(1\right),z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(2\right),...,z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(n\right)\};$

所述的紧邻均值生成序列 $Z_{\mathrm{avg}}\left({\mathrm{Xf}}_u\right)=\{z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(1\right),z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(2\right),...,z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(n\right)\};$ 中

表示任意一个超标频段Xf_u下的第1个频率点上的紧邻均值生成值、

表示任意一个超标频段Xf_u下的第2个频率点上的紧邻均值生成值、

表示任意一个超标频段Xf_u下的第n个频率点上的紧邻均值生成值，也称为在任意一个超标频段Xf_u下的任意一个频率点上的紧邻均值生成值，n表示任意一个超标频段Xf_u下的频率标识号；

步骤303：令

为GM(1，1)模型的原始形式，其中a表示发展系数的负数，b为灰色作用量；将a和b用参数列表示为AA＝[a，b]^T，T表示坐标转置；采取最小二乘法对参数列AA＝[a，b]^T进行求解得到AA＝(Z_add ^T(Xf_u)×Z_add(Xf_u))^-1×Z_add ^T(Xf_u)×X_extract(Xf_u)，所述X_extract(Xf_u)表示列向量，即 $X_{\mathrm{extract}}\left({\mathrm{Xf}}_u\right)=\left[\begin{array}{c}{x}_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(2\right)\\ x_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(3\right)\\ .\\ .\\ .\\ x_{\mathrm{test}}^{{\mathrm{Xf}}_u}\left(n\right)\end{array}\right]$ 为从频谱场强序列X_testt(Xf_u)中提取的第2个到第n个频谱场强值组成的列向量，矩阵 $Z_{\mathrm{add}}\left({\mathrm{Xf}}_u\right)=\left[\begin{array}{cc}-z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(2\right)& 1\\ -z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(3\right)& 1\\ .& .\\ .& .\\ .& .\\ -z_{\mathrm{avg}}^{{\mathrm{Xf}}_u}\left(n\right)& 1\end{array}\right]$ 的第一列为从紧邻均值生成序列Z_avg(Xf_u)中提取的第2个到第n个元素取负组成的列向量，第二列为n-1个1组成的单位列向量，因此可求得灰色GM(1，1)模型的时间响应序列为

e^-a(n-1)表示第n-1个频率上的指数分量，经过一次累减还原得到总体趋势序列 $X_{\mathrm{trend}}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{trend}}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{trend}}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{trend}}^{{\mathrm{Xf}}_u}\left(n\right)\},$ 其中

表示任意一个超标频段Xf_u下的第1个频率点上的总体趋势值、

表示任意一个超标频段Xf_u下的第2个频率点上的总体趋势值、

表示任意一个超标频段Xf_u下的第n个频率点上的总体趋势值，也称为任意一个频率点上的总体趋势值，n表示任意一个超标频段Xf_u下的频率标识号；

步骤四：提取超标频谱中的谐波分量

步骤401：计算频率上的场强值与总体趋势序列在频率上的模拟值之差，得到第一次残差序列 $X_{\mathrm{diff}1}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(n\right)\};$

所述第一次残差序列 $X_{\mathrm{diff}1}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(n\right)\}$ 中

表示任意一个超标频段Xf_u下的第1个频率点上的第一次残差值，表示任意一个超标频段Xf_u下的第2个频率点上的第一次残差值，

表示任意一个超标频段Xf_u下的第n个频率点上的第一次残差值，也称为任意一个频率点上的第一次残差值；

步骤402：对第一次残差序列 $X_{\mathrm{diff}1}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{diff}1}^{{\mathrm{Xf}}_u}\left(n\right)\}$ 采用方差分析法进行优势周期的提取，并将提取出来的优势周期记为m；然后采用均值生成公式

计算序列X_diff1(Xf_u)中每个元素的均值生成函数；对均值生成函数

按优势周期m作周期性延拓，得到周期延拓函数

最后将GM(1，1)模型和周期延拓函数

叠加，得到灰色单周期模型序列 $X_{\mathrm{gp}1}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{gp}1}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{gp}1}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{gp}1}^{{\mathrm{Xf}}_u}\left(n\right)\};$

其中，n为频率标识号，l为小于n/m的最大整数，j为从0到l-1的整数，m为周期长度，t为均值生成函数的序号，其取值范围是[1，m]；周期延拓函数

中k为延拓函数序列的序列号，取值范围是[1，n]，ta为对应的k对周期m取模后得到的值；

表示任意一个超标频段Xf_u下的第1个频率点上的灰色单周期模型拟合值；

表示任意一个超标频段Xf_u下的第2个频率点上的灰色单周期模型拟合值；

表示任意一个超标频段Xf_u下的第n个频率点上的灰色单周期模型拟合值，也称为任意一个频率点上的灰色单周期模型拟合值；

步骤403：计算曲线中的其他谐波分量，求灰色单周期模型序列X_gp1(Xf_u)与频谱场强序列X_test(Xf_u)的第二次残差序列 $X_{\mathrm{diff}2}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{diff}2}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{diff}2}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{diff}2}^{{\mathrm{Xf}}_u}\left(n\right)\};$

表示任意一个超标频段Xf_u下的第1个频率点上的第二次残差值；

表示任意一个超标频段Xf_u下的第2个频率点上的第二次残差值；

表示任意一个超标频段Xf_u下的第n个频率点上的第二次残差值，也称为任意一个频率点上的第二次残差值；

采用步骤402中对残差进行周期性建模的方法，对第二次残差序列X_diff2(Xf_u)进行周期建模，分别得到第二次残差序列X_diff2(Xf_u)的优势周期记为m₂、周期延拓函数

灰色双周期模型序列 $X_{\mathrm{gp}2}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{gp}2}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{gp}2}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{gp}2}^{{\mathrm{Xf}}_u}\left(n\right)\};$

步骤404：若m₂与m相同，则谐波分量的提取结束；

若m₂与m不相同，则继续进行第三次残差序列的周期性建模，直至与上一次提取的优势周期相同，则结束谐波分量的提取；

最后将具有总体趋势的GM(1，1)模型序列与所有具有周期趋势的谐波分量叠加，得到任意一个超标频段Xf_u的灰色多周期拟合模型序列 $X_{\mathrm{gps}}\left({\mathrm{Xf}}_u\right)=\{x_{\mathrm{gps}}^{{\mathrm{Xf}}_u}\left(1\right),x_{\mathrm{gps}}^{{\mathrm{Xf}}_u}\left(2\right),...,x_{\mathrm{gps}}^{{\mathrm{Xf}}_u}\left(n\right)\},x_{\mathrm{gps}}^{{\mathrm{Xf}}_u}\left(1\right)$ 表示任意一个超标频段Xf_u下的第1个频率点上的灰色多周期模型拟合值，

表示任意一个超标频段Xf_u下的第2个频率点上的灰色多周期模型拟合值，

表示任意一个超标频段Xf_u下的第n个频率点上的灰色多周期模型拟合值，也称为任意一个频率点上的灰色多周期模型拟合值，其表达式为 $x_{\mathrm{gps}}^{{\mathrm{Xf}}_u}\left(n\right)=x_{\mathrm{trend}}^{{\mathrm{Xf}}_u}\left(n\right)+E_m^{{\mathrm{Xf}}_u}\left(n\right)+E_{m2}^{{\mathrm{Xf}}_u}\left(n\right)+...+E_{\mathrm{ms}}^{{\mathrm{Xf}}_u}\left(n\right)$ 其中

为任意一个超标频段Xf_u下的具有周期值m的周期延拓函数，

为任意一个超标频段Xf_u下的具有周期值m₂的周期延拓函数，为任意一个超标频段Xf_u下的具有周期值m_s的周期延拓函数；不同周期的谐波分量提取出来的顺序即定义为谐波序，它表示具有不同间隔周期的谐波按从大到小的顺序对整个干扰发射频谱曲线的贡献；

步骤五：谐波分量与器件频率的匹配确定干扰源

将通过步骤四分析得出的频谱曲线表达式 $x_{\mathrm{gps}}^{{\mathrm{Xf}}_u}\left(n\right)=x_{\mathrm{trend}}^{{\mathrm{Xf}}_u}\left(n\right)+E_m^{{\mathrm{Xf}}_u}\left(n\right)+E_{m2}^{{\mathrm{Xf}}_u}\left(n\right)+...+E_{\mathrm{ms}}^{{\mathrm{Xf}}_u}\left(n\right)$ 中的谐波分量

的谐波的周期m、m₂、…m_s与干扰设备中所含的元器件的频率Freq进行比对；

(A)首先选取谐波周期m进行比对，若存在某个元器件的频率与其相等(即m＝Freq)，则认为该谐波分量

是由此元器件产生的干扰，确定该元器件为一个干扰源；

(B)若不存在元器件的频率与m相等(即m≠Freq)，则需进一步分析干扰设备中可能产生此频率间隔的元器件；

按(A)和(B)的对比法找出每个周期分量对应的的各个元器件，并按其对应的谐波提取顺序，得到具有不同间隔周期的谐波分量对整个干扰发射频谱曲线的贡献大小，按照所述谐波提取顺序依次采取干扰抑制措施，直至干扰设备的电磁干扰测试满足军标要求；

重复步骤三、步骤四和步骤五遍历超标频段FX＝{Xf₁，Xf₂，…，Xf_u}中所有频段，得到所有超标频段的灰色多周期拟合模型序列FX_gps＝{X_gps(Xf₁)，X_gps(Xf₂)，…，X_gps(Xf_u)}，X_gps(Xf₁)表示第1个超标频段Xf₁的灰色多周期拟合模型序列，X_gps(Xf₁)表示第2个超标频段Xf₂的灰色多周期拟合模型序列，X_gps(Xf_u)表示第u个超标频段Xf_u的灰色多周期拟合模型序列，也称任意一个超标频段Xf_u的灰色多周期拟合模型序列。

2.根据权利要求1所述的采用灰色多周期模型提取电磁频谱中谐波分量的方法，其特征在于：利用电磁兼容性检测软件EMC2003获取电磁频谱数据。

3.根据权利要求1所述的采用灰色多周期模型提取电磁频谱中谐波分量的方法，其特征在于：利用Matlab软件进行仿真建模。

Images