CN103245832B - 基于快速s变换的谐波时频特性参数估计方法及分析仪 - Google Patents

Abstract

基于快速S变换的谐波时频特性参数估计方法及分析仪，本发明之基于快速S变换的谐波时频特性参数估计方法包括以下步骤：（1）模数转换；（2）低通滤波；（3）快速S变换；得到谐波时频特性参数结果；本发明还包括分析仪；本发明之基于快速S变换的谐波时频特性参数估计方法，计算量少，效率高、实用性强；具有从复杂畸变信号中准确提取谐波时频特征参数的能力；本发明之分析仪，结构简单，布局合理。

Description

基于快速S变换的谐波时频特性参数估计方法及分析仪

技术领域

本发明涉及一种电力系统谐波时频分析方法，尤其是涉及一种基于快速S变换的谐波时频特性参数估计方法

背景技术

谐波是电能质量问题的主要根源，其时频特征参数的准确分析是研究电网谐波问题的重要内容，是谐波潮流计算、谐波电能计量、绝缘设备监测、谐波补偿与抑制的最关键依据，对确保电网安全、经济运行具有重要意义。

现有谐波频域分析方法主要是FFT及其改进算法，但不能实现谐波时域、频域参数的全面分析，无法满足暂态或突变等特性的非平稳信号谐波信号快速、准确分析的要求。且，时域分析法主要是基于各种矢量变换（d-q变换、对称分量变换等）和瞬时无功功率理论的参数估计方法，其信息量少，计算量大。

谐波时频参数包括：各次谐波幅值、相位、起止时刻、持续时间等。目前现有的谐波时域、频域联合特征分析和估计方法有：

(1) 短时傅里叶变换 Dennis Gabor于1946年建立了短时傅里叶变换（Short Time Fourier Transform，STFT），STFT算法简单，在电网谐波谱分析、电压暂降扰动检测等方面取得了应用。但存在不可克服的缺陷：①同一信号不同窗函数的STFT结果相差很远；②信号分析的时、频分辨率在窗函数选定后无法改变，只能固定分辨率，无法兼顾高频信息和低频信息；③STFT的离散形式没有正交展开，很难实现高效算法。

（2）小波变换（Wavelet Transform，WT）是具有多分辨率特性的时频局部化分析方法。优势在于：①可以对信号的任意细节聚焦，适合处理突变信号②根据频带进行信息处理，结果较为稳定可靠。但仍存在诸多局限：①频域分辨率粗糙，可能导致各频带间严重的频率混叠；②引入尺度因子后，其结果仍不是真正意义上的时频谱；③只能定性分析而不能准确检测信号的幅值与谐波分量；④小波基函数无统一选择方法，分析结果受小波基函数的影响较大；⑤算法复杂，不适合嵌入式系统实现。

(3) S变换:1996年，Stockwell等人将连续WT和STFT结合，提出了S变换时频分析方法，将一维时间信号变成二维包含局部时频谱性质的时间-频率谱。S变换利用Gauss窗的高度与频率成正比、宽度随频率的增大而变小的特点，克服了STFT窗固定的缺陷，可获得某一时刻的频率信息和某一频率点的信号幅值信息，实现动态信号的时频特性分析。近年来，S变换备受重视，但处于发展中的S变换的不足也日渐显现。其中， S变换及其改进算法的时频矩阵信息量非常大，计算繁琐，难以嵌入式实现。

发明内容

本发明所要解决的技术问题是，克服现有技术存在的上述缺陷，提供一种信息量

大，计算简单，易嵌入式的基于快速S变换的谐波时频特性参数估计方法及分析仪。

本发明解决其技术问题所采用的技术方案是，

本发明之基于快速S变换的谐波时频特性参数估计方法，包括以下步骤：

（1）低通滤波：将被测信号输入电力系统分析仪内，根据模数转换速率和谐波分析要求，确定被测信号进行谐波分析的最高次数和低通数字滤波器的电阻电容取值，将被测信号低通滤波处理；

（2）模数转换：将步骤（1）低通滤波好的被测信号经模数转换器高速模数转换为数字信号；

（3）快速S变换：将步骤（2）转换好的被测数字信号送入数字信号处理器内，完成构建离散多高斯窗，进行快速S变换，进行时频谐波幅值、相角、起止、突变时刻参数分析；

得到谐波时频特性参数结果。

进一步，步骤（3）中，所述快速S变换的方法包括以下步骤：

（1）FFT运算：将采样好的被测数字信号x(t)通过FFT运算求得到信号频谱;

信号的连续S变换S(τ，f )定义如下

$\left\{\begin{array}{c}S(\mathrm{\τ},t)={\∫}_{-\∞}^{\∞}x\left(t\right)w(\mathrm{\τ}-t,f)e^{-i2\mathrm{\πft}}\mathrm{dt}\\ w(\mathrm{\τ}-t,f)=\frac{\left|f\right|}{\sqrt{2\mathrm{\π}}}{e}^{-\frac{f^2{(\mathrm{\τ}-t)}^2}{2}}\end{array}\right.---\left(1\right)$

式中，w(τ-t，f )为高斯窗函数，其宽度频率自适应调整，t、τ为时移因子；f为频率，x(t)为被测信号。

由式（1）可知，令τ=mT_s，f=k/NT_s，则x(nT_s)的一维离散S变换为

$\left\{\begin{array}{cc}S(m{T}_s,\frac{k}{{\mathrm{NT}}_s})=\underset{r=0}{\overset{N-1}{\mathrm{\Σ}}}X\left[\frac{r+k}{{\mathrm{NT}}_s}\right]G(r,k)e^{j\frac{2\mathrm{\πmr}}{N}}& k\≠0\\ S({\mathrm{mT}}_s,0)=\underset{r=0}{\overset{N-1}{\mathrm{\Σ}}}X\left(\frac{r}{{\mathrm{NT}}_s}\right)& k=0\end{array}\right.---\left(2\right)$

式中，N为采样点数；

k、m、n、r取值分别为0,1,2,…,N-1；

T_s为采样时间间隔;

X(k/NT_s)为信号x(nT_s)的离散傅里叶频谱;

 G(r, k)为高斯窗的傅里叶频谱;

且

$\left\{\begin{array}{c}X\left(\frac{k}{{\mathrm{NT}}_s}\right)=\frac{1}{N}\underset{n=0}{\overset{N-1}{\mathrm{\Σ}}}x\left({\mathrm{nT}}_s\right)e^{-(i2\mathrm{\πkn}/N)}\\ G(r,k)=e^{-\frac{2{\mathrm{\π}}^2{r}^2}{k^2}}\end{array}\right.---\left(3\right)$

式中，N为采样点数；

k、m、n、r取值分别为0,1,2,…,N-1；

T_s为采样时间间隔;

X(k/NT_s)为信号x(nT_s)的离散傅里叶频谱;

G(r, k)为高斯窗的傅里叶频谱;

（2）特征频率点判断：因X(k/NT_s)各谱线对应频率与特征频率存在偏差，对频谱进行逐求导X '(k/NT_s)，求得频谱各极大值点k^’，极大值对应频率最接近实际特征频率，应为信号的特征频率；为消除频谱泄漏引起的误判断，极大值点需满足|X(k_i'/NT_s)|>ε，ε为设定阈值；

令符合条件极大值点构成的一维矩阵=[k₁ ', k₂ '…… k_q-1 ',k_q ' ](q为特性频率点数) 。

（3）窗函数自适应调整：步骤（1）中w(τ-t，f )为高斯窗函数，其宽度频率自适应调整；

（4）频谱平移频谱相乘：将步骤（2）满足条件的特定频率点k '，平移频谱X(k'/NT_s)到X(k'/NT_s+r)，计算X(k'/NT_s+r) G(r, k'/NT_s)；

（5）IFFT运算：将步骤（4）相乘的频谱X(k'/NT_s+r) G(r, k'/NT_s)进行IFFT运算，得到S(m, k')、A(m, k')和φ(m, k')，计算总量为O(rNlog₂N)；

S(m, k')为S变换矩阵，记为S矩阵；S矩阵的行对应采样时刻，列对应频率，S变换是复数变换，S矩阵是复数矩阵，S矩阵可表示为

；

式中，m,k分别代表S矩阵的行和列， A(m, k')为S矩阵的幅值矩阵，即S模矩阵；φ(m, k')为S矩阵的相位矩阵；A(m, k')和φ(m,k')的行向量分别表示信号某一采样时刻的幅值和相位随频率变化的分布，其列向量分别表示信号某一频率处的幅值和相位随时间变化的分布；

（6）复数时频矩阵分析：将步骤（4）运算得的S(m, k')、A(m, k')和φ(m, k')构成复数时频矩阵，根据相位和幅值信息，对各谐波进行分段，获取各时间段谐波幅值和相位信息，得到谐波的时频特性参数；

重复步骤（2）-（6），直到计算出所有谐波频率点对应的快速S变换。

步骤（2）中特性频率点判断方法：

特性频率点需满足公式（5）（6）

$X\left(\frac{k_{i-1}}{{\mathrm{NT}}_s}\right)<X\left(\frac{k_i}{{\mathrm{NT}}_s}\right)>X\left(\frac{k_{i+1}}{{\mathrm{NT}}_s}\right)---\left(5\right)$

$\left|X\left(\frac{k_i}{{\mathrm{NT}}_s}\right)\right|>\mathrm{\&epsiv;}---\left(6\right)$

式中ε为设定阈值，Ts为采样周期，ki为谱线号，X(k/NT_s)为信号x(nT_s)的离散傅里叶频谱，N为采样点数。

本发明之谐波特性参数分析仪，包括电源模块，所述电源模块分别与模数转换器（ADC）和数字信号处理器(DSP)电连接，所述模数转换器（ADC）分别连有电压测量信号变换电路和电流测量信号变换电路及数字信号处理器(DSP)，所述数字信号处理器(DSP)分别与复位模块、调试JTAG接口、同步动态随机存储器（SDRAM）、FLASH和有源振晶电连接。

进一步，所述模数转换器（ADC）优选采用TI公司生产的ADS8364。

进一步，所述数字信号处理采用TI公司生产的TMS320VC6745处理器，主要参数选择如下：

（1）采样速率：8kHz

（2）电网基波频率：在50Hz左右波动；

（3）FFT运算数据长度：N=2048。

进一步，所述信号数据长度N是进行一次快速S变换的数据长度，根据谐波检测精度和计算机或嵌入式系统设备的运行速度综合考虑确定，

由谐波次数h确定谐波分析采样频率f_s，

$f_s\&GreaterEqual;2\&times;h\&times;50$

由采样频率和频谱分辨率B确定采样长度N，其具体公式为

$N=2^y\&GreaterEqual;f_s/B$ （式中y为正整数,5<y<20）；

保证DSP计算时间小于一个工频电压周期，DSP计算时间为

$t=rN\&times;cycl{e}_{FFT}/f_{core}$

（式中q为特征频率点数，cycle_FFT为DSP进行一次FFT的时钟周期，f_core为DSP内核主频）。

本发明之基于快速S变换的谐波时频特性参数估计方法，基于频谱极大值检测和阈值判断的快速S变换算法，可对特征频率点进行高斯加窗和时域幅值曲线求取，计算量少，效率高和实用性强；快速S变换，算法直观、能表达各次谐波的时频特征，对于区段稳定的谐波信号，其频率特性曲线可分别准确得到幅值、相位、频率等特征参数；对于谐波幅值突变信号，时域特性曲线可分别表明相应时域幅值的变化过程，具有从复杂畸变信号中准确提取谐波时频特征参数的能力；本发明之分析仪，结构简单，布局合理。

附图说明

图1为本发明一实施例快速S变换步骤（3）窗函数自适应调整示意图；

图2为本发明一实施例的方法流程图；

图3为图2所示基于快速S变换的谐波时频特性参数估计方法中的快速S变换的流程示意图；

图4为本发明一实施例稳定谐波信号时频分析：

（a）原始信号；（b）时频三维分布；（c）时频能量分布；（d）基波幅值曲线；

图5为本发明一实施例时变谐波信号时频分析：

（a）原始信号；（b）时频三维分布；（c）时频能量分布；（d）基波幅值曲线；

图6为本发明分析仪结构示意图。

具体实施方式

实施例

参照图1，图2，本实施例包括以下步骤：

（1）低通滤波：将被测信号输入电力系统分析仪内，根据模数转换速率和谐波分析要求，确定被测信号进行谐波分析的最高次数和低通数字滤波器的电阻电容取值，将被测信号低通滤波处理；

（2）模数转换：将步骤（1）低通滤波好的被测信号经模数转换器高速模数转换为数字信号；

（3）快速S变换：将步骤（2）转换好的被测数字信号送入数字信号处理器内，完成构建离散高斯窗，进行快速S变换，进行时频谐波幅值、相角、起止、突变时刻参数分析；

得到谐波时频特性参数结果。

参照附图3，步骤（3）中，快速S变换的方法包括以下步骤：

（1）FFT运算：将采样好的被测数字信号x(t)通过FFT运算求得到信号频谱：

信号的连续S变换S(τ，f )定义如下

$\left\{\begin{array}{c}S(\mathrm{\&tau;},t)={\&Integral;}_{-\&infin;}^{\&infin;}x\left(t\right)w(\mathrm{\&tau;}-t,f)e^{-i2\mathrm{\&pi;ft}}\mathrm{dt}\\ w(\mathrm{\&tau;}-t,f)=\frac{\left|f\right|}{\sqrt{2\mathrm{\&pi;}}}{e}^{-\frac{f^2{(\mathrm{\&tau;}-t)}^2}{2}}\end{array}\right.$

式中，w(τ-t，f )为高斯窗函数，其宽度频率自适应调整， τ为时移因子；f为频率。

令τ=mT_s，f=k/NT_s，则x(nT_s)的一维离散S变换为

$\left\{\begin{array}{cc}S(m{T}_s,\frac{k}{{\mathrm{NT}}_s})=\underset{r=0}{\overset{N-1}{\mathrm{\&Sigma;}}}X\left[\frac{r+k}{{\mathrm{NT}}_s}\right]G(r,k)e^{j\frac{2\mathrm{\&pi;mr}}{N}}& k\&NotEqual;0\\ S({\mathrm{mT}}_s,0)=\underset{r=0}{\overset{N-1}{\mathrm{\&Sigma;}}}X\left(\frac{r}{{\mathrm{NT}}_s}\right)& k=0\end{array}\right.$

式中，N为采样点数；

k、m、n取值分别为0,1,2,…,N-1；

T_s为采样时间间隔，;

X(k/NT_s)为信号x(nT_s)的离散傅里叶频谱;

G(r, k)为高斯窗的傅里叶频谱;

且

$\left\{\begin{array}{c}X\left(\frac{k}{{\mathrm{NT}}_s}\right)=\frac{1}{N}\underset{n=0}{\overset{N-1}{\mathrm{\&Sigma;}}}x\left({\mathrm{nT}}_s\right)e^{-(i2\mathrm{\&pi;kn}/N)}\\ G(r,k)=e^{-\frac{2{\mathrm{\&pi;}}^2{r}^2}{k^2}}\end{array}\right.$

式中，N为采样点数；

k、m、n、r取值分别为0,1,2,…,N-1；

T_s为采样时间间隔;

X(k/NT_s)为信号x(nT_s)的离散傅里叶频谱;

 G(r, k)为高斯窗的傅里叶频谱;

（2）特征频率点判断：因X(k/NT_s)各谱线对应频率与特征频率存在偏差，对频谱进行逐求导X '(k/NT_s)，求得频谱各极大值点k^’，极大值对应频率最接近实际特征频率，应为信号的特征频率；为消除频谱泄漏引起的误判断，极大值点需满足|X(k_i'/NT_s)|>ε，ε为设定阈值；令符合条件极大值点构成的一维矩阵=[k₁ ', k₂ '… k_q-1 ',k_q'] (q为特性频率点数)

（3）窗函数自适应调整：步骤（1）中w(τ-t，f )为高斯窗函数，其宽度频率自适应调整，如图1所示；

（4）频谱平移频谱相乘：将步骤（2）满足条件的特定频率点k '，平移频谱X(k'/NT_s)到X(k'/NT_s+r)，计算X(k'/NT_s+r) G(r, k'/NT_s)；

（5）IFFT运算：将步骤（4）相乘的频谱X(k'/NT_s+r) G(r, k'/NT_s)进行IFFT运算，得到S(m, k')、A(m, k')和φ(m, k')，计算总量为O(rNlog₂N)；

S(m, k')为S变换矩阵，记为S矩阵；S矩阵的行对应采样时刻，列对应频率，S变换是复数变换，S矩阵是复数矩阵，S矩阵可表示为

；

式中，m,k分别代表S矩阵的行和列， A(m, k')为S矩阵的幅值矩阵，即S模矩阵；φ(m, k')为S矩阵的相位矩阵；A(m, k')和φ(m,k')的行向量分别表示信号某一采样时刻的幅值和相位随频率变化的分布，其列向量分别表示信号某一频率处的幅值和相位随时间变化的分布；

（6）复数时频矩阵分析：将步骤（4）运算得的S(m, k')、A(m, k')和φ(m, k')构成复数时频矩阵，根据相位和幅值信息，对各谐波进行分段，获取各时间段谐波幅值和相位信息，得到谐波的时频特性参数；

重复步骤（2）-（6），直到计算出所有谐波频率点对应的快速S变换。

步骤（2）中特性频率点判断方法：

特性频率点需满足公式以下两式：

$X\left(\frac{k_{i-1}}{{\mathrm{NT}}_s}\right)<X\left(\frac{k_i}{{\mathrm{NT}}_s}\right)>X\left(\frac{k_{i+1}}{{\mathrm{NT}}_s}\right)$  ;

$\left|X\left(\frac{k_i}{{\mathrm{NT}}_s}\right)\right|>\mathrm{\&epsiv;}$ ;

式中ε为设定阈值，Ts为采样周期，ki为谱线号，X(k/NT_s)为信号x(nT_s)的离散傅里叶频谱，N为采样点数。

离散信号x(n)进行S变换需要进行FFT和IFFT运算，且每个频率点的复数乘加运算次数和FFT运算量相同，运算量大。已知N点的FFT计算量为O(Nlog₂N)，则N点的S变换计算量为O(N²log₂N)，执行一次完整的N点S变换需要的时间为

T=2×N×cycles_ftt/f_主频;

以德州仪器(TI)的6000系列DSP(主频720MHz)实现为例，执行一次长度分为512、1024的32位FFT和IFFT，分别需要13214、26286个时钟周期，完成一次S变换分别需62.23 ms、144.76ms，计算时间远远超过工频电压周期（20ms），不利于信号的实时快速分析。

若特征频率个数为q（谐波的实际发生个数），则快速S变换的计算量可降至O(qNlog₂N)，执行一次快速S变换的时间为

T=2×N×cycles_ftt/f_主频;

设q=20，对于长度分别为512、1024的定点16位FFT和IFFT，快速S变换的计算量仅为S变换的3.94%、1.97%，计算时间分别为2.269ms、3.074ms，远小于工频电网基波额定周期（20ms）便于嵌入式系统实现和在线实时检测应用。

利用DSP仿真器SEED-XDS560PLUS，在DSP开发环境CCS4.2的CLOCK功能下测试各算法的运行时间，处理器选用TMS320VC6745，实际测试结果如表1所示：

表1   各种算法计算时间

由表1可知，采用传统S变换分析谐波信号，计算量大，而采用本实施例提出的快速S变换算法，运算时间和FFT在一个数量级，且不影响谐波的测试精度。此外，FFT算法不具备时域分析能力，无法获取谐波起止时刻、时域幅值曲线等信息。

仿真实现：

采样频率f_s=8kHz，采样长度N=2048，信号基波频率为f₀=50Hz，仿真实验结果如图4、图5所示。

参照图4，输入信号稳定谐波信号包含基波和7、49次谐波，幅值、频率保持不变。其时频能量分布图（图中的横坐标为时间，纵坐标为归一化频率）显示辨识结果近似为三条直线，能量聚集度高，基波幅值曲线为常数。

参照图5，时变谐波信号包含基波和7、15和49次谐波，且其幅值在采样阶段发生波动。快速S变换原理算法能清楚地区分基波、各次谐波信号，消除了各次谐波间的相互影响，并能准确定位基波幅值突变和谐波发生的起止时刻。

仿真实验结果表明，快速S变换原理算法直观、有效地表达了各次谐波的时频特征。对于区段稳定的谐波信号，其频率特性曲线可分别准确得到幅值、相位、频率等特征参数。对于谐波幅值突变信号，时域特性曲线可分别表明相应时域幅值的变化过程。显然，快速S变换原理算法具有从复杂畸变信号中准确提取谐波时频特征参数的能力。

实现本实施例的电力系统谐波分析仪，包括电源模块1，所述电源模块1分别与模数转换器2和数字信号处理器3电连接，所述模数转换器1分别连有电压测量信号变换电路4和电流测量信号变换电路5及数字信号处理器3，所述数字信号处理器3分别与复位模块6、调试JTAG接口7、同步动态随机存储器8、FLASH9和有源振晶10电连接。

所述模数转换器1采用TI公司生产的ADS8364。

所述数字信号处理器3的CPU采用TI公司生产的TMS320VC6745，主要参数选择如下：

（1）采样速率：8kHz

（2）电网基波频率： 50Hz波动；

（3）FFT运算数据长度：N=2048。

模数转换器2用于对被测信号（电压或电流信号）进行高速模数转换，将被测信号转换为数字信号。

电压测量信号变换电路4和电流测量信号变换电路5能将被测信号中高频信号进行滤除的滤波算法，根据模数转换速率和谐波分析要求，确定被测信号进行谐波分析的最高次数和低通数字滤波器的各参数，对被测信号进行低通滤波处理，可以消除高次谐波对谐波分析结果的影响。

信号数据长度N是进行一次快速S变换的数据长度，N根据谐波检测精度和计算机或嵌入式系统设备的运行速度综合考虑确定，

由谐波次数h确定谐波分析采样频率f_s，

$f_s\&GreaterEqual;2\&times;h\&times;50$ ;

由采样频率和频谱分辨率B确定采样长度N，其具体公式为

$N=2^y\&GreaterEqual;f_s/B$ （式中y为正整数,5<y<20）；

保证DSP计算时间小于一个工频电压周期，DSP计算时间为

$t=qN\&times;cycl{e}_{FFT}/f_{core}$ ；

（式中q为特征频率点数，cycle_FFT为DSP进行一次FFT的时钟周期？f_core为DSP内核主频）。

基波与谐波参数检测误差结果如下：

表2   基波~21次谐波幅值相对误差(%)

谐波次数	1	2	3	4	5	6	7
								误  差	2.1 E-3	3.7 E-4	9.3 E-3	2.7 E-4	6.1 E-2	9.8 E-5	2.1 E-4
谐波次数	8	9	10	11	12	13	14
								误  差	1.8 E-4	3.1E-4	6.0 E-4	8.1 E-4	1.4 E-4	2.7 E-3	2.3 E-5
谐波次数	15	16	17	18	19	20	21
								误  差	6.8 E-3	2.9 E-3	2.1 E-4	6.8 E-3	4.4 E-5	6.9 E-4	1.8E-3

表3   基波~21次谐波初相位相对误差(%)

谐波次数	1	2	3	4	5	6	7
								误  差	2.1-5	1.2 E-8	2.9 E-6	2.8 E-6	2.1 E-7	2.3E-8	2.9E-7
谐波次数	8	9	10	11	12	13	14
								误  差	1.6-7	5.3E-6	1.9 E-7	1.3 E-5	7.7 E-6	1.5E-5	1.3 E-6
谐波次数	15	16	17	18	19	20	21
								误  差	2. 7E-6	3.5 E-6	1.5E-6	2.3E-6	4.4E-5	-2.4E-5	6.9 E-7

Claims (6)

1.基于快速S变换的谐波时频特性参数估计方法，其特征在于，包括以下步骤：

(1)低通滤波：将被测信号输入电力系统分析仪内，根据模数转换速率和谐波分析要求，确定被测信号进行谐波分析的最高次数和低通数字滤波器的电阻电容取值，将被测信号低通滤波处理；

(2)模数转换：将步骤(1)低通滤波好的被测信号经模数转换器高速模数转换为数字信号；

(3)快速S变换：将步骤(2)转换好的被测数字信号送入数字信号处理器内，完成构建离散多高斯窗，进行快速S变换，进行时频谐波幅值、相角、起止、突变时刻参数分析；

得到谐波时频特性参数结果；

步骤(3)中，所述快速S变换的方法包括以下步骤：

(1)FFT运算：将采样好的被测数字信号x(t)通过FFT运算求得到信号频谱，信号的连续S变换S(τ，f)定义如下

式中，w(τ-t，f)为高斯窗函数，其宽度频率自适应调整，τ为时移因子,f为频率，t为时间因子，x(t)为被测信号；

由上式可知，令τ＝mT_s，t＝k/NT_s，则x(nT_s)的一维离散S变换为

式中，N为信号采样长度；

k、m、n、r取值分别为0,1,2,…,N-1；

T_s为采样时间间隔；

X(k/NT_s)为信号x(nT_s)的离散傅里叶频谱；

G(r,k)为高斯窗的傅里叶频谱；

且

式中，N为信号采样长度；

k、m、n、r取值分别为0,1,2,…,N-1；

T_s为采样时间间隔；

X(k/NT_s)为信号x(nT_s)的离散傅里叶频谱；

G(r,k)为高斯窗的傅里叶频谱；

(2)特征频率点判断：因X(k/NT_s)各谱线对应频率与特征频率存在偏差，对频谱进行逐求导X'(k/NT_s)，求得频谱各极大值点k’，极大值对应信号的特征频率；为消除频谱泄漏引起的误判断，极大值点需满足|X(k_i'/NT_s)|>ε，ε为设定阈值，令符合条件极大值点构成的一维矩阵＝[k₁',k₂'……k_q-1',k_q']，q为特性频率点个数；

(3)窗函数自适应调整：步骤(1)中w(τ-t，f)为高斯窗函数，其宽度频率自适应调整；

(4)频谱平移频谱相乘：将步骤(2)满足条件的特定频率点k'，平移频谱X(k'/NT_s)到X(k'/NT_s+r)，计算X(k'/NT_s+r)G(r,k'/NT_s)；其中，X(k/NT_s)为信号x(nT_s) 的离散傅里叶频谱，G(r,k'/NT_s)为Gauss窗的离散傅里叶频谱T_s为采样时间间隔；

(5)IFFT运算：将步骤(4)相乘的频谱X(k'/NT_s+r)G(r,k'/NT_s)进行IFFT运算，得到S(m,k')、A(m,k')和计算总量为O(rNlog₂N)；

S(m,k')为S变换矩阵，记为S矩阵；S矩阵的行对应采样时刻，列对应频率，S变换是复数变换，S矩阵是复数矩阵，S矩阵可表示为

式中，m,k分别代表S矩阵的行和列，A(m,k')为S矩阵的幅值矩阵，即S模矩阵；为S矩阵的相位矩阵；A(m,k')和的行向量分别表示信号某一采样时刻的幅值和相位随频率变化的分布，其列向量分别表示信号某一频率处的幅值和相位随时间变化的分布；

(6)复数时频矩阵分析：将步骤(5)运算得的S(m,k')、A(m,k')和构成复数时频矩阵，根据相位和幅值信息，对各谐波进行分段，获取各时间段谐波幅值和相位信息，得到谐波的时频特性参数；

重复步骤(2)-(6)，直到计算出所有谐波频率点对应的快速S变换。

2.根据权利要求1所述的基于快速S变换的谐波时频特性参数估计方法，其特征在于，步骤(2)中，所述特性频率点的判断方法：

特性频率点需满足以下两个公式

式中ε为设定阈值，Ts为采样周期，ki为谱线号，X(k/NT_s)为信号x(nT_s)的离散傅里叶频谱，N为信号采样长度。

3.实现权利要求1或2所述基于快速S变换的谐波时频特性参数估计方法的分析仪，包括电源模块，所述电源模块分别与模数转换器和数字信号处理器电连接，其特征在于，所述模数转换器分别连有电压测量信号变换电路和电流测量信号变换电路及数字信号处理器，所述数字信号处理器分别与复位模块、调试JTAG接口、同步动态随机存储器、FLASH和有源振晶电连接。

4.根据权利要求3所述的分析仪，其特征在于，所述模数转换器采用TI公司生产的ADS8364。

5.根据权利要求3所述的的分析仪，其特征在于，所述数字信号处理器采用TI公司生产的TMS320VC6745，主要参数选择如下：

(1)采样速率：8kHz；

(2)电网基波频率：50Hz波动；

(3)FFT信号采样长度：N＝2048。

6.根据权利要求3所述的的分析仪，其特征在于，所述信号数据长度N是进行一次快速S变换的数据长度，根据谐波检测精度和计算机或嵌入式系统设备的运行速度综合考虑确定，由谐波次数h确定谐波分析采样频率f_s，

f_s≥2×h×50；

由采样频率和频谱分辨率B确定信号采样长度N，其具体公式为

N＝2^y≥f_s/B；式中y为正整数,5<y<20；

保证DSP计算时间小于一个工频电压周期，DSP计算时间为

t＝qN×cycle_FFT/f_core；

式中q为特征频率点数，cycle_FFT为DSP进行一次FFT的时钟周期，f_core为DSP内核主频。