一种可变采样率的傅里叶分析方法与流程

本发明属于高频信号测试领域，特别是涉及一种可变采样率的傅里叶分析方法。

背景技术：

傅里叶分析最常见的即为离散傅里叶变换(dft)和快速傅里叶变换(fft)。dft是重要的谐波分析工具，可以针对复杂信号的采样序列进行数学变换，分离基波信号与各谐波信号。通常为了保证频率分辨率，需要增加采样序列长度n，当n值较大时，dft需要进行n2次复数乘法运算，所需的时间过长，若要保证良好的实时性则会对大幅增加对硬件的要求。

fft算法采用蝶形运算方式，能在较短的时间内实现谐波检测，但频率分辨率低，且要求同步采样和整数周期截断。若能保证精准的同步采样，fft针对谐波的测量可达到的精度很高；但对于间谐波的测量则取决于频率分辨率。通常，商业功率分析仪遵从iec标准而设置的窗长度与采样率只能满足(1～10)hz的频率分辨率，对于频率非整数的间谐波测量往往不够。非同步采样条件下，傅里叶方法固有的频谱泄漏和栅栏效应会大幅增加幅值和频率测量的误差，而在实际工程应用中实现严格的同步采样是难以实现的，因此如何减小频谱泄漏和栅栏效应是国内外学者研究的重点。

傅里叶方法的分析结果受频谱泄漏和栅栏效应影响较大，而且二者是相辅相成的，只有当被测频率分量恰好与频率轴单元重合时，才能得到精确的分析结果。通常，现有的谐波分析仪器和分析方法在低频领域会使用同步采样和栅栏效应的双重局限来实现该前提，但实则是为采样条件增加了强约束，缩小了采样率的有效作用空间，因此在高频领域的分析会受到采样率的限制。

为了使被测频率尽量与频率轴单元点重合，通常的谐波分析算法会要求满足同步采样和整周期截断两个条件，其中严格的同步采样是采样频率是所有频率分量的整数倍，否则一旦采样频率对某一次谐波形成非同步采样，则会对傅里叶分析的所有谱线结果产生影响。而实际工程中由于信号中包含谐波与间谐波，频率分量种类繁多，且各分量频率未知，很难实现严格的同步采样；同时整周期截断需要前期仪器的滤波来确定基波周期来计算截断时间窗的长度，在采取pwm或存在噪声的条件下，波形的过零点存在大幅振动，会导致基波周期的测量失准，而若加入滤波器的使用又会对主分量幅值产生一定的削减。

不发生频谱泄漏的充分必要条件即为被测频率与频率单元重合，同步采样与整周期截断都是为实现该前提而衍生出来的条件。对于工作频率达到几百甚至数千赫兹的高速电机，如果为实现同步采样求取各谐波次数的最小公倍数，则可能会超出目前硬件采样率的极限。

因此本发明专利针对这种情况，提出一种新型的可变采样率的傅里叶分解方法，从而提高傅里叶分解方法的精度，有效减小对于高频信号尤其是谐波信号的频谱泄露和栅栏效应。

技术实现要素：

本发明为解决高速电机的高频信号的分析采集与分析问题，提出一种可变采样率的傅里叶分析方法，从而提高傅里叶分解方法的精度，有效减小对于高频信号尤其是谐波信号的频谱泄露和栅栏效应。

本发明通过以下技术方案实现：一种可变采样率的傅里叶分析方法，所述傅里叶分析方法包括以下步骤：

s100通过对比初始设置的采样频率对所要分析的信号进行初步采样，进而采用傅里叶分析分析其基波频率及基波幅值；

s200通过分析得到的基波幅值进行初步判定，确定满足整周期截断的采样频率，重新对信号进行采样；

s300采样后的信号对其进行傅里叶分析，利用三谱线法确定最佳采样频率，即，既满足整周期截断，又满足无频谱泄露的最佳采样频率；

s400重新对信号进行采样，得到各次谐波的频率及幅值组成。

进一步的，在步骤s100中，具体的，通过公式f＝np/60直接计算运行频率f，以基频估计值的2n倍为初始采样频率进行采样，

进一步的，所述步骤s100前还包括以下步骤：

步骤s000获得被测高速电机反馈的转速信号。

进一步的，所述最佳采样频率具有多个数值，且呈现周期性变化。

进一步的，所述最佳采样频率的周期性变化，随频率升高而周期逐渐增加。

本发明的有益效果在于：本发明设计了一种可变采样率的傅里叶分析方法，能够根据转速信号快速估计初始采样率，随后根据计算算法能够快速确定最优化采样率，能够有效减小高频信号的采集误差及分析误差，能够将栅栏效应和频谱泄露减小至0，有效提高傅里叶分析算法的分析精度。同时，采用该种算法可以采用较小的采样率实现高精度、高频信号的采集及分析，减小了采集系统的成本及加快了分析速度。

附图说明

图1为可变采样率fft分解的计算流程图；

图2为最优化采样率搜索流程图；

图3为采样频率变化与基波幅值与频率的关系示意图；

图4为采样频率的变化与基频谱线之间的关系；

图5为采样频率变化时所呈现的周期性规律；

图6为主旁瓣分布的一种情况的示意图；

图7为主旁瓣分布的另一种情况的示意图；

图8为基于二分法的搜索原理示意图。

具体实施方式

下面将结合本发明实施例中的附图对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

参照图1所示，本发明提供了一种可变采样率的傅里叶分析方法，傅里叶分析方法包括以下步骤：

s100通过对比初始设置的采样频率对所要分析的信号进行初步采样，进而采用傅里叶分析分析其基波频率及基波幅值；

s200通过分析得到的基波幅值进行初步判定，确定满足整周期截断的采样频率，重新对信号进行采样；

s300采样后的信号对其进行傅里叶分析，利用三谱线法确定最佳采样频率，即，既满足整周期截断，又满足无频谱泄露的最佳采样频率；

s400重新对信号进行采样，得到各次谐波的频率及幅值组成。

在本部分优选实施例中，在步骤s100中，具体的，通过公式f＝np/60直接计算运行频率f，以基频估计值的2n倍为初始采样频率进行采样，

在本部分优选实施例中，所述步骤s100前还包括以下步骤：

步骤s000获得被测高速电机反馈的转速信号。

在本部分优选实施例中，在步骤s300中，采样频率根据所需要判别的谐波频率及谐波幅值对象进行最优化调整，从而确定最佳采样频率。

在本部分优选实施例中，最佳采样频率具有多个数值，且呈现周期性变化。

在本部分优选实施例中，最佳采样频率的周期性变化，随频率升高而周期逐渐增加。

具体的，参照图1所示，根据被测高速电机反馈的转速信号，由公式f＝np/60直接计算运行频率f，为了能够同时保证整周期截断和序列长度的要求，可以以基频估计值的2n倍为初始采样频率进行采样，同时最终的采样序列长度也要控制为n＝2p(p>n)，即信号周期数为2p-n，记作np；而对于m次谐波，在序列中的周期数为m×np。其中需要注意的是，基波在一个周期中被采样了2n个点，则m次谐波一个周期只被采样了2n/m个点，当谐波次数较大时，可能会无法满足采样定理，因此在设定采样频率时，需要根据电机特点对n进行适当调节，保证满足下式：

式中，mmax为所需分析的最高谐波次数。

若则实际基频应该对应第2^p-n条谱线，则此时m次谐波的频率则对应第m×2^p-n条谱线，即满足了与频率单元重合的条件。但由于转速反馈的估计值存在误差，一定会出现频谱泄漏结果，为了最大限度地减少频谱泄漏，本方法通过改变采样率来获得不同的信号样本序列，利用这些样本序列进行傅里叶分析，并寻找幅值结果的最大值。当幅值达到最大值时，频谱泄漏几乎全部消除，而经过不断修正的采样率，也逐渐逼近了最优采样率fsop＝f1×2^p-n。

结合图1来说明，图1中三谱线分析是用来取代fft的低计算量方法，由前文的分析可知，仅需针对基频实现优化，即可最大程度完成所有谐波分量的频谱泄漏的消除。因此，在每次更改采样率之后，仅需关注基频分量的幅值与频率情况，而无需进行fft完成全频域的计算。三谱线分析会基于之前获得的基频信息来估计新分析结果中的谱线位置。

式中[]取整数，可以四舍五入。求得对应谱线的幅值a^*(k)，其计算方法如下所示：

然而与实际基频f1的差异可能会导致取整函数的结果出现差错，使得对应谱线位置的次序相差1，即恰好计算的是主瓣两侧的谱线。为了避免出现这种差错，通常会计算第k条谱线左右相邻的谱线幅值，并找到三者中的最大值作为基频幅值的估计值。后续步骤中若再次改变采样率时，计算的谱线位置未必出错，因此并不是每次都需要计算三条谱线的幅值。为此，设置一个三谱线分析的开启阈值：

式中at为开启阈值，当计算得到的谱线幅值小于at，即其旁瓣与主瓣几乎相等，则说明该谱线一定不是主瓣，需要另行计算寻找主瓣。而主瓣对应的频率单元即为新得到的基频分析结果:

根据新得到的基频值f1^*，按照之前的设定定义新的采样率2ⁿ×f1^*。

参照图3所示，图3以8000hz的采样频率采样一组基波频率为500hz，幅值为10v；含有三次谐波(频率1500hz,幅值为3.3v)、五次谐波(频率为2500hz，1.7v)；取n＝4，p＝9，对其进行分析的结果。图3可见各频率分量的幅值都会随着采样率的变化展现出一种类抛物线的波动规律，且谐波的波动周期明显较基波小很多；且当基波取得峰值时，各谐波也一定可以达到峰值，反之则不成立。若继续扩大采样率的变化范围，可以发现幅值的波动周期会随着采样率的增加而逐渐增大。图4说明利用8khz采样基频为500hz的信号，序列长度为512对应的是32个周期，此时基频对应的为频率轴上的第32根谱线；而其相邻的峰值由于也同样消除了频谱泄漏的影响，必定也使得基频与频率单元重合，而由于采样率的改变，频率单元也会发生相应的变化，故左右两侧的峰值对应的基频恰好为第33条谱线和第31条谱线。同理，基于以上论述，我们也可以得知,与8khz相邻的左侧波谷值所对应的基频，应当恰好处第32与第33条谱线中央，而与其相邻的左侧波谷值所对应的基频，应当恰好处第31与第32条谱线中央。图5说明了不同的采样频率所处的位置。

参照图2、图6和图7所示，图6中左侧谱线幅值大于右侧谱线幅值，说明基频的真实值在主瓣与左侧旁瓣之间；同理，图7中情况对应的则是基频的真实值在主瓣与右侧旁瓣之间，则通过这种方式，可以判断出2ⁿ×f1^*与fsop的方位关系。除此之外，主瓣与旁瓣的幅值还能用来进一步修正得到的基频测量值，对于图6中的情况，修正判断为：

对于图7中的情况，则变更频率为：

式中f1^**为经过修正的基频值，a^*(k)，a^*(k-1)，a^*(k+1)分别为主瓣和左右侧旁瓣的谱线幅值，λl，λr分别为向左侧和向右侧修正采样率时的修正系数。

参照图2和图8所示，经过三谱线法的修正后，f1^**和基频真实值间的偏差已经很小了，但若继续使用三谱线法优化采样率，最后极有可能会不收敛，因为三谱线法的优化较为粗略，幅度较大，会导致最终优化值在实际值周围震荡，但无法收敛。因此，考虑使用二分法来快速锁定最佳采样率的精确位置。

二分法是一种适合在大数据量区间中搜索的快速方法，借助二分区域的原理可以以指数级速度缩小区域，如示意图8所示，其计算次数n0取决于搜索区间l＝[fa，fb]和计算精度e：

式中的取整函数需要正向取整。由于经过两次修正后的采样率与最优采样率极为接近，因此可以大幅压缩二分法的搜索区间，从而减少搜索次数。如前，经过修正的2ⁿ×f1^**必定分布在fs^c与fs^d之间，因此可以通过使fs^c与fs^d二者相减，计算搜索区间。但该搜索区间过大，而2ⁿ×f1^**与最佳采样率fsop的具体差值未知，不宜盲目缩小区间，而利用三谱线修正基频的原理可以进一步缩小搜索区间。再次使用三谱线法分析以2ⁿ×f1^**采样得到的新序列，同样按照上述步骤，根据旁瓣幅值判断采样率修正的方向，然后利用主瓣、旁瓣的幅值比例求取基频的修正系数λ：

式中a**(k)，a**(k-1)，a**(k+1)为再次三谱线分析的主瓣、旁瓣幅值结果。为确保最优采样率fsop在区间中，将基频的修正量扩大至两倍。

则对应的另一区间端点为：

搜索区间确定好后，根据二分法原理取两端点采样率和的二分之一为新采样率，并利用三谱线方法计算基于新采样率的新信号序列的主瓣幅值aop，并计算aop与两端点采样率对应的主瓣幅值的差值为：

δa＝aop-aa(16)

δb＝aop-ab(17)

根据δa和δb的大小，可以确定更新区间的端点，并重新赋值：

设置计算精度为10^-4级，当新区间两端点对应采样率的主瓣幅值结果间的差值δab＝|aa-ab|小于该精度时，即可认为二分法搜索完成。