2. 北京化工大学 信息科学与技术学院, 北京 100029
2. College of Information Science and Technology, Beijing University of Chemical Technology, Beijing 100029, China
风电、光伏发电等新能源有效降低了化石能源消耗对环境的污染。新能源技术一方面促进了电力电子器件在电网系统中的应用，使得系统控制、保护等更为准确及稳定；另一方面，电力电子器件设备也为电网注入了大量的谐波[1-2]，使得对系统谐波监测、电能计量中的电压与电流有效值及有功、无功电参量的测量更为复杂，即需要慎重考虑谐波及间谐波影响。其中对电压基波分量的准确测量是其他电参量测量的基础和关键。
作为常用的谐波分析方法，离散傅里叶变换(DFT)方法需考虑频谱泄露(长泄露和短泄露)及栅栏效应对测量结果的影响，特别是当电网信号为非平稳情况时[3]。作为减少频谱泄露的方法，多种余弦窗插值DFT方法利用频域谱线之间的关系，通过窗函数推导对应插值公式、估计被测信号频率[4-7]。基于迭代DFT方法也被用于信号参数估计，代价是计算量大[8]；而采用递归DFT(RDFT)可以有效降低计算量[9]，缺点是会产生较大的累积误差。一般通过编码等算法来减小误差[10]，但这种方法会导致计算度复杂，增加了额外的计算负担[11]。采用固定窗口长度的RDFT算法[12]可将信号分解为正交的频率分量，再利用高阶FIR低通滤波器(LPF)消除谐波引起的负面效应[13-14]，然而却未考虑直流分量对算法所带来的误差影响。
众所周知，非同步采样是引起DFT频谱泄露的主要原因[15-16]。一方面，DFT无法规避窗函数主瓣宽度和旁瓣衰减速度对谐波分析的干扰与影响，以及主瓣宽度带来的栅栏效应；另一方面，存在部分插值算法推导复杂、计算量大的情况；此外，为进一步避免低阶谐波所引起的泄漏干扰，需要采用3个DFT谱线的拉格朗日插值方法[17]，这样又增加了计算量。
为解决非同步采样情况下，数据截断造成DFT频谱泄露并导致误差增大，从而不能有效地从电压基波分量中自适应地去除稳态或时变谐波的问题，本文首先建立递归DFT与归一化带通滤波器的传递函数；其次，推导非整数延迟滤波器的DFT插值公式，设计了DFT插值的非整数延迟带通滤波器；最后，对设计的非整数延迟带通滤波器的性能进行仿真分析，并给出电压基波(谐波)幅度测量算法。
1 递归DFT延迟滤波器的传递函数以采样频率fs对电压信号x(t)进行离散采样，采样点数为N，对采样后的信号x(n)进行离散傅里叶变换，在第nT采样时刻，信号频谱可表示为[18]
| $ {X_n}\left( m \right) = \sum\limits_{l = n - N + 1}^n {x\left( l \right){{\rm{e}}^{ - {\rm{j2 \mathsf{ π} }}ml/N}}} $ | (1) |
式中，m= 0, 1, 2，…, N-1为频率指数，N为截断窗口内的电压样本数据个数，l为采样信号x(n)的采样时刻。在DFT中，设N随基波周期k变化，则DFT时间截断窗口大小fs/N是频率自适应的，且分辨率Δf=fs/kN。
类似于式(1)，在第(n－1)T采样时刻，电网电压信号的频谱可以表示为
| $ {X_{n - 1}}\left( m \right) = \sum\limits_{l = n - N}^{n - 1} {x\left( l \right){{\rm{e}}^{ - {\rm{j2 \mathsf{ π} }}lm/N}}} $ | (2) |
式(1)减去式(2)得到以下DFT递归关系
| $ {X_n}\left( m \right) = {X_{n - 1}}\left( m \right) + \left[ {x\left( n \right) - x\left( {n - N} \right)} \right]{{\rm{e}}^{ - {\rm{j}}\frac{{{\rm{2 \mathsf{ π} }}mn}}{N}}} $ | (3) |
式(3)为RDFT，用于估计电压波形中的频谱分量。对Xn(m)进行RDFT逆变换(IRDFT)，得到时域电压信号频率为mΔf分量的瞬时值
| $ {x_m}\left( n \right) = \frac{1}{N}{X_n}\left( m \right){{\rm{e}}^{{\rm{j}}\frac{{{\rm{2 \mathsf{ π} }}mn}}{N}}} $ | (4) |
通过式(3)和式(4)给出的RDFT和IRDFT，可以推导中心频率为mΔf的带通滤波器(BF)的离散传递函数。
首先，将式(3)代入式(4)得到时域信号的递推公式
| $ {x_m}\left( n \right) = \frac{1}{N}{X_{n - 1}}\left( m \right){{\rm{e}}^{{\rm{j}}\frac{{{\rm{2 \mathsf{ π} }}mn}}{N}}} + \frac{1}{N}\left[ {x\left( n \right) - x\left( {n - N} \right)} \right] $ | (5) |
其次，利用式(4)，将第(n－1)T采样时刻的时域电网电压信号表示为
| $ {x_m}\left( {n - 1} \right) = \frac{1}{N}{X_{n - 1}}\left( m \right){{\rm{e}}^{{\rm{j}}\frac{{{\rm{2 \mathsf{ π} }}m\left( {n - 1} \right)}}{N}}} $ | (6) |
然后将式(6)带入式(5)，得到以下递归关系
| $ {x_m}\left( n \right) = {x_m}\left( {n - 1} \right){{\rm{e}}^{{\rm{j}}\frac{{{\rm{2 \mathsf{ π} }}m}}{N}}} + \frac{1}{N}\left[ {x\left( n \right) - x\left( {n - N} \right)} \right] $ | (7) |
最后，对式(7)双边作z变换，得到BF离散传递函数为
| $ {H_m}\left( z \right) = \frac{{{X_m}\left( z \right)}}{{X\left( z \right)}} = \frac{1}{N}\frac{{1 - {z^{ - N}}}}{{1 - {z^{ - 1}}{{\rm{e}}^{{\rm{j2 \mathsf{ π} }}m/N}}}},m = 1,2 \cdots N $ | (8) |
式(8)中，m=1时，传递函数H1(z)的波特图如图 1所示，H1(z)的中心频率Δf为电网电压基频50 Hz。由图 1看出，H1(z)有别于低通滤波，可以滤除各次偶次和奇次谐波，还包括直流分量。
| 图 1 传递函数H1(z) (m=1)波特图 Fig.1 Bode plot of the transfer function H1(z) (m=1) |
在被测信号频率波动情况下，要保持截断窗口与基波周期同步，理想截断窗口N值可能为整数或非整数。当N为整数值时，图 1中H1(z)可以有效滤除谐波；当N为非整数值时，L < N < L+1，L为正整数，由于离散采样限制，可采用向上(ceil)或向下(floor)的近似方法获得近似整数N值。由于函数的舍入误差，H1(z)不能有效滤除谐波干扰，因此本文通过DFT内插公式，设计接近分数延迟的FIR滤波器，以解决BF完整滤除谐波干扰的问题。
2 DFT插值非整数延迟滤波器输入与输出信号关系的推导为了解决N为非整数值时，H1(z)不能有效滤除谐波干扰的问题，以三角基底插值设计偶数或奇数采样点下的离散加权值，来有效逼近所需要得到的非整数延迟，以提高非同步采样下RDFT滤波输出准确度。
对于长度为偶数N的实数序列x(n)，令其DFT为X(k)，设ξ为整数插值因子，则内插补零长度为M= ξN时的插值DFT函数为[18]
| $ {X_{{\rm{zp}}}}\left( k \right) = \left\{ \begin{array}{l} \xi X\left( k \right)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;0 \le k \le \frac{N}{2} - 1\\ 0.5\xi X\left( {\frac{N}{2}} \right)\;\;\;\;\;\;\;\;\;\;\;\;k = \frac{N}{2}\\ 0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{N}{2} + 1 \le k \le M - \frac{N}{2}\\ 0.5\xi X\left( {\frac{N}{2}} \right)\;\;\;\;\;\;\;\;\;\;\;\;k \le M - \frac{N}{2}\\ \xi X\left( {k - M + N} \right)\;\;\;\;\;\;M - \frac{N}{2} + 1 \le k \le M \end{array} \right. $ | (9) |
对式(9)中长度为M的Xzp(k)进行DFT逆变换，可得内插序列xzp(n)为
| $ \begin{array}{l} {x_{{\rm{zp}}}}\left( n \right) = \frac{1}{M}\sum\limits_{k = 0}^{M - 1} {{X_{{\rm{zp}}}}\left( k \right)W_M^{ - kn}} = \frac{1}{N}\left( {\sum\limits_{k = 0}^{N/2 - 1} {{X_{{\rm{zp}}}}\left( k \right)} } \right.\\ \left. {W_M^{ - kn} + \sum\limits_{k = M - N/2 + 1}^{M - 1} {{X_{{\rm{zp}}}}\left( k \right)W_M^{ - kn}} } \right) + \frac{1}{{2N}}\left( {{X_{{\rm{zp}}}}\frac{N}{2}W_M^{ - \left( {\frac{N}{2}} \right)n} + } \right.\\ \left. {{X_{{\rm{zp}}}}\frac{N}{2}W_M^{ - \left( {M - \frac{N}{2}} \right)n}} \right) \end{array} $ | (10) |
式(10)中，WM=e－j2π/M。令k′=M－k，且WM－Mn=1，可证得式(9)满足共轭对称Xzp(M－k)=Xzp(M－k)*, 1≤k≤N/2。此时，式(10)可化简为
| $ \begin{array}{l} {x_{{\rm{zp}}}}\left( n \right) = \frac{1}{N}\sum\limits_{m = 0}^{N - 1} {x\left( m \right)\left\{ {1 + 2\sum\limits_{k = 1}^{N/2 - 1} {\cos \left( {\frac{{2{\rm{ \mathsf{ π} }}nk}}{M} - } \right.} } \right.} \\ \left. {\left. {\frac{{2{\rm{ \mathsf{ π} }}mk}}{N}} \right) + {{\left( { - 1} \right)}^m}\cos \left( {\frac{{{\rm{ \mathsf{ π} }}Nk}}{M}} \right)} \right\} \end{array} $ | (11) |
由于补零长度满足条件M= ξN，故式(11)满足
| $ {x_{{\rm{zp}}}}\left( {n\xi } \right) = {x_{{\rm{zp}}}}\left( n \right) $ | (12) |
进一步，由于序列xzp(n)是序列x(n)以ξ为插值因子的内插样本，所以在满足0≤τ≤ ξ－1，0≤n≤N－1的条件下，根据三角变换关系可得近似关系如式(13)
| $ \begin{array}{l} {x_{{\rm{zp}}}}\left( {n + \frac{\tau }{\xi }} \right) \approx {x_{{\rm{zp}}}}\left( {n\xi + \tau } \right) = \frac{1}{N}\sum\limits_{m = 0}^{N - 1} {x\left( m \right)\left\{ {1 + } \right.} \\ 2\sum\limits_{k = 1}^{N/2 - 1} {\cos \left( {\frac{{2{\rm{ \mathsf{ π} }}\left( {n + \frac{\tau }{\xi } - m} \right)k}}{N}} \right) + {{\left( { - 1} \right)}^m}\cos \left( {{\rm{ \mathsf{ π} }}\left( {n + } \right.} \right.} \\ \left. {\left. {\left. {\frac{\tau }{\xi }} \right)} \right)} \right\} = \sum\limits_{m = 0}^{N - 1} {x\left( m \right)\left\{ {\frac{1}{N} + \sum\limits_{k = 0}^{N/2} {{b_k}} } \right.} \\ \left. {\cos \left( {\frac{{2{\rm{ \mathsf{ π} }}\left( {n + \tau /\xi - m} \right)k}}{N}} \right)} \right\} \end{array} $ | (13) |
式(13)可用插值形式表示为
| $ {x_{{\rm{zp}}}}\left( l \right) \approx \sum\limits_{m = 0}^{N - 1} {x\left( m \right){B_{\rm{e}}}\left( {m,l} \right)} $ | (14) |
式(14)中，l=n+τ/ ξ，Be(m, l)表示偶数插值基函数，Be(m, l)的权重系数为
| $ {b_k} = \left\{ \begin{array}{l} 1\;\;\;\;k = 0,\frac{N}{2}\\ 2\;\;\;k = 1,2, \cdots ,\frac{N}{2} - 1 \end{array} \right. $ | (15) |
图 2为N=10且l取不同值时偶数插值基函数Be(m, l)曲线。从图 2可以明显看出
| 图 2 m取不同值时的偶数插值基函数Be(m, l)曲线 Fig.2 Basis function of even interpolation Be(m, l) for different values of m |
| $ {B_{\rm{e}}}\left( {m,l} \right) = \left\{ \begin{array}{l} 1\;\;\;\;l = m:m \in Z\;整数集合\\ 0\;\;\;\;l \ne m,l = 0,1, \cdots ,N - 1 \end{array} \right. $ | (16) |
式(16)是DFT插值对偶数插值奇函数Be(m, l)的基本要求。类似地，可得对应奇数插值基函数为
| $ {B_{\rm{o}}}\left( {m,l} \right) = \frac{1}{{N + 1}}\sum\limits_{k = 0}^{\left( {N + 1} \right)/2} {{c_k}\cos \left( {\frac{{2{\rm{ \mathsf{ π} }}\left( {n + \frac{\tau }{\xi } - m} \right)k}}{{N + 1}}} \right)} $ | (17) |
式(17)中，权重系数ck为
| $ {c_k} = \left\{ \begin{array}{l} 1\;\;\;\;k = 0,\frac{{N + 1}}{2}\\ 2\;\;\;k = 1,2, \cdots ,\frac{{N + 1}}{2} - 1 \end{array} \right. $ | (18) |
式(14)给出了DFT插值非整数延迟滤波器输入与输出信号之间关系，据此关系可以设计DFT插值非整数延迟滤波器的传递函数。
3 DFT插值非整数延迟滤波的传递函数设计式(8)中，基波频率偏移使得准同步采样点数N为非整数时，取z-N=z-(floor(N)+d)，d∈[0, 1)。理想情况下，当信号ν(n)经过非整数延迟z-N滤波后，输出为y(n)=ν(n－floor(N)－d)；非理想准同步采样情况下，输出为y(n)=ν(n－floor(N) or ceil(N))，此时由于函数的舍入误差，滤波器不能达到有效滤除谐波干扰的效果。因此在离散采样限制的情况下，非整数延迟滤波器的冲击响应h(r)与延迟输出满足以下关系[18]
| $ v\left( {n - L - d} \right) = \sum\limits_{k = 0}^{N - 1} {h\left( k \right)v\left( {n - k} \right)} $ | (19) |
式中，L=floor(N)。为了推导非整数延迟滤波器的传递函数，在式(14)中，令xzp(l)=v(n－(N－1)+l)，得
| $ v\left( {n - \left( {N - 1} \right) + l} \right) \approx \sum\limits_{k = 0}^{N - 1} {v\left( {n - \left( {N - 1} \right) + k} \right)}B(k,l) $ | (20) |
进一步取l=N-1-L-d，k=N-1-r，式(20)可变换为
| $ v\left( {n - L - d} \right) \approx \sum\limits_{r = 0}^{N - 1} {v\left( {n - r} \right)B\left( {N - 1 - r,N - 1 - L - d} \right)} $ | (21) |
对比式(19)和式(21)，在偶数点DFT插值时，可得出整数滤波器离散冲击响应he(r)为
| $ \begin{array}{l} {h_{\rm{e}}}\left( r \right) = B\left( {N - 1 - r,N - 1 - L - d} \right) = \\ \frac{1}{N}\sum\limits_{k = 0}^{N/2} {{b_k}\cos \left( {\frac{{2{\rm{ \mathsf{ π} }}\left( {r - L - d} \right)k}}{N}} \right)} \end{array} $ | (22) |
同理，在奇数点DFT插值时，可得出整数滤波器离散冲击响应ho(r)为
| $ \begin{array}{l} {h_{\rm{o}}}\left( r \right) = B\left( {N - 1 - r,N - 1 - L - d} \right) = \\ \frac{1}{{N + 1}}\sum\limits_{k = 0}^{\left( {N + 1} \right)/2} {{b_k}\cos \left( {\frac{{2{\rm{ \mathsf{ π} }}\left( {r - L - d} \right)k}}{{N + 1}}} \right)} \end{array} $ | (23) |
对于冲击响应h(r)，其离散z域DFT插值的非整数延迟滤波器的传递函数为
| $ {H_{{\rm{DFTI}}}} = {z^{ - N}} = \sum\limits_{r = 0}^{K - 1} {h\left( r \right){z^{ - r}}} $ | (24) |
因此，在离散整数采样下，将式(24)代入式(8)，简化后DFT插值非整数延迟递归滤波器HRDFTI的传递函数为
| $ {H_{{\rm{RDFTI}}}}\left( z \right) = \frac{1}{{L + d}} \times \frac{{1 - \sum\limits_{r = 0}^{K - 1} {h\left( r \right){z^{ - r}}} }}{{1 - {z^{ - 1}}{{\rm{e}}^{{\rm{j2 \mathsf{ π} }}/\left( {L + d} \right)}}}} $ | (25) |
在非同步采样情况下，为了提高式(8)中H1(z)的滤波性能，在信号采样点v(n-L)和v(n-L-1)之间进行线性插值，则对应离散z域线性插值滤波器的传递函数为
| $ {H_{{\rm{LI}}}} = {z^{ - N}} = \left( {1 - d} \right){z^{ - L}} + d{z^{ - \left( {L + 1} \right)}} $ | (26) |
式中，线性插值权值因子d=N－L，L=floor(N)。设基波频率fd=49.8 Hz，采样频率fs=5 000 Hz，理想采样点数N=5 000 Hz/49.8 Hz≈100.4，则式(26)中L=100，d=0.4。此时，式(26)传递函数可表示为
| $ {H_{{\rm{LI}}}} = \left( {1 - 0.4} \right){z^{ - L}} + 0.4{z^{ - L - 1}} $ | (27) |
式(24)中h(r)取相同参数L=100，d=0.4时，分析得到式(24)、(26)对应的幅度-频率特性曲线，如图 3所示。
| 图 3 两种插值方法的幅频特性曲线 Fig.3 Frequency response of the two interpolation methods |
由图 3可知，DFT插值非整数延迟滤波器HDFTI在0~1归一化宽频带范围内与理想滤波器Hfd幅频特性曲线吻合；而线性插值滤波器HLI仅在0~0.2频带范围内与Hfd基本吻合，在0.2~1频带范围内与Hfd幅频特性偏离较大。所以HDFTI滤波器明显优于HLI。
4.1.2 滤波性能为了分析非整数延迟滤波器的滤波性能，定义归一化均方根误差
| $ {E_{{\rm{NRMS}}}} = \sqrt {\left( {\frac{{\int_0^{\varepsilon {\rm{ \mathsf{ π} }}} {{{\left| {{H_{{\rm{in}}}} - {H_{{\rm{fd}}}}} \right|}^2}{\rm{d}}\omega } }}{{\int_0^{\varepsilon {\rm{ \mathsf{ π} }}} {{{\left| {{H_{{\rm{fd}}}}} \right|}^2}{\rm{d}}\omega } }}} \right)} $ | (28) |
式(28)中，Hfd为理想滤波器传递函数，Hin为被估计滤波器传递函数，ENRMS越小意味着滤波器的性能越好。
在归一化频带0~0.2范围内，利用式(28)分析HDFTI和HLI两种滤波器的归一化均方根误差。电网电压频率波动0.2 Hz、延迟量d∈[0, 0.5]时，得到HDFTI和HLI的归一化均方根误差如图 4所示。图 4表明，HDFTI的归一化均方根误差比HLI小3个数量级，说明HDFTI性能明显优于HLI。
| 图 4 不同延迟参数下两种插值方法归一化均方根误差 Fig.4 Normalized root mean square errors of two kinds of interpolation methods for different delay parameters |
将线性插值和DFT插值方法的非整数延迟滤波器h(r)与式(8)带通滤波器联合使用分析信号基波，则式(8)带通滤波器可改写为
| $ {H_1}\left( z \right) = \frac{1}{{L + d}} \times \frac{{1 - {z^{ - \left( {L + d} \right)}}}}{{1 - {z^{ - 1}}{{\rm{e}}^{{\rm{j2 \mathsf{ π} }}/\left( {L + d} \right)}}}} $ | (29) |
设采样频率fs=5 000 Hz，则滤波输入的多频电压信号为
| $ \begin{array}{l} {V_{{\rm{sig}}}}\left( t \right) = \cos \left( {2{\rm{ \mathsf{ π} }}{f_{\rm{d}}}t} \right) + 0.1\cos \left( {2{\rm{ \mathsf{ π} 2}}{f_{\rm{d}}}t} \right) + \\ 0.5\cos \left( {2{\rm{ \mathsf{ π} 3}}{f_{\rm{d}}}t} \right) \end{array} $ | (30) |
对于工频50 Hz基波信号，单周期采样点数为100时，式(29)中对应参数L=100，d=0。对Vsig(t)进行滤波，得到滤波前输入的谐波信号Vsig(t)和基波信号波形分别如图 5(a)和图 5(b)所示，滤波后输出的基波波形和基波输出绝对误差分别如图 5(c)和图 5(d)所示。由图 5(d)可以看出，绝对误差过渡带的时间非常短，滤波器响应延迟仅为一个工频周期；在过渡带以后，滤波器性能非常优越。这说明式(29)所示的带通滤波器误差小、响应快，非常适用于电网信号快速波动情况下的谐波分析。
| 图 5 多频电压信号的滤波及误差 Fig.5 Filtering and errors of the multi-frequency voltage signal |
为了进一步定量评估DFT插值非整数延迟滤波器的谐波分析误差，保持fs不变，将基波频率偏离到49.8 Hz时的非理想状态，引入谐波分析误差评估公式
| $ {E_{{\rm{RMS}}}} = \sqrt {\left( {\int_{t = {t_1}}^{{t_{{\rm{end}}}}} {{{\left| {\left( {{v_{{\rm{fo}}}} - {v_{{\rm{fd}}}}} \right)} \right|}^2}{\rm{d}}t} } \right)} $ | (31) |
式中，vfo为滤波输出，vfd为理想滤波输出。滤波器输出仅考虑基波k=1，其他谐波与基波分析类似，分析中舍去了滤波器输出起始时刻一个工频周期的响应延迟。
对应电压基波频率49.8 Hz，取延迟量d=5 000/49.8-100≈0.4，DFT插值长度L=100，则线性插值递归滤波器HRLI的传递函数为
| $ {H_{{\rm{RLI}}}}\left( z \right) = \frac{1}{{L + d}} \times \frac{{1 - \left( {L - d} \right){z^{ - L}} - d{z^{ - \left( {L + 1} \right)}}}}{{1 - {z^{ - 1}}{{\rm{e}}^{{\rm{j2 \mathsf{ π} }}/\left( {L + d} \right)}}}} $ | (32) |
根据误差评估公式(31)计算得到3种递归滤波器输出误差如表 1所示。由表 1可知，DFT插值非整数延迟递归滤波器HRDFTI的滤波准确度比理想递归滤波器HRfd和线性插值递归滤波器HRLI分别提高66.2%和23.3%。
| 下载CSV 表 1 不同滤波输出误差 Table 1 Output errors of different filters |
通过DFT插值非整数延迟滤波器对谐波进行滤波后，滤波器输出仅含有基波成分，因此采用幅度归一化后的基波电压估计频率可进一步提高频率估计准确度。根据式(3)，可得到n时刻基波幅值
| $ {V_{{\rm{aof}}}}\left( n \right) = \frac{2}{N}\sqrt {{{\left[ {{\mathop{\rm Re}\nolimits} \left( {{V_1}\left( n \right)} \right)} \right]}^2} + {{\left[ {{\mathop{\rm Im}\nolimits} \left( {{V_1}\left( n \right)} \right)} \right]}^2}} $ | (33) |
其中
| $ \left\{ \begin{array}{l} {\mathop{\rm Re}\nolimits} \left( {{V_1}\left( n \right)} \right) = {\mathop{\rm Re}\nolimits} \left( {{V_1}\left( {n - 1} \right)} \right) + \\ \;\;\;\;\;\;\;\left[ {\nu \left( n \right) - \nu \left( {n - N} \right)} \right]\cos \left( {\frac{{2{\rm{ \mathsf{ π} }}n}}{N}} \right)\\ {\mathop{\rm Im}\nolimits} \left( {{V_1}\left( n \right)} \right) = {\mathop{\rm Im}\nolimits} \left( {{V_1}\left( {n - 1} \right)} \right) + \\ \;\;\;\;\;\;\;\left[ {\nu \left( {n - N} \right) - \nu \left( n \right)} \right]\sin \left( {\frac{{2{\rm{ \mathsf{ π} }}n}}{N}} \right) \end{array} \right. $ |
同样，可得归一化电压基波为
| $ v_1^{{\rm{naof}}} = \frac{{{v_1}\left( n \right)}}{{{V_{{\rm{aof}}}}\left( n \right)}} = \cos \left( {2{\rm{ \mathsf{ π} }}n{f_{\rm{d}}}/{f_{\rm{s}}} + {\phi _{\rm{v}}}} \right) $ | (34) |
式中ϕv为基波电压初始相位。
为了评估基波测量和频率估计准确度，定义基波测量和频率估计误差评估公式
| $ {E_{{\rm{SM}}}} = \frac{1}{R}\sum\limits_1^R {\frac{{\left| {{{\hat x}_{\rm{d}}} - {x_{\rm{d}}}} \right|}}{{{x_{\rm{d}}}}}} $ | (35) |
式(35)中，R为对应总的相位点数，xd为原始值(基波幅值或频率)，
频率估计采用加汉宁窗插值FFT的方法，基波频率fd=49.8 Hz，相位设置的范围为0°~180°，变化间隔为5°，干扰谐波电压次数为3、5、7、9次，相位设置为0°。取所有相位点下的基波电压测量和频率估计误差平均值，由公式(34)仿真得出基波测量和频率估计误差如表 2所示。
| 下载CSV 表 2 不同滤波下基波测量和频率估计误差 Table 2 The errors of different filters for fundamental wave measurement and frequency estimation |
表 2结果表明，本文所设计的DFT插值非整数延迟递归滤波器HRDFTI有效降低了谐波干扰对电压测量与频率估计的影响，基波测量和频率估计的误差分别为0.53%和0.91%，明显优于另外两种滤波器。
6 结论(1) 通过推导RDFT带通滤波器的传递函数，设计了一种非整数延迟滤波器，仿真表明，同步采样情况下理想带通滤波器不产生理论测量误差，而非同步采样情况下，由于基波频率偏移，非整数延迟滤波器的滤波准确度比理想滤波器提高66.2%，比线性插值滤波器提高23.3%。
(2) 利用DFT插值非整数延迟递归滤波器滤除谐波之后，基波电压测量和频率估计的准确度较理想滤波器分别提高60.5%和72.5%，较线性插值滤波器分别提高44.8%和41.7%，表明非整数延迟滤波器用于基波电压测量与频率估计可以有效降低谐波干扰。
| [1] |
CHAUDHARY S K, LASCU C, TEODORESCU R, et al. Voltage feedback based harmonic compensation for an offshore wind power plant[C]//2016 IEEE International Conference on Power Electronics, Drives and Energy Systems. Trivandrum, 2016: 1-5.
|
| [2] |
WANG W, RUAN X. A modified reference of an intermediate bus capacitor voltage-based second-harmonic current reduction method for a standalone photovoltaic power system[J]. IEEE Transactions on Power Electronics, 2016, 31(8): 5562-5573. DOI:10.1109/TPEL.2015.2497314 |
| [3] |
AL-FEILAT E A A, EL-AMIN I, BETTAYEB M. Power system harmonic estimation:a comparative study[J]. Electric Power Systems Research, 1994, 29(94): 91-97. |
| [4] |
DUDA K. DFT interpolation algorithm for Kaiser-Bessel and Dolph-Chebyshev windows[J]. IEEE Transactions on Instrumentation and Measurement, 2011, 60(3): 784-790. DOI:10.1109/TIM.2010.2046594 |
| [5] |
GRANDKE T. Interpolation algorithms for discrete Fourier transforms of weighted signals[J]. IEEE Transactions on Instrumentation & Measurement, 1983, 32(2): 350-355. |
| [6] |
牛胜锁, 梁志瑞, 张建华, 等. 基于四项余弦窗三谱线插值FFT的谐波检测方法[J]. 仪器仪表学报, 2012, 33(9): 2002-2008. NIU S S, LIANG Z R, ZHANG J H, et al. Harmonic detection approach based on 4-term cosine window triple-spectrum-line interpolation FFT[J]. Chinese Journal of Scientific Instrument, 2012, 33(9): 2002-2008. (in Chinese) |
| [7] |
许珉, 刘凌波. 基于三次样条函数的加Blackman-harris窗插值FFT算法[J]. 电力自动化设备, 2009, 29(2): 59-63. XU M, LIU L B. Blackman-harris window interpolated FFT algorithm based on cubic spline function[J]. Electric Power Automation Equipment, 2009, 29(2): 59-63. (in Chinese) |
| [8] |
LIN H C, LEE C S. Enhanced FFT-based parametric algorithm for simultaneous multiple harmonics analysis[J]. IEE Proceedings-Generation, Transmission and Distribution, 2001, 148(3): 209-214. DOI:10.1049/ip-gtd:20010278 |
| [9] |
VARKONYI-KOCZY A R. A recursive fast Fourier transformation algorithm[J]. IEEE Transactions on Circuits and Systems Ⅱ:Analog and Digital Signal Processing, 1995, 42(9): 614-616. DOI:10.1109/82.466641 |
| [10] |
ROSCOE A J, BURT G M, MCDONALD J R. Frequency and fundamental signal measurement algorithms for distributed control and protection applications[J]. IET Generation, Transmission & Distribution, 2009, 3(5): 485-495. |
| [11] |
DARWISH H A, FIKRI M. Practical considerations for recursive DFT implementation in numerical relays[J]. IEEE Transactions on Power Delivery, 2007, 22(1): 42-49. DOI:10.1109/TPWRD.2006.874642 |
| [12] |
FUNAKI T, MATSUURA K, TANAKA S. Error correction for phase detection by recursive algorithm real time DFT[J]. Electrical Engineering in Japan, 2002, 141(1): 8-17. DOI:10.1002/(ISSN)1520-6416 |
| [13] |
BORKOWSKI J, KANIA D, MROCZKA J. Interpolated-DFT-based fast and accurate frequency estimation for the control of power[J]. IEEE Transactions on Industrial Electronics, 2014, 61(12): 7026-7034. DOI:10.1109/TIE.2014.2316225 |
| [14] |
高云鹏, 李峰, 陈婧, 等. 改进能量算子Nuttall窗频谱校正的闪变测量方法[J]. 仪器仪表学报, 2014, 35(4): 745-753. GAO Y P, LI F, CHEN J, et al. Voltage flicker measurement method using improved Teager-Kaiser energy operator based on Nuttall window spectral correction[J]. Chinese Journal of Scientific Instrument, 2014, 35(4): 745-753. (in Chinese) |
| [15] |
MITRA S K. Digital signal processing:a computer-based approach[M]. 4th ed. New York: McGraw Hill Higher Education, 2011.
|
| [16] |
SMITH S. Digital signal processing:a practical guide for engineers and scientists[M]. 3rd ed. Boston: Newnes, 2002.
|
| [17] |
PEI S C, WANG P H, LIN C H. Design of fractional delay filter, differintegrator, fractional Hilbert transformer, and differentiator in time domain with Peano kernel[J]. IEEE Transactions on Circuits and Systems Ⅰ:Regular Papers, 2010, 57(2): 391-404. DOI:10.1109/TCSI.2009.2027806 |
| [18] |
OPPENHEIM A V, SCHAFER R W. Discrete-time signal processing[M]. 3rd ed. New York: Pearson Education Limited, 2013.
|