CN110068729A

CN110068729A - 一种信号相量计算方法

Info

Publication number: CN110068729A
Application number: CN201910322409.8A
Authority: CN
Inventors: 庞吉耀
Original assignee: NANJING PANENG ELECTRIC POWER TECHNOLOGY CO LTD
Current assignee: NANJING PANENG ELECTRIC POWER TECHNOLOGY CO LTD
Priority date: 2019-04-22
Filing date: 2019-04-22
Publication date: 2019-07-30
Anticipated expiration: 2039-04-22
Also published as: CN110068729B

Abstract

本发明公开了一种信号相量计算方法，包括如下步骤：通过等间隔采样获取信号样本数据；通过时域插值构建虚拟采样点获得完整信号周期的采样信息；利用变数据窗长度的全周期数据参与DFT变换计算，结合梯形积分推导信号各次谐波DFT变换的实部和虚部的计算公式；根据预编制的查找表计算DFT变换系数；根据DFT变换系数计算信号各次谐波的实部和虚部；根据信号各次谐波的实部和虚部计算信号的幅值和相角，确定信号相量。本发明能够快速准确获得信号的各次谐波的相量参数，能够同时满足计算量小、跟踪速度快和计算精度高等要求。

Description

一种信号相量计算方法

技术领域

本发明涉及信号处理技术领域，尤其涉及一种信号相量计算方法。

背景技术

信号的幅值、相角和频率这三个特征量是反映供电系统运行特性的重要参数。及时、准确地测量系统频率及相量可以预测系统是否稳定运行，从而触发继电保护动作来保证系统的安全运行。而通过对信号的时域采样，用DFT对信号进行谱分析就可以获得信号的各次谐波的相量参数。当电网处于工频时，基于定间隔采样技术的DFT算法具有良好的谐波滤波特性，测量结果十分精确，但当电网频率偏离50hz时，由于非同步采样带来频域泄漏导致传统测量算法难以同时满足计算量小、跟踪速度快和计算精度高等要求。

为提高传统DFT算法精度，目前已提出了多种改进方法，修正传统的DFT算法，减少频谱泄漏的影响。这些算法除去在中间的步骤推导时可能会出现一些舍入误差，计算量比较大外，还因为使用多个信号周期的数据导致算法实时性下降，限制了其应用。

发明内容

本发明的目的在于克服现有技术中的不足，提供一种信号相量计算方法，解决现有技术中信号相量计算存在频谱泄漏、计算量大、实时性差的技术问题。

为解决上述问题，本发明是采用下述技术方案实现的：一种信号相量计算方法，所述方法包括如下步骤：

通过等间隔采样获取信号样本数据；

通过时域插值构建虚拟采样点获得完整信号周期的采样信息；

利用变数据窗长度的全周期数据参与DFT变换计算，结合梯形积分推导信号各次谐波DFT变换的实部和虚部的计算公式；

根据预编制的查找表计算DFT变换系数；

根据DFT变换系数计算信号各次谐波的实部和虚部；

根据信号各次谐波的实部和虚部计算信号的幅值和相角，确定信号相量。

进一步地，所述虚拟采样点的构建方法包括如下步骤：

计算信号每周期所需的最大整数采样点数N；

根据信号在(N-1)T_s、NT_s、(N+1)T_s时刻的实际采样值，计算获得时域信号x(t)在T_x处的虚拟采样值x(N_x)；

其中：x(t)为时域信号；t为时间；T_s为采样间隔；T_x为信号周期，T_x＝NT_s+Δt；Δt为时域整数倍采样泄漏时间；N_x为广义上的信号每周期采样点数，N_x＝T_x/T_s。

进一步地，所述最大整数采样点数N采用下述公式计算获取：

其中：INT()表示向下取整函数。

进一步地，所述虚拟采样值x(N_x)通过二阶内插法计算获取，具体计算公式如下：

其中：x(N-1)为时域信号x(t)在(N-1)T_s时刻的实际采样值；x(N)为时域信号x(t)在NT_s时刻的实际采样值；x(N+1)为时域信号x(t)在(N+1)T_s时刻的实际采样值。

进一步地，信号各次谐波DFT变换的实部和虚部的计算公式推导方法包括如下步骤：

根据傅立叶变换原理，若信号具有离散谱和可变信号周期T_x，则其k次谐波实部和虚部分别表示如下：

式中：为k次谐波傅里叶变换的实部；为k次谐波傅里叶变换的虚部；x(t)为时域信号；t为时间；k为谐波次数；为k次谐波角频率；

T_s为采样间隔，存在最大整数采样点数N使得下式成立：

NT_s≤T_x≤(N+1)T_s (3)

将公式(1)和公式(2)变形为：

对公式(4)和公式(5)进行离散化，进一步依据梯形积分思路将定积分使用离散求和替代，将公式(4)和公式(5)进行变形推导，并记x(nT_s)为x(n)如下：

式中：x(n)为时域信号x(t)在n*T_s时刻的实际采样值，且x(n)对应x(nT_s)，n为数据窗索引且n∈[0,N]；x(n+1)为时域信号x(t)在(n+1)*T_s时刻的实际采样值；x(0)为时域信号x(t)在0*T_s时刻的实际采样值；x(N_x)为时域信号x(t)在N_xT_s(即T_x，T_x＝N_xT_s＝NT_s+Δt)处的虚拟采样值；Δt＝T_x-NT_s，0≤Δt≤T_s；

当满足同步采样条件时，T_x＝NT_s且Δt＝0时x(N_x)＝x(N)，对公式(6)和公式(7)进行简化：

进一步地，所述查找表的编制方法包括如下步骤：

a、确定查找表每行间的频率间距f_Δ、最低索引频率f_min和最高索引频率f_max，设定谐波次数k＝1，将k作为查找表页的页地址索引，变量m作为查找表的列索引；

b、设定初始索引频率为f_b＝f_min，变量L＝0；

c、计算索引频率f_b处的最大整数采样点数

d、按下式计算当前条件下的正弦表和余弦表的一行数据，其中正弦表和余弦表当前行中每一项的计算公式如下：

其中：f_s为采样频率，f_b为索引频率，m作为查找表的列索引，相同频率下其值和信号数据窗索引对应相等；

进一步地，上述和的计算结果作为查找表项放在第k页的第L行的第m列；

e、将变量L自增为L＝L+1，重新计算索引频率f_b＝f_min+L×f_Δ转步骤(d)计算下一个频点的一行数据，直至f_b＝f_max时结束；

f、增加谐波次数k＝k+1，确定新的查找表页，转步骤(b)计算下一个谐波对应的DFT变换系数页表，直至期望的谐波查找表计算完成。

进一步地，查找表的每一页为某次谐波的DFT变换系数，最大行数为(f_max-f_min)/f_Δ，特定页每一行对应一个索引频率f_b的个三角函数值，表的列数为N_max＝f_s/f_min，当某行的实际系数数量时，将不足的部分补0。

进一步地，基于查找表的DFT变换系数和计算方法包括如下步骤：

取当前计算的谐波次数k作为页索引确定查找表的页地址table_base，根据当前信号频率f_x计算查找表的行索引row_index，row_index＝INT[(f_x-f_min)/f_Δ]，其中：f_x为当前信号频率；

取当前数据窗索引n为查找表的列索引col_index，并根据table_base和row_index通过查表获得索引频率f_b处的正余弦表项，以此作为和的计算结果；f_b＝row_index×f_Δ+f_min；

确定索引频率f_b处的频率微小增量Δf，Δf＝f_x-f_b；

计算当前的DFT系数和

1/T_x＝f_b+Δf。

与现有技术相比，本发明所达到的有益效果是：采用等间隔采样，不需要动态改变采样间隔以适应信号频率变化，可以使用传统的数字滤波技术同时保证采样系统的稳定；泄漏补偿在时域进行，计算精度高，概念清晰明确，有效地简化以往频域插值的运算，通过时域插值构建虚拟采样点得到完整信号周期的采样信息，再利用变窗长度的全周期DFT计算，快速获得信号的各次谐波的相量参数，能够同时满足计算量小、跟踪速度快和计算精度高等要求；计算仅需要一个周期左右的采样数据，内存消耗少。

附图说明

图1是适用于本发明实施例的信号时域采样和插值示意图。

具体实施方式

本发明实施例提供的信号相量计算方法，能够克服信号频率变化因素对离散谱信号的相量计算影响，适用于对计算精度和实时性要求较高，但信号频率变化的场合。本发明实施例首先通过等间隔采样获得信号一个周期左右的样本数据，通过时域插值构建虚拟采样点获得完整信号周期的采样信息，再利用变数据窗长度的全周期数据参与DFT计算，并结合梯形积分和构建基于查找表的DFT系数修正算法，快速准确获得信号的各次谐波的相量参数，能够同时满足计算量小、跟踪速度快和计算精度高等要求。

以下为本发明实施例的基本实现步骤：

1).根据信号中含有的最高次谐波频率f_{h_max}依采样定理确定采样频率f_s(f_s＞2f_{h_max})，并以固定采样间隔T_s＝1/f_s获取信号的样本数据；

2).计算信号基波频率f_x(或周期T_x)，可采用过零检测法或傅氏测频法等方法；

3).根据信号周期T_x按下面的公式计算信号每周期所需最大整数采样点数N和信号时域整数倍采样泄漏时间Δt；

T_x＝NT_s+Δt

NT_s≤T_x≤(N+1)T_s

其中：T_s为采样间隔；T_x为信号周期；Δt为时域整数倍采样泄漏时间；

4).从当前采样点开始向前选取信号N+2点采样数据；

5).由信号在(N-1)T_s、NT_s、(N+1)T_s时刻的采样值通过二阶内插获得x(t)在T_x(＝NT_s+Δt)处的虚拟采样值x(N_x)；

6).若记信号x(t)在nT_s处的采样值为x(n)，利用下面的公式求取信号指定谐波的实部虚部：

其中：k为谐波次数；N为信号每周期所需的最大整数采样点数；x(0)～x(N)为信号x(t)当前实际采样点的前N+1个值；n为数据窗索引，n∈[0,N]；x(N_x)为信号x(t)在N_xT_s(即T_x，T_x＝N_xT_s＝NT_s+Δt)处的虚拟采样值，由信号在(N-1)T_s、NT_s、(N+1)T_s时刻的值通过二阶内插获得；为k次谐波傅里叶变换的实部；为k次谐波傅里叶变换的虚部；

7).利用和计算信号各次谐波的幅值和瞬时相位。

本发明利用时域信号x(t)完整信号周期T_x上N_x个采样数据(N_x为广义上的信号每周期采样点数，且N_x＝T_x/T_s)，即包含N+1个实际采样点和一个虚拟采样点x(N_x)，x(N_x)通过时域信号x(t)在(N-1)T_s、NT_s、(N+1)T_s时刻的值按下面的公式获得：

其中：x(N-1)为信号(N-1)T_s时刻的采样值；x(N)为信号NT_s时刻的采样值；x(N+1)为信号(N+1)T_s时刻的采样值。

进一步地，DFT变换系数也使用了完整信号周期T_x。

为减少运算量和适应定点计算需要，本发明实施例进一步推演出基于查找表的DFT变换系数快速计算方法。所述查找表基于频域编制，首先计算每一个索引频率f_b处的最大整数采样点数再以下面的公式作为查找表的计算核心：

m作为查找表的列索引，相同频率下其值和数据窗索引对应相等；

进一步地在构造出查找表之后，本发明利用下面的公式计算DFT变换系数和

1/T_x＝f_b+Δf。

其中:f_b为查找表索引频率；Δf为索引频率f_b处的频率微小增量；和为对应f_b和n的正余弦计算结果，本发明利用前述查找表以f_b为索引频率，同时以n为列索引，并取对应查找表项的值作为和的计算结果。

进一步地，利用修正的DFT变换系数和N+1点的实际采样数据及虚拟采样点x(N_x)，按照下面公式分别计算信号k次谐波的实部和虚部再根据实部和虚部在利用传统方法获得信号的幅值和相角。

下面结合附图对本发明作进一步描述。以下实施例仅用于更加清楚地说明本发明的技术方案，而不能以此来限制本发明的保护范围。

1.基础计算公式推导

1.1.全周期DFT变换公式推演

依傅立叶变换原理，若信号具有离散谱和可变周期T_x，则其k次谐波实部和虚部可表示如下：

参考附图1，若以采样间隔T_s对信号进行等间隔采样，将不能保证T_x＝NT_s，但一定存在整数N使的下式成立

NT_s≤T_x≤(N+1)T_s (3)

因此可以将式1a和式2a变形为:

对式1b和式2b进行离散化，将定积分使用离散求和替代，得到全周期DFT计算表达式

1.2.基于梯形积分的改进

为提升计算精度，本发明进一步推导使用梯形积分替代传统DFT的计算公式，参考附图1，将式4和式5进行变形推导(也可基于式1b和2b直接推导)并记x(nT_s)为x(n)如下：

其中：x(N_x)对应为x(t)在T_x(＝NT_s+Δt)处的虚拟采样值。

1.3.虚拟采样值计算

对照附图1可知，在非同步采样条件下，x(t)在T_x(＝NT_s+Δt)处的值x(N_x)不能通过直接采样获得，但可以通过信号在(N-1)T_s、NT_s、(N+1)T_s时刻的采样值利用二阶内插法获得，在误差允许的情况下，甚至可直接使用(N+1)T_s和NT_s时刻的值通过线性内插得到；精度要求更高的可以使用高阶内插。内插方法取决于计算量和精度的平衡，本发明实施例仅仅展示了通过二阶内插获得虚拟采样值的一种思路，二阶以上的内插对精度贡献不是很大。

若记t₀＝(N-1)T_s、t₁＝NT_s、t₂＝(N+1)T_s，并记x(nT_s)为x(n)，则推导二阶内插的计算表达式如下：

x(NT_s+Δt)＝k_ax[(N-1)T_s]+k_bx(NT_s)+k_cx[(N+1)T_s]

按式8，当满足同步采样条件时，T_x＝NT_s且Δt＝0时x(N_x)＝x(N)，式6和式7简化为下面的式6a和式7a。

不难验证，对于直流分量式7a的计算结果为0，而式6a的计算结果将等于直流分量且不损失精度。再考虑到在同步采样条件满足时周期信号x(0)＝x(N)，此时式6a和式7a的计算结果与传统算法一致。因此，式6和式7为实部和虚部计算的通式。

但在非同步采样条件下由于x(0)≠x(N)、T_x≠NT_s、且Δt≠0，本发明使用梯形积分优势会体现出来，其计算精度将进一步提高。

1.4.周期采样点数N计算

本发明假定信号周期已经事先测得，设信号周期为T_x(对应频率f_x＝1/T_x)，采样间隔为T_s(对应采样率f_s＝1/T_s)，则每周期得到的采样数据点数N_x为:

N_x＝T_x/T_s＝f_s/f_x (9)

对非同步采样来说N_x为小数，但用下面简单方法可以确定每个基频周期最大整数采样点数N:

N＝INT(T_x/T_s)＝INT(f_s/f_x) (N≤N_x≤N+1) (10)

上式中，函数INT(x)表示向下取整。显然满足同步采样条件时N_x＝N为整数。

1.5.基于查找表的三角函数计算

在上述表达式6和7中包含大量的三角函数计算，为减少计算量，本发明进一步推演出满足前述计算公式特点的DFT变换系数的快速计算方法和构建查找表项的计算核。设f_b为查找表的索引频率，是频率为f_b时的最大整数采样点数，DFT变换系数和的推演变形过程如下所示

Δf＝f_x-f_b (13)

其中：k为谐波次数，Δf为索引频率f_b处微小增量；

在保证足够小的前提下，进一步对式11和式12展开，并近似为：

式14和式15表明，DFT变换系数和可以利用事先计算出的和通过乘加及乘减运算获得，说明构建查找表是可行的，而且减少频率间距还可以提高公式14和公式15的精度。若以查找表每行间的频率间距为f_Δ＝1/M HZ(M为整数)等分频率区间[f_min,f_max]，对于任意频率f_x∈[f_min,f_max]，则存在整数L使得下式成立：

f_x＝f_b+Δf＝f_min+L×f_Δ+Δf Δf∈[0,f_Δ] (16)

进一步地，以＜f_b,L,m＞三元组为索引构建查找表，其中索引频率f_b和整数L相对应，具体实施时，L对应为查找表的行索引，m对应查找表的列索引，而则对应频率为f_b时的最大整数采样点数，查找表的表项对应采用和计算。

进一步地，为了计算谐波的相量，本发明为谐波构建和基波结构一样的查找表，称之为查找表页，使用以谐波次数k作为查找表的页索引。

作为本发明的另一种实施例是：不采用DFT变换系数查找表，而直接在线计算公式11和公式12，此种办法运算量增加，但没有查找表的存储开销。当查找表的频率间距足够小时，使用在线计算DFT系数将没有精度方面的优势，本发明在实际实施时需要平衡考虑系统目标，“存储优先”的离线计算可以采用本发明的公式11和公式12，“速度优先”的在线计算可以采用本发明的公式14和公式15。

1.6.三角函数查找表制作

查找表编制时以m为查找表的列索引变量，结合式14和式15，得到正余弦表项的计算核和具体的查找表的编制方法如下：

(1)、根据查找表每行间的频率间距f_Δ，并确定最低索引频率f_min和最高索引频率f_max，设定谐波次数k＝1，k作为查找表页的页地址索引；

(2)、设定初始索引频率为f_b＝f_min，变量L＝0；

(3)、计算索引频率f_b处的最大整数采样点数

(4)、按下式计算当前条件下的正弦表和余弦表的一行数据，其中正弦表和余弦表当前行中每一项的计算公式如下：

其中：f_s为采样频率，f_b为索引频率，m作为查找表的列索引；

(5)、将变量L自增为L＝L+1，重新计算索引频率f_b＝f_min+L×f_Δ转步骤(4)计算下一个频点的一行数据，直至f_b＝f_max时结束；

(6)、增加谐波次数k＝k+1，确定新的查找表页，转步骤(2)计算下一个谐波对应的DFT变换系数页表，直至期望的谐波查找表计算完成。

2.应用举例

以市电交流信号为例，以采样率f_s＝4000次/秒(T_s＝250μS)对信号进行等间隔采样时，算法操作步骤如下：

(1).取频率步距f_Δ＝0.03125Hz，并设系统最低索引频率为：f_min＝45Hz，最高索引频率为：f_max＝55Hz，构造DFT变换基波(对应k＝1)系数查找表；

(2).根据采样信号计算频率，假设所得频率为f_x＝49.5hz；

(3).获得信号每周期最大整数采样点个数及所需序列长度：

f_s/f_x＝4000/49.5＝80.8080，向下取整得N＝INT(f_s/f_x)＝80；

(4).计算虚拟采样点位置，Δt_x＝1/f_x-NT_s＝0.8080μS；

(5).从当前采样点向前，取82点数据构建数据窗，分别取前79，80、81三点数据按式8计算虚拟采样点80.808位置的虚拟采样值x(N_x)；

(6).计算查找表行索引，row＝INT[(f_x-f_min)/f_Δ]＝144，并在144行提取三角函数表中基准值；

(7).按式(14)和式(15)修正实际频率f_x下正余弦函数的准确值；

(8).按式(6)和式(7)计算信号基波分量的实部和虚部，并由此计算幅值和相角；

(9).类似的方法计算所需要的谐波实部及虚部，并导出谐波的幅值和相角。

3.应用分析

3.1.误差分析

本发明通过在时域虚拟采样和经过改进的计算公式6和公式7，实现全周期信息参与变换，成功解决时域采样泄漏带来的精度损失，进一步通过梯形积分取代传统DFT算法提高计算精度，但以下因素仍然是计算误差的主要来源：

(1).数据序列的{x(n)}采样精度，取决于ADC的精度；

(2).采样间隔T_s(对应采样率f_s＝1/T_s)，取决于ADC和采样系统性能；

(3).虚拟采样值的精度，取决于内插方式；

(4).DFT变换系数，取决于系统的计算能力；

上述几个方面因素直接决定了计算结果精度，但这些因素和系统成本和复杂度相关，系统设计时应综合考虑，如提高ADC的解析度会直接增加系统成本，提高采样率会导致运算量成倍增加。

另外，对DFT变换系数采用在线计算可以减少查找表的存储空间，但和本发明基于查找表的计算方法相比，对计算精度帮助不大但会导致运算量大幅增加。另一方面，阶数更高的内插方式同样对计算精度贡献有限。

3.2.开销分析

目前的同类算法为了达到所需要的精度至少需要两个周期以上的数据，也即动态响应能力更差。同目前的频率修正DFT算法相比，本发明需要的数据存储开销小，在同步采样条件下只需要保留信号的一个周期采样数据，即使在非同步采样条件下也只需要额外增加两个采样数据。和同样采用短数据窗变长DFT算法相比运算量要小，此前的算法需要分别计算N点DFT和N+1点DFT再进行频域内插，而且本发明的计算精度提高1个数量级。

此外本发明采用的是矩形窗直接截断，DFT计算系数利用查找表并进行精度修正，内插算法仅仅对最后虚拟采样点x(N_x)进行，因此本发明给出的算法在精度、运算量和动态响应方面取得比较好的平衡。

本发明从采用短矩形窗截取采样序列的一个周期数据，通过变长度DFT和时域虚拟采样补偿实现频域泄漏修正，和同类方法相比具有明显的优势：具体如下：1.泄漏补偿在时域进行，计算精度高，概念清晰明确；2.计算仅需要一个周期左右的采样数据，内存消耗少；3.采用短数据窗，运算量和常规DFT算法接近，时延比同类算法好；4.全周期数据参与DFT变换，但是精度较常规算法高1个数量级；5.通过增加采样率(减小采样间隔)，可以增加计算精度但不影响动态特性；6.计算结果不受信号频率的影响；7.采用等间隔采样，不需要动态改变采样间隔以适应信号频率变化，可以使用传统的数字滤波技术同时保证采样系统的稳定。

以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明技术原理的前提下，还可以做出若干改进和变形，这些改进和变形也应视为本发明的保护范围。

Claims

1.一种信号相量计算方法，其特征在于，所述方法包括如下步骤：

通过等间隔采样获取信号样本数据；

根据预编制的查找表计算DFT变换系数；

根据DFT变换系数计算信号各次谐波的实部和虚部；

2.根据权利要求1所述的信号相量计算方法，其特征在于，所述虚拟采样点的构建方法包括如下步骤：

计算信号每周期所需的最大整数采样点数N；

3.根据权利要求2所述的信号相量计算方法，其特征在于，所述最大整数采样点数N采用下述公式计算获取：

其中：INT()表示向下取整函数。

4.根据权利要求2所述的信号相量计算方法，其特征在于，所述虚拟采样值x(N_x)通过二阶内插法计算获取，具体计算公式如下：

5.根据权利要求1所述的信号相量计算方法，其特征在于，信号各次谐波DFT变换的实部和虚部的计算公式推导方法包括如下步骤：

T_s为采样间隔，存在最大整数采样点数N使得下式成立：

NT_s≤T_x≤(N+1)T_s (3)

将公式(1)和公式(2)变形为：

式中：x(n)为时域信号x(t)在n*T_s时刻的实际采样值，且x(n)对应x(nT_s)，n为数据窗索引且n∈[0,N]；x(n+1)为时域信号x(t)在(n+1)*T_s时刻的实际采样值；x(0)为时域信号x(t)在0*T_s时刻的实际采样值；x(N_x)为时域信号x(t)在N_xT_s处的虚拟采样值；T_x＝N_xT_s,Δt＝T_x-NT_s，0≤Δt≤T_s；

6.根据权利要求1所述的信号相量计算方法，其特征在于，所述查找表的编制方法包括如下步骤：

a、确定查找表每行间的频率间距f_Δ、最低索引频率f_min和最高索引频率f_max，设定谐波次数k＝1，将k作为查找表页的页地址索引；

b、设定初始索引频率为f_b＝f_min，变量L＝0；

c、计算索引频率f_b处的最大整数采样点数

7.根据权利要求6所述的信号相量计算方法，其特征在于，查找表的每一页为某次谐波的变换系数，最大行数为(f_max-f_min)/f_Δ，特定页每一行对应一个索引频率f_b的个三角函数值，表的列数为N_max＝f_s/f_min，当某行的实际系数数量时，将不足的部分补0。

8.根据权利要求6所述的信号相量计算方法，其特征在于，DFT变换系数的计算方法包括如下步骤：

取当前述数据窗索引n为查找表的列索引col_index，并根据table_base和row_index通过查表获得索引频率f_b处的正余弦表项，分别作为和的计算结果；f_b＝row_index×f_Δ+f_min；

确定索引频率f_b处的频率微小增量Δf，Δf＝f_x-f_b；

计算当前的DFT系数和

1/T_x＝f_b+Δf。