CN102967761B - 频率自适应信号分析方法 - Google Patents

Abstract

本发明公开了一种频率自适应信号分析方法，以T秒为采样周期对被测信号进行定时采样和量化得到采样数据；设定交流分量个数N的值，通过迭代直接获得直流分量x₀，通过离散频率估计算法直接获得估计频率ω₁、ω₂、…、ω_N，交流分量x₁、x₂、…、x_N和正交分量z₁、z₂、…、z_N，并计算估计幅值a₁、a₂、…、a_N。本发明的频率自适应信号分析方法由于采用迭代方法进行信号分析，易于在计算机系统中应用，并且具有较高的精度和较快的收敛速度。

Description

频率自适应信号分析方法

技术领域

本发明属于信号分析技术领域，具体涉及一种用于把信号分解为直流分量与多个交流分量，并获得每个交流分量的正交分量、频率和幅值的信号分析方法。

背景技术

在工程中，经常需要把电压、电流、声音、振动等信号分解为直流分量与多个交流分量，计算直流分量的数值以及各个交流分量的频率和幅值，例如电力系统的谐波与间谐波分析以及机械振动的振动分量分析等。

为此，现有技术中已经提出了几种基于频率估计的信号分析方法，这些方法虽然在一定程度上能够获得较好的分解结果，但都还存在一定的不足，例如：

基于微分方程形式的针对连续时间信号的信号分析方法，不能直接应用于计算机系统中；

采用一阶无限冲激响应（IIR）算法的信号分析方法，虽然具有较高的精度和较快的收敛速度，但是由于每次迭代的运行时间可能不相等，因而难以确定定时采样周期的数值，软件编程实现复杂；

采用一阶有限冲激响应（FIR）算法的方法虽然软件编程实现简单，但是仅具有一阶的收敛速度与精度。

发明内容

本发明的目的是提供一种频率自适应信号分析方法，以解决现有技术中的信号分析方法无法同时满足易于在计算机系统中应用、具有较高的精度和较快的收敛速度、并且易于软件编程实现的问题。

为了实现上述目的，本发明提供了一种频率自适应信号分析方法，以T秒为采样周期对被测信号进行定时采样和量化得到采样数据；设定交流分量个数N的值，通过迭代直接获得直流分量x₀，通过离散频率估计算法直接获得估计频率ω₁、ω₂、…、ω_N，交流分量x₁、x₂、…、x_N和正交分量z₁、z₂、…、z_N，并计算估计幅值a₁、a₂、…、a_N。

作为优选，针对被测信号的采样数据u[s]，按下述第一处理方案和第二处理方案中的一种执行，连续的两个采样数据分别执行不同的处理方案；其中，β₁、β₂、…、β_2N-1、β_2N为依次递增且均不大于2π/T的正数，0＜μ<2π/T，0<γ<10⁶，分母函数b为大于零的二元函数；

第一处理方案依次包括以下步骤：

S101：利用式(1)获得直流分量增量h₀[2]；

$\left\{\begin{array}{c}e\left[2\right]=u\left[s\right]-\underset{m=0}{\overset{N}{\mathrm{\&Sigma;}}}(x_m+T\&CenterDot;h_m[1\left]\right)\\ h_0\left[2\right]=\mathrm{\&mu;}\&CenterDot;e\left[2\right]\end{array}\right.---\left(1\right)$

S102：令下标n依次取值为1,2,…,N，循环执行式(2)，获得交流分量增量h_n[2]、正交分量增量g_n[2]和估计频率增量w_n[2]；

$\left\{\begin{array}{c}{h}_n\left[2\right]=\mathrm{\&mu;}\&CenterDot;e\left[2\right]+({\mathrm{\&omega;}}_n+T\&CenterDot;w_n[1\left]\right)(z_n+T\&CenterDot;g_n[1\left]\right)\\ g_n\left[2\right]=-({\mathrm{\&omega;}}_n+T\&CenterDot;w_n[1\left]\right)(x_n+T\&CenterDot;h_n[1\left]\right)\\ w_n\left[2\right]=\mathrm{\&gamma;}\&CenterDot;\mathrm{\&mu;}\&CenterDot;e\left[2\right]\&CenterDot;(z_n+T\&CenterDot;g_n[1\left]\right)/b((x_n+T\&CenterDot;h_n[1\left]\right),(z_n+T\&CenterDot;g_n[1\left]\right))\end{array}\right.---\left(2\right)$

S103：利用式(3)获得直流分量增量h₀[3]；

$\left\{\begin{array}{c}e\left[3\right]=u\left[s\right]-\underset{m=0}{\overset{N}{\mathrm{\&Sigma;}}}(x_m+T\&CenterDot;h_m[2\left]\right)\\ h_0\left[3\right]=\mathrm{\&mu;}\&CenterDot;e\left[3\right]\end{array}\right.---\left(3\right)$

S104：令下标n依次取值为1,2,…,N，循环执行式(4)，获得交流分量增量h_n[3]、正交分量增量g_n[3]和估计频率增量w_n[3]；

$\left\{\begin{array}{c}{h}_n\left[3\right]=\mathrm{\&mu;}\&CenterDot;e\left[3\right]+({\mathrm{\&omega;}}_n+T\&CenterDot;w_n[2\left]\right)(z_n+T\&CenterDot;g_n[2\left]\right)\\ g_n\left[3\right]=-({\mathrm{\&omega;}}_n+T\&CenterDot;w_n[2\left]\right)(x_n+T\&CenterDot;h_n[2\left]\right)\\ w_n\left[3\right]=\mathrm{\&gamma;}\&CenterDot;\mathrm{\&mu;}\&CenterDot;e\left[3\right]\&CenterDot;(z_n+T\&CenterDot;g_n[2\left]\right)/b((x_n+T\&CenterDot;h_n[2\left]\right),(z_n+T\&CenterDot;g_n[2\left]\right))\end{array}\right.---\left(4\right)$

第二处理方案依次包括以下步骤：

S201：先利用式(5)获得直流分量增量h₀[4]，再依据式(6)对直流分量x₀进行迭代处理；

$\left\{\begin{array}{c}e\left[4\right]=u\left[s\right]-\underset{m=0}{\overset{N}{\mathrm{\&Sigma;}}}(x_m+2T\&CenterDot;h_m[3\left]\right)\\ h_0\left[4\right]=\mathrm{\&mu;}\&CenterDot;e\left[4\right]\end{array}\right.---\left(5\right)$

$x_0\&LeftArrow;x_0+\frac{T}{3}\left(h_0\right[4]+2h_0[3]+2h_0[2]+h_0[1\left]\right)---\left(6\right)$

S202：令下标n依次取值为1,2,…,N，循环执行式(7)和式(8)；其中，每次循环，先利用式(7)获得交流分量增量h_n[4]、正交分量增量g_n[4]和估计频率增量w_n[4]，再利用式(8)分别对交流分量x_n、正交分量z_n、估计频率ω_n进行迭代处理，并对迭代后的估计频率ω_n进行限幅处理以满足条件β_2n-1≤ω_n＜β_2n；依据迭代后的交流分量x_n和正交分量z_n计算估计幅值a_n的值；

$\left\{\begin{array}{c}{h}_n\left[4\right]=\mathrm{\&mu;}\&CenterDot;e\left[4\right]+({\mathrm{\&omega;}}_n+2T\&CenterDot;w_n[3\left]\right)(z_n+2T\&CenterDot;g_n[3\left]\right)\\ g_n\left[4\right]=-({\mathrm{\&omega;}}_n+2T\&CenterDot;w_n[3\left]\right)(x_n+2T\&CenterDot;h_n[3\left]\right)\\ w_n\left[4\right]=\mathrm{\&gamma;}\&CenterDot;\mathrm{\&mu;}\&CenterDot;e\left[4\right]\&CenterDot;(z_n+2T\&CenterDot;g_n[3\left]\right)/b((x_n+2T\&CenterDot;h_n[3\left]\right),(z_n+2T\&CenterDot;g_n[3\left]\right))\end{array}\right.---\left(7\right)$

$\left\{\begin{array}{c}{x}_n\&LeftArrow;x_n+\frac{T}{3}\left(h_n\right[4]+2h_n[3]+2h_n[2]+h_n[1\left]\right)\\ z_n\&LeftArrow;z_n+\frac{T}{3}\left(g_n\right[4]+2g_n[3]+2g_n[2]+g_n[1\left]\right)\\ {\mathrm{\&omega;}}_n\&LeftArrow;{\mathrm{\&omega;}}_n+\frac{T}{3}\left(w_n\right[4]+2w_n[3]+2w_n[2]+w_n[1\left]\right)\\ a_n=\sqrt{{x_n}^2+{z_n}^2}\end{array}\right.---\left(8\right)$

S203：依据迭代处理后的直流分量x₀和交流分量x₁、x₂、…、x_N的当前值，利用式(9)获得直流分量增量h₀[1]；

$\left\{\begin{array}{c}e\left[1\right]=u\left[s\right]-\underset{m=0}{\overset{N}{\mathrm{\&Sigma;}}}{x}_m\\ h_0\left[1\right]=\mathrm{\&mu;}\&CenterDot;e\left[1\right]\end{array}\right.---\left(9\right)$

S204：令下标n依次取值为1,2,…,N，循环执行式(10)，

$\left\{\begin{array}{c}{h}_n\left[1\right]=\mathrm{\&mu;}\&CenterDot;e\left[1\right]+{\mathrm{\&omega;}}_n\&CenterDot;z_n\\ g_n\left[1\right]=-{\mathrm{\&omega;}}_n\&CenterDot;x_n\\ w_n\left[1\right]=\mathrm{\&gamma;}\&CenterDot;\mathrm{\&mu;}\&CenterDot;e\left[1\right]\&CenterDot;z_n/b(x_n,z_n)\end{array}\right.---\left(10\right)$

依据迭代处理后的交流分量x_n和正交分量z_n，获得交流分量增量h_n[1]、正交分量增量g_n[1]和估计频率增量w_n[1]。

作为进一步地优选，当采用非归一化方法时，所述分母函数b为大于0的常数。

作为进一步地优选，当采用非归一化方法时，所述分母函数b=1。

作为进一步地优选，当采用归一化方法时，所述分母函数b的函数式为b(x,y)=ε+x²+y²或者b(x,y)=max(ε,x²+y²)，其中正数ε不大于量化误差，函数max表示返回两个自变量中的最大值。

作为进一步地优选，具体包括以下步骤：

S1：设定参数T、N、β₁、β₂、…、β_2N-1、β_2N、μ、γ、ε的值，设定直流分量x₀，交流分量x₁、x₂、…、x_N，正交分量z₁、z₂、…、z_N和估计频率ω₁、ω₂、…、ω_N的初值；

S2：设定直流分量增量h₀[1]，交流分量增量h₁[1]、h₂[1]、…、h_N[1]，正交分量增量g₁[1]、g₂[1]、…、g_N[1]以及估计频率增量w₁[1]、w₂[1]、…、w_N[1]的初值，设定标志字的初值为执行第一处理方案；

S3：对被测信号的采样数据u[s]，依据标志字的值执行第一处理方案或者第二处理方案，然后修改标志字的值为执行另外一种处理方案；

S4：返回步骤S3循环执行，直至迭代完成。

作为进一步地优选，将步骤S2替换为：设定直流分量增量h₀[1]、h₀[2]、h₀[3]，交流分量增量h₁[1]、h₁[2]、h₁[3]、h₂[1]、h₂[2]、h₂[3]、…、h_N[1]、h_N[2]、h_N[3]，正交分量增量g₁[1]、g₁[2]、g₁[3]、g₂[1]、g₂[2]、g₂[3]、…、g_N[1]、g_N[2]、g_N[3]以及估计频率增量w₁[1]、w₁[2]、w₁[3]、w₂[1]、w₂[2]、w₂[3]、…、w_N[1]、w_N[2]、w_N[3]的初值，设定标志字的初值为执行第二处理方案。

与现有技术相比，本发明具有以下有益效果：

1、本发明针对经定时采样后得到的离散时间数据，采用迭代方法进行信号分析，便于利用计算机实现；

2、本发明以迭代方法直接得到直流分量，消除了信号的直流分量对计算交流分量及其幅值和频率的不良影响，拓宽了应用范围；

3、本发明以迭代方法直接得到交流分量和正交分量，不需要计算正弦函数和余弦函数，结构简单、运算量小；

4、本发明直接估计交流分量的频率，适于分析交流分量频率未知的概周期信号，不要求被测信号各交流分量的频率之间保持特定的数值关系，便于分析频率未知的谐波和间谐波成分；

5、本发明属于四阶有限脉冲响应算法，具有四阶精度和四阶收敛速度，较无限冲击响应算法更加易于实现，比一阶有限脉冲响应算法具有更高的精度和更快的收敛速度；

6、本发明的归一化频率自适应信号分析方法，在计算估计频率增量时进行了除法运算，分母与估计幅值的平方有关，抵消了包含在分子中的幅值的平方，使得估计频率的收敛速度基本不受交流分量幅值大小的影响，增强了鲁棒性。

附图说明

图1为本发明实施例一的频率自适应信号分析方法的流程示意图；

图2为本发明实施例二的频率自适应信号分析方法的流程示意图；

图3为归一化方法在计算机中仿真运行后所得信号分量跟踪误差的变化曲线；

图4为非归一化方法在计算机中仿真运行后所得信号分量跟踪误差的变化曲线；

图5为在参数γ和参数μ分别取不同数值的情况下归一化方法仿真运行所得的交流幅值跟踪误差变化曲线；

图6为在参数γ和参数μ分别取不同数值的情况下非归一化方法仿真运行所得的交流幅值跟踪误差变化曲线；

图7为在参数γ和参数μ分别取不同数值的情况下归一化方法仿真运行所得的交流频率跟踪误差变化曲线；

图8为在参数γ和参数μ分别取不同数值的情况下非归一化方法仿真运行所得的交流频率跟踪误差变化曲线。

具体实施方式

下面结合附图对本发明的具体实施例进行详细说明。

本发明的频率自适应信号分析方法包括离线分析与在线分析两种实施方式，其应用包括分析电压和电流的谐波与间谐波成分，以及分析不同频段的机械振动分量。

实施例一：

图1为本发明实施例一的频率自适应信号分析方法的流程示意图，实施例一的频率自适应信号分析方法为离线分析实施方式，其适用于分析保存在存储器件中的经周期采样得到的被测信号的离散时间序列。

离线分析实施方式的特点是采样与分析计算分开进行。首先以T秒为采样周期，定时采样被测信号，把所得采样数据存储在存储器件中，形成有K个数据的离散时间序列，然后再对离散时间序列进行分析计算。

如图1所示，在实施例一的频率自适应信号分析方法中，分析计算过程包括初始化步骤和数据处理步骤，其中数据处理步骤包括增量计算及数据迭代步骤、循环控制步骤两部分。

在初始化步骤中，首先设定离散时间序列中数据的个数K，设定当前处理的数据在离散时间序列中的位置号k为1，设定当前迭代处理后的数据在输出序列中的位置号m为1；接着设定参数N、β₁、β₂、…、β_2N-1、β_2N、μ、γ、ε的值，设定直流分量x₀，交流分量x₁、x₂、…、x_N，正交分量z₁、z₂、…、z_N和估计频率ω₁、ω₂、…、ω_N的初值；再设定直流分量增量h₀[1]，交流分量增量h₁[1]、h₂[1]、…、h_N[1]，正交分量增量g₁[1]、g₂[1]、…、g_N[1]以及估计频率增量w₁[1]、w₂[1]、…、w_N[1]的初值，设定标志字FLAG的初值为执行第一处理方案（将在下文中进行详细说明）的数值1；然后进入数据处理步骤循环执行。

在增量计算及数据迭代步骤中，先读取离散时间序列中的第k个数据作为当前处理的采样数据u[s]；接着依据标志字FLAG的值执行第一处理方案或者执行第二处理方案（将在下文中进行详细说明），在执行第二处理方案时，为查看数据变化趋势，把迭代处理后的数据存储到输出序列中，并把输出序列中的位置号m增加1；再对标志字FLAG进行逻辑取反运算，确定下次循环执行另外一种处理方案；然后进入循环控制步骤。

在循环控制步骤中，先把位置号k增加1，再依据位置号k和离散时间序列中数据的个数K的值判断是否返回执行增量计算及数据迭代步骤。若k≤K，返回执行增量计算及数据迭代步骤；若k＞K，表示离散时间序列中所有采样数据都已处理完毕，应该终止运行，结束离线分析过程。

实施例二：

图2为本发明实施例二的频率自适应信号分析方法的流程示意图，实施例二的频率自适应信号分析方法为在线分析实施方式，其适用于对被测信号的每个采样数据进行实时分析处理，在线分析实施方式的特点是一边采样一边分析计算，即每采样一个数据，就进行一次分析计算。

如图2所示，实施例二的频率自适应信号分析方法包括主步骤和定时中断步骤。

主步骤包括初始化步骤和实现其它的功能的主循环步骤。在初始化步骤中，首先设定参数N、β₁、β₂、…、β_2N-1、β_2N、μ、γ、ε的值，设定直流分量x₀，交流分量x₁、x₂、…、x_N，正交分量z₁、z₂、…、z_N和估计频率ω₁、ω₂、…、ω_N的初值；接着设定直流分量增量h₀[1]，交流分量增量h₁[1]、h₂[1]、…、h_N[1]，正交分量增量g₁[1]、g₂[1]、…、g_N[1]以及估计频率增量w₁[1]、w₂[1]、…、w_N[1]的初值，设定标志字FLAG的初值为执行第一处理方案（将在下文中进行详细说明）的数值1；再设定定时器的定时时间为满足香农采样定理要求的T秒（等于采样周期），并开放系统的定时中断。

在定时中断步骤中，首先保存中断现场的各个寄存器的当前值；接着对被测信号进行采样获得采样数据u[s]，依据标志字FLAG的值决定执行第一处理方案计算增量数值或者执行第二处理方案计算增量数值并进行迭代处理，然后对标志字FLAG进行逻辑取反运算，使下次中断执行另外一种处理方案；再恢复中断现场的各个寄存器的值；然后中断返回到主步骤执行。

在线分析实施方式通过定时器循环产生定时中断事件，引起定时中断步骤循环执行。

实施例一和实施例二的频率自适应信号分析方法中所述的第一处理方案和第二处理方案分别如下。不论是离线分析还是在线分析，针对被测信号的采样数据u[s]，均按下述第一处理方案和第二处理方案中的一种执行。

第一处理方案依次包括：

S101：利用式(1)获得直流分量增量h₀[2]；

$\left\{\begin{array}{c}e\left[2\right]=u\left[s\right]-\underset{m=0}{\overset{N}{\mathrm{\&Sigma;}}}(x_m+T\&CenterDot;h_m[1\left]\right)\\ h_0\left[2\right]=\mathrm{\&mu;}\&CenterDot;e\left[2\right]\end{array}\right.---\left(1\right)$

S102：令下标n依次取值为1,2,…,N，循环执行式(2)，获得交流分量增量h_n[2]、正交分量增量g_n[2]和估计频率增量w_n[2]；

$\left\{\begin{array}{c}{h}_n\left[2\right]=\mathrm{\&mu;}\&CenterDot;e\left[2\right]+({\mathrm{\&omega;}}_n+T\&CenterDot;w_n[1\left]\right)(z_n+T\&CenterDot;g_n[1\left]\right)\\ g_n\left[2\right]=-({\mathrm{\&omega;}}_n+T\&CenterDot;w_n[1\left]\right)(x_n+T\&CenterDot;h_n[1\left]\right)\\ w_n\left[2\right]=\mathrm{\&gamma;}\&CenterDot;\mathrm{\&mu;}\&CenterDot;e\left[2\right]\&CenterDot;(z_n+T\&CenterDot;g_n[1\left]\right)/b((x_n+T\&CenterDot;h_n[1\left]\right),(z_n+T\&CenterDot;g_n[1\left]\right))\end{array}\right.---\left(2\right)$

S103：利用式(3)获得直流分量增量h₀[3]；

$\left\{\begin{array}{c}e\left[3\right]=u\left[s\right]-\underset{m=0}{\overset{N}{\mathrm{\&Sigma;}}}(x_m+T\&CenterDot;h_m[2\left]\right)\\ h_0\left[3\right]=\mathrm{\&mu;}\&CenterDot;e\left[3\right]\end{array}\right.---\left(3\right)$

S104：令下标n依次取值为1,2,…,N，循环执行式(4)，获得交流分量增量h_n[3]、正交分量增量g_n[3]和估计频率增量w_n[3]。

$\left\{\begin{array}{c}{h}_n\left[3\right]=\mathrm{\&mu;}\&CenterDot;e\left[3\right]+({\mathrm{\&omega;}}_n+T\&CenterDot;w_n[2\left]\right)(z_n+T\&CenterDot;g_n[2\left]\right)\\ g_n\left[3\right]=-({\mathrm{\&omega;}}_n+T\&CenterDot;w_n[2\left]\right)(x_n+T\&CenterDot;h_n[2\left]\right)\\ w_n\left[3\right]=\mathrm{\&gamma;}\&CenterDot;\mathrm{\&mu;}\&CenterDot;e\left[3\right]\&CenterDot;(z_n+T\&CenterDot;g_n[2\left]\right)/b((x_n+T\&CenterDot;h_n[2\left]\right),(z_n+T\&CenterDot;g_n[2\left]\right))\end{array}\right.---\left(4\right)$

第二处理方案依次包括：

S201：先利用式(5)获得直流分量增量h₀[4]，再依据式(6)对直流分量x₀进行迭代处理；

$\left\{\begin{array}{c}e\left[4\right]=u\left[s\right]-\underset{m=0}{\overset{N}{\mathrm{\&Sigma;}}}(x_m+2T\&CenterDot;h_m[3\left]\right)\\ h_0\left[4\right]=\mathrm{\&mu;}\&CenterDot;e\left[4\right]\end{array}\right.---\left(5\right)$

$x_0\&LeftArrow;x_0+\frac{T}{3}\left(h_0\right[4]+2h_0[3]+2h_0[2]+h_0[1\left]\right)---\left(6\right)$

S202：令下标n依次取值为1,2,…,N，循环执行式(7)和式(8)，其中，每次循环，先利用式(7)获得交流分量增量h_n[4]、正交分量增量g_n[4]和估计频率增量w_n[4]，再利用式(8)分别对交流分量x_n、正交分量z_n、估计频率ω_n进行迭代处理，并对迭代后的估计频率ω_n进行限幅处理以满足条件β_2n-1≤ω_n＜β_2n；依据迭代后的交流分量x_n和正交分量z_n计算估计幅值a_n的值；

$\left\{\begin{array}{c}{h}_n\left[4\right]=\mathrm{\&mu;}\&CenterDot;e\left[4\right]+({\mathrm{\&omega;}}_n+2T\&CenterDot;w_n[3\left]\right)(z_n+2T\&CenterDot;g_n[3\left]\right)\\ g_n\left[4\right]=-({\mathrm{\&omega;}}_n+2T\&CenterDot;w_n[3\left]\right)(x_n+2T\&CenterDot;h_n[3\left]\right)\\ w_n\left[4\right]=\mathrm{\&gamma;}\&CenterDot;\mathrm{\&mu;}\&CenterDot;e\left[4\right]\&CenterDot;(z_n+2T\&CenterDot;g_n[3\left]\right)/b((x_n+2T\&CenterDot;h_n[3\left]\right),(z_n+2T\&CenterDot;g_n[3\left]\right))\end{array}\right.---\left(7\right)$

$\left\{\begin{array}{c}{x}_n\&LeftArrow;x_n+\frac{T}{3}\left(h_n\right[4]+2h_n[3]+2h_n[2]+h_n[1\left]\right)\\ z_n\&LeftArrow;z_n+\frac{T}{3}\left(g_n\right[4]+2g_n[3]+2g_n[2]+g_n[1\left]\right)\\ {\mathrm{\&omega;}}_n\&LeftArrow;{\mathrm{\&omega;}}_n+\frac{T}{3}\left(w_n\right[4]+2w_n[3]+2w_n[2]+w_n[1\left]\right)\\ a_n=\sqrt{{x_n}^2+{z_n}^2}\end{array}\right.---\left(8\right)$

S203：依据迭代处理后的直流分量x₀和交流分量x₁、x₂、…、x_N的当前值，利用式(9)获得直流分量增量h₀[1]；

$\left\{\begin{array}{c}e\left[1\right]=u\left[s\right]-\underset{m=0}{\overset{N}{\mathrm{\&Sigma;}}}{x}_m\\ h_0\left[1\right]=\mathrm{\&mu;}\&CenterDot;e\left[1\right]\end{array}\right.---\left(9\right)$

S204：令下标n依次取值为1,2,…,N，循环执行式(10)，依据迭代处理后的交流分量x_n和正交分量z_n，获得交流分量增量h_n[1]、正交分量增量g_n[1]和估计频率增量w_n[1]。

$\left\{\begin{array}{c}{h}_n\left[1\right]=\mathrm{\&mu;}\&CenterDot;e\left[1\right]+{\mathrm{\&omega;}}_n\&CenterDot;z_n\\ g_n\left[1\right]=-{\mathrm{\&omega;}}_n\&CenterDot;x_n\\ w_n\left[1\right]=\mathrm{\&gamma;}\&CenterDot;\mathrm{\&mu;}\&CenterDot;e\left[1\right]\&CenterDot;z_n/b(x_n,z_n)\end{array}\right.---\left(10\right)$

上述分母函数b在非归一化方法中，取为大于零的常数，优选地，b(x,y)=1；在归一化方法中，函数式为b(x,y)=ε+x²+y²或者b(x,y)=max(ε,x²+y²)，其中函数max表示返回两个自变量中的最大值，正数ε不大于量化误差，其目的是为在除法运算时保证除数大于零。

在上述两种实施方式中，采样周期T、交流分量个数N以及参数β₁、β₂、…、β_2N-1、β_2N的数值，均可依据被测信号的先验知识与信号分析要求设定。采样周期T首先要满足香农采样定理要求，在满足香农采样定理要求下，采样周期越小，分析精度越高。对于在线分析实施方式，受实时性的限制，交流分量个数N和采样周期T还应该满足(N+1)Δt＜T，其中Δt表示执行一遍定时中断步骤所需要的最大时间。

由于采用本发明的频率自适应信号分析方法得到的估计频率是局部收敛的，需要把每个估计频率限制在一定变化范围内，因而对估计频率进行限幅处理，这使得估计频率的物理量纲为弧度/秒。依次递增的β₁、β₂、…、β_2N保证任意两个估计频率的变化范围不会有交集。由于计算机系统中数据是有限字长的，为避免迭代过程中出现饱和，对直流分量、交流分量、正交分量也可进行限幅处理。

参数μ的物理意义相当于通频带的带宽，依据香农采样定理，限定其数值为0＜μ＜2π/T。μ数值大小对估计幅值的收敛速度具有主要影响，μ的值越大，估计幅值越快地收敛到实际值，但同时加大了干扰对幅值估计精度的不良影响。

参数γ的数值对估计频率的收敛速度具有主要影响，γ的值越大，估计频率越快地收敛到实际值，但同时加大了干扰对频率估计精度的不良影响。优选地，限定0＜γ＜10⁶。

当所有估计频率ω₁、ω₂、…、ω_N均收敛到实际值以后，直流分量x₀，交流分量x₁、x₂、…、x_N和正交分量z₁、z₂、…、z_N总能分别收敛到各自的实际值。各个估计频率的数值在经过第一次迭代处理及限幅处理后，即分别被限定在设定的容许区间内。因此对于这些迭代变量的初值没有特别限制，估计频率的初值可依据被测信号的先验知识与信号分析要求设定。

对于直流分量增量h₀[1]，交流分量增量h₁[1]、h₂[1]、…、h_N[1]，正交分量增量g₁[1]、g₂[1]、…、g_N[1]和估计频率增量w₁[1]、w₂[1]、…、w_N[1]的初值，没有特别限制。优选地，均设定为0。

设定标志字FLAG是为了区分两种不同的处理方案，只要能体现两种数值即可。优选地，标志字FLAG的值取为0和1。

若被测信号表示为u=U₀+U₁sin(Ω₁t+δ₁)+U₂sin(Ω₂t+δ₂)+…+U_Nsin(Ω_Nt+δ_N)，且对于下标n分别取值为1,2,…,N，被测信号中最多只有一个频率Ω_n总满足β_2n-1≤Ω_n＜β_2n，则经过上述方法分析之后，直流分量x₀收敛到U₀，交流分量x₁、x₂、…、x_N分别收敛到U₁sin(Ω₁t+δ₁)、U₂sin(Ω₂t+δ₂)、…、U_Nsin(Ω_Nt+δ_N)，正交分量z₁、z₂、…、z_N分别收敛到U₁cos(Ω₁t+δ₁)、U₂cos(Ω₂t+δ₂)、…、U_Ncos(Ω_Nt+δ_N)，估计频率ω₁、ω₂、…、ω_N分别收敛到Ω₁、Ω₂、…、Ω_N，估计幅值a₁、a₂、…、a_N分别收敛到U₁、U₂、…、U_N，获得高精度的分析结果。

以下结合实例说明本发明的频率自适应信号分析方法的有效性。

例如：被测信号为u=U₀+U₁sin(2πf₁t+δ₁)+U₂sin(2πf₂t+δ₂)+U₃sin(2πf₃t+δ₃)，其中各个参数随时间t的变化如表1所示：

表1被测信号各个参数随时间t变化（其中π为圆周率）

为通过变化曲线说明本发明的频率自适应信号分析方法的功效，定义信号分量跟踪误差err_x为：

err_x＝|U₀-x₀|+|U₁sin(2πf₁t+δ₁)-x₁|+|U₂sin(2πf₂t+δ₂)-x₂|+|U₃sin(2πf₃t+δ₃)-x₃|

定义交流幅值跟踪误差err_a为：

err_a＝|U₁-a₁|+|U₂-a₂|+|U₃-a₃|

定义交流频率跟踪误差err_f为：

err_f＝(|f₁-ω₁/2/π|+|f₂-ω₂/2/π|+|f₃-ω₃/2/π|)15

其中分母中的数值“15”等于三个交流分量频率的跳变量绝对值的累加。

首先以T＝0.05毫秒的采样周期对被测信号进行定时采样，获得采样数据形成离散时间序列，再按照图1所示的离线分析实施方式，编写程序在计算机中仿真运行。设定N＝3，所述正数β₁、β₂、β₃、β₄、β₅、β₆的值依次为80π、120π、180π、250π、280π、350π，设定直流分量x₀、交流分量x₁、x₂、x₃的初值均为0，设定正交分量z₁、z₂、z₃的初值均为0，设定估计频率ω₁、ω₂、ω₃的初值依次为100π、200π、300π。设定直流分量增量h₀[1]，交流分量增量h₁[1]、h₂[1]、h₃[1]，正交分量增量g₁[1]、g₂[1]、g₃[1]和估计频率增量w₁[1]、w₂[1]、w₃[1]的初值均为0。选择归一化方法的分母函数为b(x,y)=ε+x²+y²，其中ε=10^-5。设定μ和γ的值，归一化方法所得结果显示在图3、图5、图7中，非归一化方法所得结果显示在图4、图6、图8中。

图3和图4分别示出归一化方法和非归一化方法的信号分析性能。当信号分量跟踪误差err_x等于零时，说明直流分量与各个交流分量都分别收敛到其实际值，实现对被测信号的精确分析。图3说明归一化方法的信号分量收敛速度基本与被测信号的直流分量和交流幅值的大小无关，图4说明非归一化方法的信号分量收敛速度与被测信号的直流分量和交流幅值的大小有较大关系，幅值越小，收敛速度越慢。

图5和图6分别示出归一化方法和非归一化方法的幅值估计性能。若交流幅值跟踪误差err_a等于零，说明所有估计幅值都收敛到其实际值。在归一化和非归一化方法中，估计幅值的收敛速度均主要决定于参数μ的数值，参数μ的数值越大，估计幅值收敛速度越快，受参数γ的影响小。

图7和图8分别示出归一化方法和非归一化方法的频率估计性能。若交流频率跟踪误差err_f等于零，说明所有估计频率都收敛到其实际值。图7说明归一化方法的估计频率的收敛速度主要决定于受参数γ的数值，参数γ的数值越大，估计频率收敛速度越快，受参数μ的影响较小，并且基本与被测信号的直流分量和交流幅值的大小无关。图8说明非归一化方法的估计频率的收敛速度受参数γ和参数μ的影响都较大，还受到被测信号幅值大小的影响。交流频率跟踪误差err_f的收敛速度在第2秒～第3秒时之间比在第4秒～第5秒时之间要慢，这是因为被测信号的幅值在前一时段比在后一时段小。

以上实施例仅为本发明的示例性实施例，不用于限制本发明，本发明的保护范围由权利要求书限定。本领域技术人员可以在本发明的实质和保护范围内，对本发明做出各种修改或等同替换，这种修改或等同替换也应视为落在本发明的保护范围内。

Claims (6)

1.一种频率自适应信号分析方法，其特征在于，以T秒为采样周期对被测信号进行定时采样和量化得到采样数据；设定交流分量个数N的值，通过迭代直接获得直流分量x₀，通过离散频率估计算法直接获得估计频率ω₁、ω₂、…、ω_N，交流分量x₁、x₂、…、x_N和正交分量z₁、z₂、…、z_N，并计算估计幅值a₁、a₂、…、a_N，其中，针对被测信号的采样数据u[s]，按下述第一处理方案和第二处理方案中的一种执行，连续的两个采样数据分别执行不同的处理方案；其中，β₁、β₂、…、β_2N-1、β_2N为依次递增且均不大于2π/T的正数，0<μ<2π/T，0<γ<10⁶，分母函数b为大于零的二元函数；

第一处理方案依次包括以下步骤：

S101：利用式(1)获得直流分量增量h₀[2]；

$\left\{\begin{array}{c}e\left[2\right]=u\left[s\right]-\underset{m=0}{\overset{N}{\mathrm{\&Sigma;}}}(x_m+T\&CenterDot;h_m[1\left]\right)\\ h_0\left[2\right]=\mathrm{\&mu;}\&CenterDot;e\left[2\right]\end{array}\right.----\left(1\right)$

S102：令下标n依次取值为1,2,…,N，循环执行式(2)，获得交流分量增量h_n[2]、正交分量增量g_n[2]和估计频率增量w_n[2]；

$\left\{\begin{array}{c}{h}_n\left[2\right]=\mathrm{\&mu;}\&CenterDot;e\left[2\right]+({\mathrm{\&omega;}}_n+T\&CenterDot;w_n[1\left]\right)(z_n+T\&CenterDot;g_n[1\left]\right)\\ g_n\left[2\right]=-({\mathrm{\&omega;}}_n+T\&CenterDot;w_n[1\left]\right)(x_n+T\&CenterDot;h_n[1\left]\right)\\ w_n\left[2\right]=\mathrm{\&gamma;}\&CenterDot;\mathrm{\&mu;}\&CenterDot;e\left[2\right]\&CenterDot;(z_n+T\&CenterDot;g_n[1\left]\right)/b((x_n+T\&CenterDot;h_n[1\left]\right),(z_n+T\&CenterDot;g_n[1\left]\right))\end{array}\right.---\left(2\right)$

S103：利用式(3)获得直流分量增量h₀[3]；

$\left\{\begin{array}{c}e\left[3\right]=u\left[s\right]-\underset{m=0}{\overset{N}{\mathrm{\&Sigma;}}}(x_m+T\&CenterDot;h_m[2\left]\right)\\ h_0\left[3\right]=\mathrm{\&mu;}\&CenterDot;e\left[3\right]\end{array}\right.----\left(3\right)$

S104：令下标n依次取值为1,2,…,N，循环执行式(4)，获得交流分量增量h_n[3]、正交分量增量g_n[3]和估计频率增量w_n[3]；

$\left\{\begin{array}{c}{h}_n\left[3\right]=\mathrm{\&mu;}\&CenterDot;e\left[3\right]+({\mathrm{\&omega;}}_n+T\&CenterDot;w_n[2\left]\right)(z_n+T\&CenterDot;g_n[2\left]\right)\\ g_n\left[3\right]=-({\mathrm{\&omega;}}_n+T\&CenterDot;w_n[2\left]\right)(x_n+T\&CenterDot;h_n[2\left]\right)\\ w_n\left[3\right]=\mathrm{\&gamma;}\&CenterDot;\mathrm{\&mu;}\&CenterDot;e\left[3\right]\&CenterDot;(z_n+T\&CenterDot;g_n[2\left]\right)/b((x_n+T\&CenterDot;h_n[2\left]\right),(z_n+T\&CenterDot;g_n[2\left]\right))\end{array}\right.---\left(4\right)$

第二处理方案依次包括以下步骤：

S201：先利用式(5)获得直流分量增量h₀[4]，再依据式(6)对直流分量x₀进行迭代处理；

$\left\{\begin{array}{c}e\left[4\right]=u\left[s\right]-\underset{m=0}{\overset{N}{\mathrm{\&Sigma;}}}(x_m+2T\&CenterDot;h_m[3\left]\right)\\ h_0\left[4\right]=\mathrm{\&mu;}\&CenterDot;e\left[4\right]\end{array}\right.----\left(5\right)$

$x_0\&LeftArrow;x_0+\frac{T}{3}\left(h_0\right[4]+2h_0[3]+2h_0[2]+h_0[1\left]\right)---\left(6\right)$

S202：令下标n依次取值为1,2,…,N，循环执行式(7)和式(8)；其中，每次循环，先利用式(7)获得交流分量增量h_n[4]、正交分量增量g_n[4]和估计频率增量w_n[4]，再利用式(8)分别对交流分量x_n、正交分量z_n、估计频率ω_n进行迭代处理，并对迭代后的估计频率ω_n进行限幅处理以满足条件β_2n-1≤ω_n<β_2n；依据迭代后的交流分量x_n和正交分量z_n计算估计幅值a_n的值；

$\left\{\begin{array}{c}{h}_n\left[4\right]=\mathrm{\&mu;}\&CenterDot;e\left[4\right]+({\mathrm{\&omega;}}_n+2T\&CenterDot;w_n[3\left]\right)(z_n+2T\&CenterDot;g_n[3\left]\right)\\ g_n\left[4\right]=-({\mathrm{\&omega;}}_n+2T\&CenterDot;w_n[3\left]\right)(x_n+2T\&CenterDot;h_n[3\left]\right)\\ w_n\left[4\right]=\mathrm{\&gamma;}\&CenterDot;\mathrm{\&mu;}\&CenterDot;e\left[4\right]\&CenterDot;(z_n+2T\&CenterDot;g_n[3\left]\right)/b((x_n+2T\&CenterDot;h_n[3\left]\right),(z_n+2T\&CenterDot;g_n[3\left]\right))\end{array}\right.---\left(7\right)$

$\left\{\begin{array}{c}{x}_n\&LeftArrow;x_n+\frac{T}{3}\left(h_n\right[4]+2h_n[3]+2h_n[2]+h_n[1\left]\right)\\ z_n\&LeftArrow;z_n+\frac{T}{3}\left(g_n\right[4]+2g_n[3]+2g_n[2]+g_n[1\left]\right)\\ {\mathrm{\&omega;}}_n\&LeftArrow;{\mathrm{\&omega;}}_n+\frac{T}{3}\left(w_n\right[4]+2w_n[3]+2w_n[2]+w_n[1\left]\right)\\ {\mathrm{\&alpha;}}_n=\sqrt{{x_n}^2+{z_n}^2}\end{array}\right.---\left(8\right)$

S203：依据迭代处理后的直流分量x₀和交流分量x₁、x₂、…、x_N的当前值，利用式(9)获得直流分量增量h₀[1]；

$\left\{\begin{array}{c}e\left[1\right]=u\left[s\right]-\underset{m=0}{\overset{N}{\mathrm{\&Sigma;}}}{x}_m\\ h_0\left[1\right]=\mathrm{\&mu;}\&CenterDot;e\left[1\right]\end{array}\right.---\left(9\right)$

S204：令下标n依次取值为1,2,…,N，循环执行式(10)，

$\left\{\begin{array}{c}{h}_n\left[1\right]=\mathrm{\&mu;}\&CenterDot;e\left[1\right]+{\mathrm{\&omega;}}_n\&CenterDot;z_n\\ g_n\left[1\right]=-{\mathrm{\&omega;}}_n\&CenterDot;x_n\\ w_n\left[1\right]=\mathrm{\&gamma;}\&CenterDot;\mathrm{\&mu;}\&CenterDot;e\left[1\right]\&CenterDot;z_n/b(x_n,z_n)\end{array}\right.---\left(10\right)$

依据迭代处理后的交流分量x_n和正交分量z_n，获得交流分量增量h_n[1]、正交分量增量g_n[1]和估计频率增量w_n[1]。

2.根据权利要求1所述的频率自适应信号分析方法，其特征在于，当采用非归一化方法时，所述分母函数b为大于0的常数。

3.根据权利要求2所述的频率自适应信号分析方法，其特征在于，当采用非归一化方法时，所述分母函数b＝1。

4.根据权利要求1所述的频率自适应信号分析方法，其特征在于，当采用归一化方法时，所述分母函数b的函数式为b(x,y)＝ε+x²+y²或者b(x,y)＝max(ε,x²+y²)，其中正数ε不大于量化误差，函数max表示返回两个自变量中的最大值。

5.根据权利要求1所述的频率自适应信号分析方法，其特征在于，具体包括以下步骤：

S1：设定参数T、N、β₁、β₂、…、β_2N-1、β_2N、μ、γ、ε的值，设定直流分量x₀，交流分量x₁、x₂、…、x_N，正交分量z₁、z₂、…、z_N和估计频率ω₁、ω₂、…、ω_N的初值；

S2：设定直流分量增量h₀[1]、交流分量增量h₁[1]、h₂[1]、…、h_N[1]，正交分量增量g₁[1]、g₂[1]、…、g_N[1]和估计频率增量w₁[1]、w₂[1]、…、w_N[1]的初值，设定标志字的初值为执行第一处理方案；

S3：对被测信号的采样数据u[s]，依据标志字的值执行第一处理方案或者第二处理方案，然后修改标志字的值为执行另外一种处理方案；

S4：返回步骤S3循环执行，直至迭代完成。

6.根据权利要求5所述的频率自适应信号分析方法，其特征在于，将步骤S2替换为：设定直流分量增量h₀[1]、h₀[2]、h₀[3]，交流分量增量h₁[1]、h₁[2]、h₁[3]、h₂[1]、h₂[2]、h₂[3]、…、h_N[1]、h_N[2]、h_N[3]，正交分量增量g₁[1]、g₁[2]、g₁[3]、g₂[1]、g₂[2]、g₂[3]、…、g_N[1]、g_N[2]、g_N[3]和估计频率增量w₁[1]、w₁[2]、w₁[3]、w₂[1]、w₂[2]、w₂[3]、…、w_N[1]、w_N[2]、w_N[3]的初值，设定标志字的初值为执行第二处理方案。