FFT Nuttall窗四谱线插值修正谐波检测算法

Modified Harmonic Detection Algorithm by Nuttall Window Four⁃Line Interpolation for FFT

摘要

利用传统傅里叶变换进行谐波信号分析时， 由于很难严格做到同步采样， 整周期截断， 从而不可避免地产生频谱泄露与栏栅效应， 降低了检测结果的精确性。为了提高检测结果精确度， 在FFT（Fast Fourier Transformation）的基础上， 对Hann窗、 Blackman窗、 Blackmanharris窗、 Nuttall窗的四谱线插值算法进行了研究， 通过多项式拟合推导出不同余旋窗四线插值对应的幅值、 频率、 相位的修正公式。并在Matlab上对不同余旋窗四线插值算法进行仿真实验， 与对应余旋窗的三谱线插值算法进行对比分析， 结果表明四项三阶Nuttall窗的四谱线插值算法幅值精确度可达到10^-5~10^-9数量级， 频率精确度可达到 10^-8~10^-9数量级， 相位精确度可达到10^-3~10^-7数量级， 与其它算法相比精确度更高， 旁瓣性能更好， 能更好地抑制频谱泄露与栏栅效应。

Abstract

When the traditional Fourier transform is used for harmonic signal analysis， it is difficult to strictly synchronize the sampling and truncate the whole period， which inevitably produces spectrum leakage and fence effect， and reduces the accuracy of the detection results. To improve the accuracy of detection results based on FFT， four-line interpolation algorithms of Hann window， Blackman window， Blackmanharris window， and Nuttall window are studied. The correction formulas of amplitude， frequency and phase corresponding to the four-line interpolation of different corotation windows are derived by polynomial fitting. In addition， the simulation experiments of different four-line interpolation algorithms for corotating windows and the corresponding three-line interpolation algorithms for corotating windows are compared and analyzed in Matlab. The simulation results show that the amplitude accuracy of the four-line interpolation algorithms for four third-order Nuttall windows can reach 10^-5~10^-9 orders， and the frequency accuracy can reach 10^-8~10^-9 orders. The phase accuracy can reach 10^-3~10^-7 orders of magnitude. Compared with other algorithms， the accuracy is higher， the sidelobe performance is better， and the spectrum leakage and fence effect can be better suppressed.

Graphical abstract

关键词

频谱泄露 / 栏栅效应 / 窗函数 / 三线谱插值 / 四线谱插值

Key words

spectrum leakage / fence effect / window function / three-line spectral interpolation / four-line spectrum interpolation

0 引　言

随着国家经济的发展， 电力电子技术逐渐成熟， 许多非线性电力电子器件越来越频繁地出现在生活中， 造成电网谐波问题愈发严重。谐波问题不但会影响电能质量及运行效率， 还可能会影响电网的运行， 造成安全隐患^［1‐2］。因此， 信号的谐波治理尤为重要。对信号谐波治理的前提是对信号谐波进行相位、 频率、 幅值的精确测量。迄今为止， 谐波分析使用最广泛的方法是傅里叶变换， 其要求整周期截断与同步采样， 但在实际分析时很难达到此条件， 会造成频谱泄露与栏栅效应， 对谐波参数检测的精确度造成很大影响^［3］。

对于傅里叶变换存在的关于频谱泄露与栏栅效应的问题主要用窗函数算法^［4］与谱线插值算法去修正快速傅里叶变换（Fast Fourier Transformation， FFT）变换对谐波检测精度造成的影响。在待测谐波信号上加性能良好的窗函数可以有效降低非正周期截断在频域造成的频谱泄露。使用谱线插值修正可以利用多条谱线加权计算， 降低栏栅效应引起的误差， 进而提高FFT对谐波信号检测的精度。加窗算法常用的窗函数有三角窗、 矩形窗、 Hamming窗、 Hanning窗等。插值算法中， 双谱线插值算法、 三谱线插值算法等使用较多。周俊等^［5］于2006年提出基于Blackman窗函数的插值FFT算法， 仿真实验表明该算法对谐波参数的检测精度都有所提高。周西峰等^［6］于2014年提出的基于Blackman-Harris窗插值分析法对频谱的泄露抑制有一定的效果。但是因为这些算法所使用的窗函数旁瓣性能不够理想， 进而造成对信号频谱泄漏的抑制能力有限。卿柏元等^［7］提出基于Nuttall窗双谱线插值修正算法， 在一定程度上提高了计算精度， 但双谱线的插值修正只考虑最高峰值谱线附近的两条次高谱线， 以此求取的偏移量去修正各个谐波参数， 并没有考虑峰值频点附近的泄露谱线信息。李威等^［8］提出基于Nuttall窗的三谱线插值修正算法， 并进行仿真实验证明相比于常见窗函数双峰谱线的插值修正算法精确度更高， 然而三谱线插值修正算法对待测谐波峰值频点两侧对称谱线所包含的信息没能进行充分利用。为提高检测精度， 本文研究了几种常见窗的四线插值修正算法， 并与不同窗函数的三谱线插值修正算法进行分析对比， 结果表明旁瓣性能良好的Nuttall窗的四谱线插值算法对谐波信号的检测精确度有极大地提高。

1 传统傅里叶变换存在的问题

利用傅里叶变换对信号进行分析时， 理论上信号为无限长序列， 但在实际测量中， 被测信号都是截取的有限长序列^［9］。对于长度有限的信号序列进行快速傅里叶变换分析时， 需要将有限的信号序列进行周期延拓至整个时域上^［10］。这种情况下， 只有对信号进行同步采样才能确保周期延拓后的序列信号与原来待测的序列信号一致， 才能保证信号谐波检测分析结果的精确性。在工程实际测量中很难做到真正意义上的同步采样， 这时所用的时间窗长度同信号周期不再一致导致采样周期延拓首末两端不连续， 产生跳变， 发生频谱泄露^［11］。除此之外， 非同步采样还会带来的另一问题是栏栅效应。信号采样分析需要将连续信号离散化， 离散谱线像栅栏一样， 如果谱线正好在采样点之间， 那就无法被识别， 对应分量的频率、 幅值、 相位都无法准确获取， 这就是栅栏效应^［12］。

2 经典余旋窗

为了解决非同步采样带来的频谱泄露问题， 采用加窗函数抑制， 在加窗函数的选择上应尽可能在截断的地方选择平缓的， 以保留更多的信息。不同特性的窗函数对采样精度有很大影响， 旁瓣衰减速度小的余旋窗函数主瓣通常比较窄， 频谱泄露主要发生在旁瓣范围内。而主瓣比较宽的余旋窗函数通常旁瓣衰减速度大， 频谱泄露集中在主瓣范围内^［13］。通常情况下选择窗函数的原则是：

1） 旁瓣衰减速度快而且旁瓣峰值足够小的余旋窗函数抑制信号频谱泄露的效果更好；

2） 主瓣宽度窄的窗函数能提高被淹没的谱线频率的分辨率^［14］。

在选择余旋窗函数时， 一般情况下无法同时考虑窗函数的主瓣特性和旁瓣特性， 可以根据情况恰当地用主瓣宽度来降低旁瓣峰值高度， 达到抑制频谱泄露的目的。

常用的窗函数在离散时域和频域的表达式如下：

1）　Hann窗

设Hann窗长度为N， 时域表达式为

w(n)=\mathrm{0.5}-\mathrm{0.5}\mathrm{c}\mathrm{o}\mathrm{s}\frac{\mathrm{2}\mathrm{\pi }n}{N}\mathrm{。}（1）

频域表达式为

W(\omega )=\frac{N\omega \mathrm{s}\mathrm{i}\mathrm{n}(\omega \mathrm{\pi })}{\mathrm{\pi }}{\mathrm{e}}^{-\mathrm{j}\omega \mathrm{\pi }}\left(\frac{\mathrm{0.5}}{{\omega }^{\mathrm{2}}}-\frac{\mathrm{0.5}}{{\omega }^{\mathrm{2}}-\mathrm{1}}\right)\mathrm{。}（2）

2）　Blackman窗

时域表达式为

w(n)=\mathrm{0.42}-\mathrm{0.5}\mathrm{c}\mathrm{o}\mathrm{s}\frac{\mathrm{2}\mathrm{\pi }n}{N}+\mathrm{0.08}\mathrm{c}\mathrm{o}\mathrm{s}\frac{\mathrm{4}\mathrm{\pi }n}{N}\mathrm{。}（3）

频域表达式为

W(\omega )=\frac{N\omega \mathrm{s}\mathrm{i}\mathrm{n}(\omega \mathrm{\pi })}{\mathrm{\pi }}{\mathrm{e}}^{-\mathrm{j}\omega \mathrm{\pi }}\left(\frac{\mathrm{0.42}}{{\omega }^{\mathrm{2}}}-\frac{\mathrm{0.5}}{{\omega }^{\mathrm{2}}-\mathrm{1}}+\frac{\mathrm{0.08}}{{\omega }^{\mathrm{2}}-\mathrm{4}}\right)\mathrm{。}（4）

3）　Blackmanharris窗

时域表达式为

\begin{array}{c}w(n)=\mathrm{0.358}\mathrm{75}-\mathrm{0.488}\mathrm{29}\mathrm{c}\mathrm{o}\mathrm{s}\frac{\mathrm{2}\mathrm{\pi }n}{N}+\\ \mathrm{0.141}\mathrm{28}\mathrm{c}\mathrm{o}\mathrm{s}\frac{\mathrm{4}\mathrm{\pi }n}{N}-\mathrm{0.011}\mathrm{68}\mathrm{c}\mathrm{o}\mathrm{s}\frac{\mathrm{6}\mathrm{\pi }n}{N}\mathrm{。}\end{array}（5）

频域表达式为

\begin{array}{c}W(\omega )=\frac{N\omega \mathrm{s}\mathrm{i}\mathrm{n}(\omega \mathrm{\pi })}{\mathrm{\pi }}·\\ {\mathrm{e}}^{-\mathrm{j}\omega \mathrm{\pi }}\left(\frac{\mathrm{0.358}\mathrm{75}}{{\omega }^{\mathrm{2}}}-\frac{\mathrm{0.488}\mathrm{29}}{{\omega }^{\mathrm{2}}-\mathrm{1}}+\frac{\mathrm{0.141}\mathrm{28}}{{\omega }^{\mathrm{2}}-\mathrm{4}}-\frac{\mathrm{0.011}\mathrm{68}}{{\omega }^{\mathrm{2}}-\mathrm{9}}\right)\mathrm{。}\end{array}（6）

4）　Nuttall窗

时域表达式为

\begin{array}{c}w(n)=\mathrm{0.312}\mathrm{5}-\mathrm{0.468}\mathrm{75}\mathrm{c}\mathrm{o}\mathrm{s}\frac{\mathrm{2}\mathrm{\pi }n}{N}+\\ \mathrm{0.187}\mathrm{5}\mathrm{c}\mathrm{o}\mathrm{s}\frac{\mathrm{4}\mathrm{\pi }n}{N}-\mathrm{0.031}\mathrm{25}\mathrm{c}\mathrm{o}\mathrm{s}\frac{\mathrm{6}\mathrm{\pi }n}{N}\mathrm{。}\end{array}（7）

频域表达式为

\begin{array}{c}W(\omega )=\frac{N\omega \mathrm{s}\mathrm{i}\mathrm{n}(\omega \mathrm{\pi })}{\mathrm{\pi }}·\\ {\mathrm{e}}^{-\mathrm{j}\omega \mathrm{\pi }}\left(\frac{\mathrm{0.312}\mathrm{5}}{{\omega }^{\mathrm{2}}}-\frac{\mathrm{0.468}\mathrm{75}}{{\omega }^{\mathrm{2}}-\mathrm{1}}+\frac{\mathrm{0.187}\mathrm{5}}{{\omega }^{\mathrm{2}}-\mathrm{4}}-\frac{\mathrm{0.031}\mathrm{25}}{{\omega }^{\mathrm{2}}-\mathrm{9}}\right)\mathrm{。}\end{array}（8）

上述时域表达式中N为采样点数， n=0， 1，2，…，N-1； 频域表达式中ω为角频率。

3 谱线插值算法原理

非同步采样带来的的栅栏效应是不可避免的， 谱线插值算法通过待测峰值谱线两端幅值最高的谱线去估计真实的谱线幅值。利用数据拟合逼近， 根据相应的修正公式求解谱线参数， 抑制栅栏效应， 减小对谐波信号检测结果的影响， 从而提高谐波信号检测的精度^［15］。本文主要对窗函数四谱线插值修正算法与对应窗函数三谱线插值算法做分析对比研究。

3.1 三谱线插值

对输入的谐波信号x(t)用频率为f_s采样离散处理后可以表示为

x(n)={\displaystyle \sum _{h=\mathrm{1}}^H}{A}_h\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{2}\mathrm{\pi }\frac{f_h}{f_s}n+{\varphi }_h\right),（9）

式中： h为信号谐波次数； n为信号采样点数； A_h为信号第h次谐波的幅值； {\varphi }_h为信号第h次谐波的相位； f_h为第h次谐波的频率。用窗函数对信号截断处理后的离散傅里叶变换为

X（k）={\displaystyle \sum _{h=\mathrm{1}}^H}\frac{A_h}{\mathrm{2}\mathrm{j}}\left[{\mathrm{e}}^{\mathrm{j}{\varphi }_h}W\left(k-\frac{f_h}{\mathrm{\Delta }f}\right)-{\mathrm{e}}^{\mathrm{j}{\varphi }_h}{W}_{}\left(k+\frac{f_h}{\mathrm{\Delta }f}\right)\right],

k=\mathrm{0,1},\mathrm{2},\cdots ,N-\mathrm{1},（10）

式中： \mathrm{\Delta }f为采样间隔， 忽略负频点的旁瓣影响， 有

X（k）={\displaystyle \sum _{h=\mathrm{1}}^H}\frac{A_h}{\mathrm{2}\mathrm{j}}\left[{\mathrm{e}}^{\mathrm{j}{\varphi }_h}W\left(k-\frac{f_h}{\mathrm{\Delta }f}\right)\right],

k=\mathrm{0,1},\mathrm{2},\cdots ,N-\mathrm{1}\mathrm{。}（11）

在实际工程中， 因为不可避免会发生非同步采样， 离散频谱点与实际的峰值频点有相应的偏移量。如图 1 所示， 设峰值频点为k， k附近的最大峰值谱线为k₁， 其左右两边的谱线分别为k₀， k₂。3条谱线对应的幅值分别为： y_{\mathrm{0}}=\left|X(k_{\mathrm{0}})\right|， y_{\mathrm{1}}=\left|X(k_{\mathrm{1}})\right|， y_{\mathrm{2}}=\left|X(k_{\mathrm{2}})\right|。由于0≤k_{\mathrm{1}}-k_{\mathrm{01}}≤1， 令\alpha =k-k_{\mathrm{1}}-\mathrm{0.5}， 即\alpha \in [-\mathrm{0.5,0.5}]， \beta =\frac{(y_{\mathrm{2}}-y_{\mathrm{0}})}{y_{\mathrm{1}}}。

由式（11）得

\beta =\frac{\left|W（\mathrm{1}-\alpha ）\right|-W\left|-\mathrm{1}-\alpha \right|}{\left|W（-\alpha ）\right|}\mathrm{。}（12）

其反函数为\alpha =g^{-\mathrm{1}}(\beta )， 通过多项式拟合方法求得修正公式， 由\beta可求出参数\alpha。频率与相位的修正公式^［16］分别为

f_h=(k_{\mathrm{1}}+\alpha +\mathrm{0.5})\mathrm{\Delta }f,

{\varphi }_h=\mathrm{a}\mathrm{r}\mathrm{g}[y(k_{\mathrm{1}})]+\frac{\mathrm{\pi }}{\mathrm{2}}-\mathrm{\pi }(\alpha +\mathrm{0.5})\mathrm{。}（13）

谐波信号的幅值修正需要对3条谱线进行加权平均计算， k_{\mathrm{1}}为距离实际峰值点最近的最大谱线， 在计算式中应该给予谱线k_{\mathrm{1}}的幅值y_{\mathrm{1}}较大的权重^{［17‐19］}， 由式（3）得

\begin{array}{c}\mathrm{2}{y}_{\mathrm{1}}+y_{\mathrm{0}}+y_{\mathrm{2}}=\mathrm{0.5}{A}_{\mathrm{0}}（\left|W（-\mathrm{1}-\alpha ）\right|+\\ \mathrm{2}\left|W（-\alpha ）\right|+\left|W（\mathrm{1}-\alpha ）\right|）\mathrm{。}\end{array}（14）

由式（11）得

A_h=\frac{（\mathrm{2}{y}_{\mathrm{1}}+{\mathrm{y}}_{\mathrm{0}}+y_{\mathrm{2}}）}{\mathrm{2}\left|W（-\alpha ）\right|+\left|W（-\mathrm{1}-\alpha ）\right|+\left|W（\mathrm{1}-\alpha ）\right|}\mathrm{。}（15）

当采样点数N足够大时， 式（7）化简为

A_h=N^{-\mathrm{1}}（\mathrm{2}{y}_{\mathrm{1}}+y_{\mathrm{2}}+y_{\mathrm{0}}）\nu (\alpha )\mathrm{。}（16）

利用Matlab多项式拟合可得各窗三谱线插值修正公式的逼近式如下：

Hann窗

\begin{array}{c}\alpha =\mathrm{0.666}\mathrm{669}\mathrm{38}\beta -\mathrm{0.073}\mathrm{628}\mathrm{64}{\beta }^{\mathrm{3}}+\\ \mathrm{0.021}\mathrm{707}\mathrm{41}{\beta }^{\mathrm{5}}+\mathrm{0.000}\mathrm{322}\mathrm{8}{\beta }^{\mathrm{7}},\end{array}

\begin{array}{c}\nu （\alpha ）=\mathrm{1.333}\mathrm{333}\mathrm{36}+\mathrm{0.526}\mathrm{620}\mathrm{90}{\alpha }^{\mathrm{2}}+\\ \mathrm{0.119}\mathrm{295}\mathrm{68}{\alpha }^{\mathrm{4}}+\mathrm{0.025}\mathrm{268}\mathrm{85}{\alpha }^{\mathrm{6}}\mathrm{。}\end{array}

Blackman窗

\begin{array}{c}\alpha =\mathrm{0.782}\mathrm{609}\mathrm{49}\beta -\mathrm{0.077}\mathrm{639}\mathrm{06}{\beta }^{\mathrm{3}}+\\ \mathrm{0.019}\mathrm{248}\mathrm{47}{\beta }^{\mathrm{5}}+\mathrm{0.000}\mathrm{670}\mathrm{99}{\beta }^{\mathrm{7}},\end{array}

\begin{array}{c}\nu （\alpha ）=\mathrm{1.492}\mathrm{537}\mathrm{32}+\mathrm{0.490}\mathrm{463}\mathrm{86}{\alpha }^{\mathrm{2}}+\\ \mathrm{0.088}\mathrm{529}\mathrm{77}{\alpha }^{\mathrm{4}}+\mathrm{0.013}\mathrm{714}\mathrm{84}{\alpha }^{\mathrm{6}}\mathrm{。}\end{array}

Balckmanharris窗

\begin{array}{c}\alpha =\mathrm{0.938}\mathrm{918}\mathrm{9}\beta -\mathrm{0.082}\mathrm{038}\mathrm{1}{\beta }^{\mathrm{3}}+\\ \mathrm{0.015}\mathrm{411}\mathrm{5}{\beta }^{\mathrm{5}}-\mathrm{0.003}\mathrm{169}\mathrm{9}{\beta }^{\mathrm{7}},\end{array}

\begin{array}{c}\nu （\alpha ）=\mathrm{1.658}\mathrm{663}\mathrm{6}+\mathrm{0.448}\mathrm{653}\mathrm{8}{\alpha }^{\mathrm{2}}+\\ \mathrm{0.064}\mathrm{787}\mathrm{3}{\alpha }^{\mathrm{4}}+\mathrm{0.006}\mathrm{938}\mathrm{3}{\alpha }^{\mathrm{6}}\mathrm{。}\end{array}

Nuttall窗

\begin{array}{c}\alpha =\mathrm{1.014}\mathrm{677}\mathrm{93}\beta -\mathrm{0.085}\mathrm{167}\mathrm{55}{\beta }^{\mathrm{3}}+\\ \mathrm{0.016}\mathrm{236}\mathrm{90}{\beta }^{\mathrm{5}}-\mathrm{0.002}\mathrm{888}\mathrm{39}{\beta }^{\mathrm{7}},\end{array}

\begin{array}{c}\nu （\alpha ）=\mathrm{1.724}\mathrm{338}\mathrm{599}\mathrm{3}+\mathrm{0.430}\mathrm{762}\mathrm{682}\mathrm{3}{\alpha }^{\mathrm{2}}+\\ \mathrm{0.056}\mathrm{560}\mathrm{313}\mathrm{4}{\alpha }^{\mathrm{4}}+\mathrm{0.005}\mathrm{697}\mathrm{239}{\alpha }^{\mathrm{6}}\mathrm{。}\end{array}

3.2 四谱线插值

如图 2 所示， 设实际峰值频点k附近的最大的4根峰值谱线分别为k₀， k₁， k₂， k₃（k_i+1=k_i +1）， 其中， k₁， k₂为k左右两侧最大的两条谱线。幅值y_{\mathrm{0}}=\left|X(k_{\mathrm{0}})\right|，y_{\mathrm{1}}=\left|X(k_{\mathrm{1}})\right|， y_{\mathrm{2}}=\left|X(k_{\mathrm{2}})\right|， y_{\mathrm{3}}=\left|X(k_{\mathrm{3}})\right|。

同三谱线插值\alpha \in [-\mathrm{0.5,0.5}]， 给k_{\mathrm{1}}， k_{\mathrm{2}}谱线幅值最大权值后有

\beta =\frac{y_{\mathrm{3}}+\mathrm{2}{y}_{\mathrm{2}}-\mathrm{2}{y}_{\mathrm{1}}-y_{\mathrm{0}}}{y_{\mathrm{3}}+\mathrm{2}{y}_{\mathrm{2}}+\mathrm{2}{y}_{\mathrm{1}}+y_{\mathrm{0}}}\mathrm{。}（17）

由式（11）得

\beta =\frac{\left|W（-\alpha +\mathrm{1.5}）\right|+\mathrm{2}\left|W（-\alpha +\mathrm{0.5}）\right|-\mathrm{2}W\left|-\alpha -\mathrm{0.5}\right|-W\left|-\alpha -\mathrm{1.5}\right|}{\left|W（-\alpha +\mathrm{1.5}）\right|+\mathrm{2}\left|W（-\alpha +\mathrm{0.5}）\right|+\mathrm{2}W\left|-\alpha -\mathrm{0.5}\right|+W\left|-\alpha -\mathrm{1.5}\right|}\mathrm{。}（18）

其反函数为\alpha =g^{-\mathrm{1}}(\beta )， 通过多项式拟合方法求得修正公式， 由\beta可求出参数\alpha。频率与相位的修正公式分别为

f_h=(k_{\mathrm{1}}+\alpha +\mathrm{0.5})\mathrm{\Delta }f,（19）

{\varphi }_h=\mathrm{a}\mathrm{r}\mathrm{g}[y(k_{\mathrm{1}})]+\frac{\mathrm{\pi }}{\mathrm{2}}-\mathrm{\pi }(\alpha +\mathrm{0.5})\mathrm{。}（20）

谐波信号的幅值修正需要对4条谱线进行加权平均计算， k_{\mathrm{1}}， k_{\mathrm{2}}为距离实际峰值谱线最大的两根谱线， 在计算式中应该给予谱线k_{\mathrm{1}}， k_{\mathrm{2}}的幅值y_{\mathrm{1}}， y_{\mathrm{2}}较大的权重， 由式（11）得

A_h=\frac{\mathrm{2}（y_{\mathrm{3}}+\mathrm{2}{y}_{\mathrm{2}}-\mathrm{2}{y}_{\mathrm{1}}-y_{\mathrm{0}}）}{\left|W（-\alpha +\mathrm{1.5}）\right|+\mathrm{2}\left|W（-\alpha +\mathrm{0.5}）\right|+\mathrm{2}\left|W（-\alpha -\mathrm{0.5}）\right|+\left|W（-\alpha -\mathrm{1.5}）\right|}\mathrm{。}（21）

当采样点数N足够大时， 上式化简为

A_h=N^{-\mathrm{1}}（y_{\mathrm{3}}+\mathrm{2}{y}_{\mathrm{2}}-y_{\mathrm{0}}-\mathrm{2}{y}_{\mathrm{1}}）\nu (\alpha )\mathrm{。}（22）

用Matlab多项式拟合可得各窗四谱线插值修正公式的逼近式如下：

1）　Hann窗

\begin{array}{c}\alpha =\mathrm{1.130}\mathrm{137}\mathrm{019}\beta +\mathrm{0.184}\mathrm{063}\mathrm{844}\mathrm{41}{\beta }^{\mathrm{3}}+\\ \mathrm{0.0729}\mathrm{659}\mathrm{84}{\beta }^{\mathrm{5}}+\mathrm{0.036}\mathrm{296}\mathrm{392}{\beta }^{\mathrm{7}}+\\ \mathrm{0.031}\mathrm{193}\mathrm{281}\mathrm{4}{\beta }^{\mathrm{9}},\end{array}

\begin{array}{c}\nu （\alpha ）=\mathrm{1.070}\mathrm{997}\mathrm{483}+\mathrm{0.275}\mathrm{923}\mathrm{44}{\alpha }^{\mathrm{2}}+\\ \mathrm{0.044}\mathrm{609}\mathrm{484}{\alpha }^{\mathrm{4}}+\mathrm{0.005}\mathrm{801}\mathrm{698}\mathrm{4}{\alpha }^{\mathrm{6}}\mathrm{。}\end{array}

2）　Blackman窗

\begin{array}{c}\alpha =\mathrm{1.446}\mathrm{489}\mathrm{078}\mathrm{8}\beta +\mathrm{0.293}\mathrm{253}\mathrm{038}\mathrm{9}{\beta }^{\mathrm{3}}+\\ \mathrm{0.139}\mathrm{233}\mathrm{113}\mathrm{02}{\beta }^{\mathrm{5}}+\mathrm{0.083}\mathrm{680}\mathrm{045}{\beta }^{\mathrm{7}}+\\ \mathrm{0.074}\mathrm{770}\mathrm{475}\mathrm{4}{\beta }^{\mathrm{9}},\end{array}

\begin{array}{c}\nu （\alpha ）=\mathrm{1.157}\mathrm{591}\mathrm{334}\mathrm{6}+\mathrm{0.268}\mathrm{164}\mathrm{009}{\alpha }^{\mathrm{2}}+\\ \mathrm{0.036}\mathrm{806}\mathrm{919}\mathrm{73}{\alpha }^{\mathrm{4}}+\mathrm{0.003}\mathrm{908}\mathrm{211}\mathrm{1}{\alpha }^{\mathrm{6}}\mathrm{。}\end{array}

3）　Balckmanharris窗

\begin{array}{c}\alpha =\mathrm{1.859}\mathrm{820}\mathrm{740}\mathrm{4}\beta +\mathrm{0.402}\mathrm{993}\mathrm{487}\mathrm{55}{\beta }^{\mathrm{3}}+\\ \mathrm{0.197}\mathrm{582}\mathrm{594}\mathrm{23}{\beta }^{\mathrm{5}}+\mathrm{0.122}\mathrm{972}\mathrm{274}\mathrm{27}{\beta }^{\mathrm{7}}+\\ \mathrm{0.100}\mathrm{804}\mathrm{903}\mathrm{6}{\beta }^{\mathrm{9}},\end{array}\begin{array}{c}\nu （\alpha ）=\mathrm{1.249}\mathrm{143}\mathrm{970}\mathrm{5}+\mathrm{0.256}\mathrm{990}\mathrm{337}{\alpha }^{\mathrm{2}}+\\ \mathrm{0.029}\mathrm{757}\mathrm{66}{\alpha }^{\mathrm{4}}+\mathrm{0.002}\mathrm{591}\mathrm{251}{\alpha }^{\mathrm{6}}\mathrm{。}\end{array}

4）　Nuttall窗

\begin{array}{c}\alpha =\mathrm{2.053}\mathrm{604}\mathrm{089}\mathrm{138}\beta +\mathrm{0.416}\mathrm{221}\mathrm{004}\mathrm{4}{\beta }^{\mathrm{3}}+\\ \mathrm{0.197}\mathrm{270}\mathrm{849}\mathrm{9}{\beta }^{\mathrm{5}}+\mathrm{0.120}\mathrm{451}\mathrm{124}\mathrm{56}{\beta }^{\mathrm{7}}+\\ \mathrm{0.094}\mathrm{324}\mathrm{459}\mathrm{76}{\beta }^{\mathrm{9}},\end{array}\begin{array}{c}\nu （\alpha ）=\mathrm{1.285}\mathrm{445}\mathrm{934}\mathrm{234}+\mathrm{0.251}\mathrm{636}\mathrm{666}\mathrm{1}{\alpha }^{\mathrm{2}}+\\ \mathrm{0.027}\mathrm{266}\mathrm{645}\mathrm{1}{\alpha }^{\mathrm{4}}+\mathrm{0.002}\mathrm{212}\mathrm{455}\mathrm{44}{\alpha }^{\mathrm{6}}\mathrm{。}\end{array}

4 Matlab仿真实验分析

4.1 Matlab工具

Matlab有丰富的功能、 工具箱及函数库， 是数学类科技软件中功能非常强大的数学仿真软件， 广泛使用于科学与工程领域， 非常适合于信号处理、 数值计算、 绘图等， 可以很好地解决仿真问题。此次仿真所用到的Matlab工具版本是Matlab 2022b。

4.2 信号模型的建立与仿真

为了测量Hann窗、 Blackman窗、 Blackmanharris窗、 Nuttall窗三谱线插值修正与四谱线插值修正对谐波信号幅值、 相位、 频率检测精度的影响， 根据电网信号的特点， 建立如下含有不同谐波的信号

x(\mathrm{t})={\displaystyle \sum _{h=\mathrm{0}}^{\mathrm{6}}}{A}_{h}cos\left(\mathrm{2}\right.\pi f_{h}t+{\varphi }_h),（23）

其中， 基波频率为50 Hz， 采样频率为3 200 Hz， 采样点数是1 024个。构建的两个信号的基波与各谐波信号的相位、 频率、 幅值等参数信息如表 1 和表 2 所示。使用Matlab仿真工具对该信号进行仿真实验， 不同余旋窗的不同谱线插值算法仿真所测得不同谐波次数的幅值、 频率、 相位相对误差如表 3~表 8 所示。

通过对信号Ⅰ与信号Ⅱ仿真结果的表格数据分析对比可知， 对于信号的幅值、 相位、 频率， 相同窗函数的四谱线插值修正结果的误差普遍小于三谱线插值修正的误差， 表明四谱线的修正算法要优于三谱线修正算法。对比不同余旋窗修正结果可以看到， Nuttall窗函数对信号的幅值、 频率、 相位检测结果的误差明显小于其它余旋窗函数， 加Nuttall窗进行四谱线插值修正后相比于其三谱线的修正算法对信号的幅值、 频率、 相位的测量精度都有所提高。从表 3、 表 4 可以看出Nuttall四线插值的幅值误差相比Nuttall三线插值都有所降低， 最低减小了5.71×10^-8（第四次谐波）， 从表 5、 表 6 可以看出Nuttall四线插值的频率误差相比Nuttall三线插值普遍减小了10^-1 ~10^-6数量级， 从表 7、 表 8 可以看出Nuttall四线插值的相位误差相比Nuttall三线插值普遍降低了10^-1量级， 结果表明Nuttall窗四谱线插值修正算法可以很好地减小非同步采样带来的问题。

5 结　语

本文对稳态谐波信号进行经典傅里叶变换谐波分析时， 对信号进行加余旋窗截断处理以减弱频谱泄露的影响， 并对其进行谱线插值修正抑制栏栅效应。通过最小二乘法及使用Matlab进行多项式拟合从而推导出不同余旋窗的不同插值算法相应的谐波幅值、 频率及相位的修正公式。利用Matlab进行仿真实验， 由检测结果可知基于Nuttall余弦窗四谱线插值修正算法相较于其它余旋窗插值修正算法的检测准确度更高，比其它算法抑制频谱泄漏以及减弱栅栏效应影响的效果更好， 可达到精确检测谐波的目的。

参考文献

[1] 李晓龙， 陈建青， 刘波. 基于汉宁窗的插值FFT算法在11号数据链测频中的应用［J］. 舰船电子工程， 2018， 38（5）： 37-40.

[2] LI Xiaolong， CHEN Jianqing， LIU Bo. Application of Link11 measuring frequency based on hanning window interpolated FFT algorithm［J］. Ship Electronic Engineering， 2018， 38（5）： 37-40. （in Chinese）

[3] WEN H， ZHANG J， MENG Z， et al. Harmonic estimation using symmetrical interpolation FFT based on triangular self-convolution window［J］. IEEE Transactions on Industrial Informatics， 2015， 11（1）： 16-26.

[4] 吕文杰. 电力系统谐波检测及其抑制方法的研究［D］. 成都： 成都理工大学， 2012.

[5] 李永军， 马立元， 王天辉， 等. 一种新型组合窗及其在模态参数辨识中的应用［J］. 测试技术学报， 2011， 25（5）： 445-449.

[6] LI Yongjun， MA Liyuan， WANG Tianhui，et al. A new combined window and its application to module parameter identification［J］. Journal of Test and Measurement Technology， 2011， 25（5）： 445-449. （in Chinese）

[7] 周俊， 王小海， 祁才君. 基于Blackman窗函数的插值FFT在电网谐波信号分析中的应用［J］. 浙江大学学报（理学版）， 2006， 33（6）： 650-653.

[8] ZHOU Jun， WANG Xiaohai， QI Caijun. Estimation of electrical harmonic parameters by using the interpolated FFT algorithm based on Blackman window［J］. Journal of Zhejiang University （Science Edition）， 2006， 33（6）： 650-653. （in Chinese）

[9] 周西峰， 赵蓉， 郭前岗. Blackman-Harris窗的插值FFT谐波分析与应用［J］. 电测与仪表， 2014， 51（11）： 81-85.

[10] ZHOU Xifeng， ZHAO Rong， GUO Qiangang. Blackman-Harris window based interpolation FFT harmonic analysis and its application［J］. Electrical Measurement & Instrumentation， 2014， 51（11）： 81-85. （in Chinese）

[11] 卿柏元， 滕召胜， 高云鹏， 等. 基于Nuttall窗双谱线插值FFT的电力谐波分析方法［J］. 中国电机工程学报， 2008， 28（25）： 153-158.

[12] QING Baiyuan， TENG Zhaosheng， GAO Yunpeng， et al. An approach for electrical harmonic analysis based on nuttall window double-spectrum-line interpolation FFT［J］. Proceedings of the CSEE， 2008， 28（25）： 153-158. （in Chinese）

[13] 李威， 李春红， 杨鑫. 基于Nuttall窗三谱线插值FFT谐波测量［J］. 价值工程， 2017， 36（4）： 79-81.

[14] LI Wei， LI Chunhong， YANG Xin. Harmonic analysis based on four-term fifth derivative nuttall window and triple-spectrum-line interpolation FFT［J］. Value Engineering， 2017， 36（4）： 79-81. （in Chinese）

[15] 杨阳. 窗函数特性及加窗插值FFT算法的研究［D］. 郑州： 郑州大学， 2010.

[16] 温和， 滕召胜， 王永， 等. 改进加窗插值FFT动态谐波分析算法及应用［J］. 电工技术学报， 2012， 27（12）： 270-277.

[17] WEN He， TENG Zhaosheng， WANG Yong， et al.Improved windowed interpolation FFT algorithm and application for power harmonic analysis［J］. Transactions of China Electrotechnical Society， 2012， 27（12）： 270-277. （in Chinese）

[18] 沈杰， 韩峰，梅秀庄. 窗函数对信号稀疏特性的影响［J］. 测试技术学报， 2009， 23（2）： 151-155.

[19] SHEN Jie， HAN Feng， MEI Xiuzhuang. Effect of window function on signal’s sparse characteristics ［J］. Journal of Test and Measurement Technology， 2009， 23（2）： 151-155. （in Chinese）

[20] 荆雪君， 智泽英. 基于FFT的间谐波优化检测算法研究［J］. 电气自动化， 2018， 40（5）： 51-54.

[21] JING Xuejun， ZHI Zeying. Interharmonics optimization detection algorithm based on FFT ［J］. Electrical Automation， 2018， 40（5）： 51-54. （in Chinese）

[22] 陆东亮. 高精度电力谐波检测算法研究［D］. 武汉： 武汉纺织大学， 2022.

[23] 耿莉莉. 基于加窗插值FFT和原子分解的间谐波检测算法［D］. 郑州： 郑州大学， 2017.

[24] 赵心迪. 一种基于FFT的电力系统谐波分析改进方法研究［J］. 煤炭工程， 2017， 49（4）： 98-100.

[25] ZHAO Xindi. An improved algorithm for harmonic analysis of power system based on FFT algorithm［J］. Coal Engineering， 2017， 49（4）： 98-100. （in Chinese）

[26] 张强， 王海云， 王维庆， 等. 基于混合卷积窗三谱线插值的介损角测量方法［J］. 电力电容器与无功补偿， 2020， 41（6）： 125-131.

[27] ZHANG Qiang， WANG Haiyun， WANG Weiqing， et al. Dielectric loss angle measurement method based on mixed convolution window with three spectral lines interpolation［J］. Power Capacitor & Reactive Power Compensation， 2020， 41（6）： 125-131. （in Chinese）

[28] 牛胜锁， 梁志瑞， 张建华， 等. 基于三谱线插值FFT的电力谐波分析算法［J］. 中国电机工程学报， 2012， 32（16）： 130-136.

[29] NIU Shengsuo， LIANG Zhirui， ZHANG Jianhua， et al. An algorithm for electrical harmonic analysis based on triple-spectrum-line interpolation FFT［J］. Proceedings of the CSEE， 2012， 32（16）： 130-136. （in Chinese）

[30] 蔡晓峰， 张鸿博， 鲁改凤. 应用三谱线插值FFT分析电力谐波的改进算法［J］. 电力系统保护与控制， 2015， 43（2）： 33-39.

[31] CAI Xiaofeng， ZHANG Hongbo， LU Gaifeng. Improvement algorithm for harmonic analysis of power system using triple-spectrum-line interpolation algorithm based on window FFT［J］. Power System Protection and Control， 2015， 43（2）： 33-39. （in Chinese）

[32] 李振华， 胡廷和， 杜亚伟， 等. 基于窗函数和谱线插值理论的谐波检测方法［J］. 电力系统保护与控制， 2019， 47（22）： 78-88.

[33] LI Zhenhua， HU Tinghe， DU Yawei， etal. Harmonic detection method based on the theory of windows and spectrum line interpolation［J］. Power System Protection and Control， 2019， 47（22）： 78-88. （in Chinese）