CN102435844A

CN102435844A - 一种频率无关的正弦信号相量计算方法

Info

Publication number: CN102435844A
Application number: CN2011103406808A
Authority: CN
Inventors: 庞吉耀
Original assignee: NANJING PANENG ELECTRIC POWER TECHNOLOGY CO LTD
Current assignee: NANJING PANENG ELECTRIC POWER TECHNOLOGY CO LTD
Priority date: 2011-11-01
Filing date: 2011-11-01
Publication date: 2012-05-02
Anticipated expiration: 2031-11-01
Also published as: CN102435844B

Abstract

一种频率无关的正弦信号相量计算方法，1)综合信号频率和计算精度确定采样频率f_s，并以固定采样间隔T_s＝1/f_s获取信号的样本数据；2)计算信号基波频率f_x(或周期)，采用过零检测法或傅氏测频法方法；3)根据信号频率(或周期)计算出信号的一个基频周期所需的最大采样点数N；N＝floor(f_s/f_x)，f_s-采样率，f_x-信号频率，floor-将小数向下取整的函数；4)分别选取信号N点和N+1点采样数据，计算信号N点DFT和N+1点DFT；5)利用两次DFT计算结果按下面两种表达式计算信号真实频点f_x处的实部和虚部；最后计算出信号各次谐波的幅值和瞬时相位。本发明适用于计算精度和实时性要求较高场合，但信号频率随机变化的场合。

Description

一种频率无关的正弦信号相量计算方法

技术领域

本发明属于数字信号处理和电力自动化领域，用来计算正弦信号的相量。本发明适用于对计算精度和实时性要求较高，但信号频率有随机变化的场合。

背景技术

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

为提高传统DFT算法精度，目前已提出了多种改进方法，修正传统的DFT算法，减少频谱泄漏的影响。这些算法除去在中间的步骤推导时可能会出现一些舍入误差，计算量比较大外，还因为使用多个信号周期的数据导致算法实时性下降，限制了其应用。本发明提出了一种基于频率的跟踪算法，有效地简化以往频域插值的运算，利用两次DFT结果进行线性内插，快速获得信号的各次谐波的相量参数，能够同时满足计算量小、跟踪速度快和计算精度高等要求。

目前研究较多的方法主要包括以下几类：

1.使用加窗函数和内插方法来减小泄漏误差^[5][6][7][8]

在非同步采样情况下，对时域信号加窗截断，将会导致频谱泄漏。为了提高计算精度，可以通过选择不同的窗函数对采样数据进行加权，使信号能量尽可能落在主瓣内从而减少谱泄漏，再对加窗后的序列进行FFT运算并进行双谱线插值可得到精度较高的时频参数。

2.非同步采样点的同步化^[4]

基于双速率的软件同步采样或在已知信号基频的情况下通过对原始采样信号进行拉格朗日插值，得到近似的同步化序列，然后应用DFT算法计算得到该信号各项时频参量

3.动态调整采样率实现同步采样

采用锁相环技术实现硬件同步采样或根据当前频率实时调整对信号的采样间隔，然后应用DFT算法计算得到该信号各项时频参量，可以大大地减小由非同步采样所导致的截断误差和周期误差。

以上提到的几种方法都可以在一定程度上减小频谱泄漏但都有不足之处：

1.加窗内插法为保证精度需要用到多个信号周期的数据，不但运算量大而且算法的实时性和动态响应性能不好

2.非同步采样点的同步化这类方法需要在时域对采样序列进行同步化插值，由于频率的随机变化，导致时变插值器的运算量比较大在当信号频率偏差过大时会发生插值点的跑位，插值公式这时会产生很大误差。而且算法以信号过零点作为理想采样的起点，当信号谐波过大出现多个过零点时将出现插值错误。

3.由于***频率在一定范围内随机变化，采用等角度的策略需要动态调整采样间隔，这对牵涉到多个设备的同步采样***来说会引发连锁反应导致***不稳定，传统的数字滤波器也无法适用。而实际情况中，A/D采样的速度无法如此准确地满足同步采样，尤其是在信号受到干扰发生频率波动的时候。

参考文献

1.刘建林，陈兵一种高精度单频信号频率估计算法无线电工程第2011卷第4期

2.闵勇丁仁杰韩英铎等自适应调整采样率的相量在线测量算法研究电力***及其自动化学报2010年4月

3.吴峻吴崇昊基于Kaiser窗的改进傅氏测频算法江苏电机工程2008年7月

4.曾泽昊余有灵许维胜基于插值同步算法的频谱泄漏分析与仿真***仿真技术2005年10月

5.黄纯，江亚群谐波分析的加窗插值改进算法中国电机工程学报2005年8月

6.李新福刘教民崔玉龙电器检测的谐波分析算法河北工业大学学报2003年10月

7.庞浩，李东霞，俎云霄，等应用FFT进行电力***谐波分析的改进算法[J].中国电机工程学报2003，23

8.梁志国，孙宇.信号周期的一种数字化测量方法[J].仪器仪表学报，2003，24(4，增刊)：195-198

9.高云鹏滕召胜卿柏元基于Kaiser窗双谱线插值FFT的谐波分析方法电子学报2000年12月

10.杨洪耕惠锦侯鹏电力***谐波和间谐波检测方法综述电力***自动化1998年10月

发明内容

本发明的目的是提出一种不受频率影响的电力***正弦信号的实时相量计算方法。本发明提出了一种基于频率的跟踪算法，有效地简化以往频域插值的运算，利用两次DFT结果进行线性内插，快速获得信号的各次谐波的相量参数，能够同时满足计算量小、跟踪速度快和计算精度高等要求。

本发明是通过这样的方案实现的：一种频率无关的电力***正弦信号的相量计算方法，

1).综合信号频率和计算精度确定采样频率f_s，并以固定采样间隔T_s＝1/f_s获取信号的样本数据；

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

3).根据信号频率(或周期)计算出信号的一个基频周期所需的最大采样点数N；

N＝floor(f_s/f_x)f_s-采样率 f_x-信号频率 floor-将小数向下取整的函数

4).分别选取信号N点和N+1点采样数据，计算信号N点DFT和N+1点DFT；

5).利用两次DFT计算结果按下面两种表达式计算信号真实频点f_x处的实部和虚部；

y_{x} = \frac{y_{2} \cdot (f_{x} - f_{1}) + y_{1} \cdot (f_{2} - f_{x})}{f_{2} - f_{1}}

或

y_{x} = y_{2} \cdot \frac{N}{N_{x}} \cdot (N + 1 - N_{x}) + y_{1} \cdot \frac{N + 1}{N_{x}} (N_{x} - N)

N_x＝f_s/f_x N＝floor(f_s/f_x)

式中：N_x-每周期采样点数，N-每周期最大的整数采样点数

f₂对应采用N点DFT计算时所得基频(f₂＝1/NT_s)

f₁对应采用N+1点DFT计算时所得基频(f₁＝1/[(N+1)·T]))

f_s-采样率 f_x-信号频率 floor-将小数向下取整的函数

y₁-对应信号N点DFT结果中基波或谐波的实部或虚部

y₂-对应信号N+1点DFT结果中基波或谐波的实部或虚部计算结果

y_x-插值所得最终基波或谐波的实部或虚部

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

本发明的有益效果是：电力***正弦信号实时相量的采集是电网稳定运行中重要和基础的参数，需要快速精确的获得实时相量是极为重要的，

1.本发明计算仅需要一个周期的采样数据，在最短时间内(一个信号周期)即可准确获得信号的相量值，比同类算法动态响应特性好；

2.算法简单，可以直接使用DFT或者通过对采样点补零后利用FFT算法；

3.内存消耗少，采样数据序列仅需要存储一个周期左右的数据；

4.通过增加采样率(减小采样间隔)，可以增加计算精度但不影响动态特性；

5.计算结果不受信号频率的影响

6.采用等间隔采样，不需要动态改变采样间隔以适应信号频率变化，可以使用传统的数字滤波技术同时保证采样***的稳定。

附图说明

图1是本发明信号相量插值示意图

图2是本发明信号相量的实部/虚部插值示意图

表1是本发明误差分析对照表，表中采样率单位为：采样点/秒，误差均为相对误差。

具体实现方式

1.算法的基本原理

按照数字信号处理原理，信号的N点DFT变换在频域分辨率为f_s＝1/(N·T_s)，若取序列长度为信号的一个周期，其对应为信号的基波频率为f_xh＝1/(N·T_s)。同理，N+1点DFT变换在频域分辨率为f_s＝1/[(N+1)·T])，对应信号频率为f_xl＝1/[(N+1)·T])。对于基波频率为f_x信号，总能找到合适的序列长度N，使得下式成立：

f_xl＜f_x＜f_xh 式中：f_xh＝1/NT_s，f_xl＝1/[(N+1)T_s]

当信号的频率在[f_xl，f_xh]区间内变化时，信号的相量值(及对应的实部/虚部)构成一个以频率f_x为变量的连续曲线，频点f_x处的相量则可以通过在频率区间[f_xl，f_xh]线性内插获得，而闭区间[f_xl，f_xh]的大小将决定插值的误差。

2.信号频率计算

本发明提出的算法需要先测得出信号频率，最简单的方法是通过检查信号序列中的两个过零点的位置来获得信号周期进而推算出信号频率，其基本原理如下：

(1).搜索信号序列x(n)，若发现有x(n)≤0且x(n+1)≥0，则在采样时刻t_n和t_n+1之间存在一个过零点；

(2).在时间轴上通过拉格朗日插值，得到当前过零点的时刻t_cur(t_n≤t_cur≤t_n+1)；

(3).通过计算当前零点和前一个零点t_pre之间的时间间距，可得信号周期T_x＝t_cur-t_pre

(4).记录当前零点位置并转步骤1进行新的零点搜索

在进行零点搜索时为减少噪声和谐波的影响可以对序列进行带通滤波。信号频率的计算当然也可以利用文献^[1][3][8]中提出的方法，如傅氏测频法。

3.最大周期采样点数N

设信号基频周期为T_x(对应频率f_x＝1/T_x)，采样间隔为T_s(对应采样率f_s＝1/T_s)，则每周期得到的采样数据点数为：

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

当满足同步采样条件时N_x为整数，没有泄漏。但由于信号频率扰动，将导致采样不满足每周期整数个点(即N_x不为整数)的要求，若仍以整数点计算DFT将至频域泄漏。

算法需要确定每个基频周期最大整数采样点数N，按下面的方法计算：

N＝floor(T_x/T_s)＝floor(f_s/f_x) (2)

上式中，函数floor(x)表示向下取整，显然下式成立：

N≤N_x≤N+1

(3)

4.DFT计算

本发明简单地使用矩形窗截取信号采样值序列的N个采样数据进行计算，依离散傅立叶变换(DFT)的定义，可直接用下面的表达式计算信号N点DFT

\begin{matrix} X (k) = \frac{2}{N} Σ_{n = 0}^{N - 1} x (n) \cdot w (k) & w (k) = e^{- j 2 kπ / N} & 0 \leq k \leq N - 1 \end{matrix} - - - (4)

当需要分析的谐波次数少运算量不大时，可直接利用上式计算。类似的方法可以计算长度N+1点DFT。当需要分析的谐波次数多时，可以在将N序列长度之后补0，使补0后的长度达到2的某个幂次再使用FFT计算。

5.线性插值

本发明采用线性内插的方法来获得频点f_x处的信号相量。如果只对信号的幅值和相角感兴趣可按附图1，直接插值获得频点f_x处的幅值和角度。若同时需要信号频点f_x处实部和虚部可按附图2，对相量的实部和虚部分别插值。

频点f_x处的相量的y_x值可以通过拉格朗日插值多项式获得，即：

y_{x} = \frac{y_{2} \cdot (f_{x} - f_{1}) + y_{1} \cdot (f_{2} - f_{x})}{f_{2} - f_{1}} - - - (5)

式中：f₂对应采用N点DFT计算时基频所对应的f_xh，f₁则对应采用N+1点DFT计算时基频所对应f_xl。

转化为采样点数的插值公式为

y_{x} = \frac{y_{2} \cdot (\frac{1}{N_{x} T_{s}} - \frac{1}{(N + 1) T_{s}}) + y_{1} \cdot (\frac{1}{{NT}_{s}} - \frac{1}{N_{x} T_{s}})}{\frac{1}{{NT}_{s}} - \frac{1}{(N + 1) T_{s}}}

T_s--采样间隔 (6)

上式化简为：

y_{x} = y_{2} \cdot \frac{N}{N_{x}} \cdot (N + 1 - N_{x}) + y_{1} \cdot \frac{N + 1}{N_{x}} (N_{x} - N) - - - (7)

式(7)表明，当满足同步采样条件时必有N_x＝N+1(对应：f_x＝f₁，y_x＝y₁)或者N_x＝N(对应：f_x＝f₂，y_x＝y₂)，即满足同步采样条件时算法等同于传统DFT。

根据式(4)，信号的k次谐波对应于DFT结果的第k根谱线，所以可以通过频域内插的方法一并获得对应于kf_x处的信号各次谐波参数。

6.幅值和相角计算

根据式(4)，由DFT计算结果是复向量，即相量的实部和虚部，且随着k的取值不同结果对应k次谐波。为得到信号在频点f_x处的相量，可以对DFT结果的实部和虚部用式(7)分别进行线性插值，然后从插值的结果中利用复数取模等方法获得幅值和相角。

7.应用举例

此处以市电交流信号为例，对其他频率信号操作步骤相同只是采样率会有所不同。设采样率f_s＝3200次/秒(采样间隔为T_s＝312.5μS，信号频率50hz时为64个样点每周期)，本发明的实际操作步骤如下：

(1).以采样间隔T_s对信号进行等间隔采样

(2).利用2所给出的方法计算信号的当前频率，此处假设所得频率为f_x＝49.5hz；

(3).获得信号每周期最大整数采样点个数及所需序列长度f_s/f_x＝3200/49.5＝64.6464，向下取整得N＝floor(f_s/f_x)＝64

(4).从当前采样点向前，分别取64点和65点数据构建数据窗，并计算64点和65点的DFT；

(5).按式(7)根据步骤4的结果计算信号基波分量的实部和虚部，并由此计算幅值和相角；

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

8.误差分析和对策

从图1和图2中可以看出，当f₁和f₂之间的间距越小则误差越小，而f₁和f₂则和序列长度有关。增加每周期采样点数(对应减小采样间隔)，可以减小f₁和f₂之间的间距，导致提高精度，但增加采样点数会增加运算量和存储开销，应用时应该根据精度和***运算能力做权衡。

在实际ADC采样率无法增加的场合可以在计算前使用Sinc插值法进行时域内插，对原始采样序列每两个点之间进行内插增加采样率，其结果相当于增加每周期采样点数，从而达到提高精度，该方法没有前文所述非同步采样点同步化的问题。

利用Matlab对本算法进行误差扫描，分析时幅度归一化为1(此时得到幅值误差即相对误差)，以频率和初相角为变量，采样间隔T_s为参变量按下式构建含有基波和2、3、5谐波正弦序列，再利用上述算法计算信号的幅度和相位与理论值进行比较误差。

x(n)＝sin(2πf_xnT_s+φ₁)+sin(4πf_xnT_s+φ₂)+sin(6πf_xnT_s+φ₃)+sin(10πf_xnT_s+φ₅)

式中：T_s-采样间隔，f_x-信号频率，φ₁，φ₂，φ₃，φ₅-谐波初相角

表1给出电力***常用的几种采样率下的误差数据，结果表明即使在采样率3200次/秒的条件下依然可以满足电力***对测控装置的精度要求，由于算法只使用信号单个周期的数据，所以动态响应实时性也能很好满足保护装置的要求。

表1 不同采样率下45hz-55hz正弦信号算法仿真误差对照表

采样率	3200	4000	6400	8000	10000	12800
							相位误差	0.001149	0.000728	0.000291	0.0001858	0.0001189	0.0000715
幅度误差	0.0011	0.0006959	0.0003	0.0001769	0.000113151	0

Claims

1.一种频率无关的正弦信号相量计算方法，其特征是步骤如下：

N＝floor(f_s/f_x)f_s-采样率 f_x-信号频率 floor-将小数向下取整的函数；

y_{x} = \frac{y_{2} \cdot (f_{x} - f_{1}) + y_{1} \cdot (f_{2} - f_{x})}{f_{2} - f_{1}}

或

y_{x} = y_{2} \cdot \frac{N}{N_{x}} \cdot (N + 1 - N_{x}) + y_{1} \cdot \frac{N + 1}{N_{x}} (N_{x} - N)

N_x＝f_s/f_x N＝floor(f_s/f_x)

式中：N_x-每周期采样点数，N-每周期最大的整数采样点数

f₂对应采用N点DFT计算时所得基频(f₂＝1/NT_s)

f₁对应采用N+1点DFT计算时所得基频(f₁＝1/[(N+1)·T]))

f_s-采样率 f_x-信号频率 floor-将小数向下取整的函数

y₁-对应信号N点DFT结果中基波或谐波的实部或虚部

y_x-插值所得最终基波或谐波的实部或虚部

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

2.根据权利要求1所述的频率无关的正弦信号相量计算方法，其特征是：信号每周期的最大整数采样点个数N按下面方法计算：

N＝floor(f_s/f_x)f_s-采样率 f_x-信号频率 floor-将小数向下取整的函数。

3.根据权利要求1所述的频率无关的正弦信号相量计算方法，其特征是：分别计算信号N点和N+1点的DFT，并将其结果按下面两种表达式计算信号真实频点f_x处的实部和虚部，再根据实部和虚部在利用传统方法获得信号的幅值和相角；

y_{x} = \frac{y_{2} \cdot (f_{x} - f_{1}) + y_{1} \cdot (f_{2} - f_{x})}{f_{2} - f_{1}}

或

y_{x} = y_{2} \cdot \frac{N}{N_{x}} \cdot (N + 1 - N_{x}) + y_{1} \cdot \frac{N + 1}{N_{x}} (N_{x} - N)

N_x＝f_s/f_x N＝floor(f_s/f_x)

式中：N_x-每周期采样点数，N-每周期最大的整数采样点数；

f₂对应采用N点DFT计算时所得基频(f₂＝1/NT_s)；

f₁对应采用N+1点DFT计算时所得基频(f₁＝1/[(N+1)·T]))；

f_s-采样率 f_x-信号频率 floor-将小数向下取整的函数；

y₁-对应信号N+1点DFT结果中基波或谐波的实部或虚部；

y₂-对应信号N点DFT结果中基波或谐波的实部或虚部计算结果

y_x-插值所得最终基波或谐波的实部或虚部。