CN1023156C - 快速高精度信号频谱分析方法 - Google Patents

Abstract

本发明在严格的理论分析基础上，提出了一种新的快速高精度信号频谱分析方法。它主要包括下列步骤：初始数据采集；奇偶排列；按现有FFT方法计算even(k)和′odd(k)；利用事先存储的实常数Uk.Rk和Zk合成Ak；将K次谐波的系数Ak显示或打印。本发明不仅可提高频谱分析精度，扩展频谱分析范围，缩短分析时间，还可降低硬件性能要求与简化硬件结构。

Description

本发明是关于一种对被测电信号进行频谱分析的方法，更具体地说它是关于一种具有抗混效果的快速富里叶频谱分析方法。

当前，数字计算机与数字信号处理器（DSP）等的迅速发展与普及，使得模拟信号也被离散成数字信号来处理。然而从信号的复原与信号的频谱分析角度考虑，模拟信号用离散信号代替后会出现以下原理误差。其一是非时限信号用有限时域的离散信号代替时产生的截尾误差（Truncatederrors），对周期信号这种误差又表现为采样周期与信号同期不同步时出现的泄漏误差（Leakage    errors）;其二是非带限信号时，以及带限信号但采样频率低于香农（Shannon）采样定理要求时出现的频谱混迭误差（Aliasing    errors）;其三是离散信号频谱周期性重复，不管原模拟信号包含多少频率分量，采样N点均只能确定N/2个带有上述误差的频谱系数，并且频率越高的分量，频谱系数的误差越大。

泄漏误差已找到一些较有效的解决办法，比如加窗技术与加窗插值技术等。混迭误差虽然已深入研究80余年，在各种介绍FFT原理的文献中都被涉及到，但还只停留在如何确定误差规模阶段;如何消除这种误差，特别是如何在信号快速处理过程中消除这种误差尚待解决。

频谱分析是信号处理的核心，电子系统的谐波分析即频谱分析之一种。分析连续信号频谱时，目前普遍采用快速富里叶变换（FFT）。现有各种版本的FFT都只是离散富里叶变换（DFT）的一种快速算法，属数字处理技术之列，存在着本发明要讨论的频谱混迭误差与频谱系数周期性重复的上述通病。混迭误差严重程度随信号形式而异。以一周期内的表达式为

x（t）＝（2/T）t-1

的锯齿波为例（见附表1），在采样周期与信号周期同步因而没有泄漏误差的情况下，一周期采样64点时，现有FFT算出的第31次谐波的正弦系数只有理论值的0.074748，一周期采样128点时，31次谐波的正弦系数也只有理论值的0.079914;余弦系数（理论值全为零）误差更大。混迭误差的影响由此可见一斑。

为了消除混迭误差，现有频谱分析仪采取的对策是，用低通滤波器将信号中高于1/2采样频率的高频分量滤除掉。这种方法虽能减少，所分析出的N/2个频谱系数中的混迭误差，但却带来其它误差。原因是低通滤波器的幅频特性不可能是矩形，相频特性不可能是直线;信号通过后，高频分量滤掉了，低频分量的幅度与相位却产生了失真，使得分析结果仍与真实情况不同。特别是低通滤波器带来的相位误差，在许多要求精确知道谐波相位的场合是不允许的。因此采用低通滤波器也不能很好解决现有FFT算法存在的精度不高，所分析的频谱范围不宽的问题。

面对这一情况，一种保留信号的高频分量仍能压制谱混迭误差，从而可大大提高频谱分析精度与扩大所分析的频谱范围，并能缩短分析时间的高精 度快速频谱分析方法，不仅对工程应用有重大意义，对数字处理技术与信号分析理论也有重大意义。

本发明的目的就是要提供一种上述的快速高精度频谱分析方法。

下面阐述一下本发明快速高精度频谱分析方法的数字原理。

满足Dirichlet条件的函数x（t），在[o，T]区间上可用Fourier级数表示为

A_K称为K次谐波的复振幅，A_K是其幅值，ψ_K是其相位，a_K，b_K分别为K次谐波的余弦系数与正弦系数。所谓频谱分析，就是求出各次谐波对应的A_K。

设△为很小的时间间隔，T＝N，△，x（t）在t_n＝n△时的值为x_n。用各个X_n代替上式中的x（t），即用求和代替积分时，就得到了现有快速频谱分析方法FFT计算A_K的公式：

采样定理从信号理论角度分析了A_KF与A_K的差别，从数学角度考虑，A_KF是求和得到的，真正的A_K应由积分得到，故用FFT方法得到的信号频谱与真实频谱之间必存在原理误差。

为求出真正的A_K，不能采用现有FFT采用的简单求和方式，本发明采用精确计算（2）式所给的积分。

A_K计算可分段进行：

一种可供选择的计算A_kn的方法是利用欧拉数值积分公式，参见安德列、安戈著、陆志刚等译，人民邮电出版社出版的“电工、电信工程师数学（下）”，第十章。该公式计算积分时多次实行分部积分法，直至后面的分部积分结果可忽略。此法直接运用于（5）式时效果不佳，原因是被积函数中有周期因子exp（-jzkπt/T），无论分部积分多少次，这个因子总原样存在，无法判定后面部分是否可忽略。针对（5）式这样一种特殊被积函数，我们对欧拉数值积分方法作了改进。我们先将被积函数中的x_n（t）部分在[t_n-△，t_n+△]区间展开成台劳级数，由于△很小，x_n（t）可用它在该区间展开式的前m+1项来逼近：

（6）

M的取值视△的取值而定，△大则逼近x_n（t）所需的M值高，△小则M值低。将（6）式代入（5）式，这时再多次使用分部积分法，直至被积函数变为零为止，结果得：

（7）

式中x^（n21） _n（t），x^（21+1）（t）分别为x_n（t）的21阶与（21+1）阶导数。现有FFT的频谱系数计算方法，可看成是（6）式在M＝0时的特例，即用x（t）在t_n处的采样值代替[T_n-△，t_n+△]区间的x_n（t），因而带来混迭误差。为消除混迭误差，必须考虑x_n（t）展开式中的高式次项。一般△＜＜T，因此取x_n（t）在该区间的级数展开式的前三项即可有相当高的精度。将（6）式中的各个系数a_m用t（n-1）△，n△，（n+1）△时的采样值x_n-1，X_n，X_n+1的函数来表示，并考虑截尾误差，然后代入（7）式，计算整理得：

（8）

式中，

Z_k＝U_k+jV_k

Z_k为Z_k的共轭复数，V_k、U_k与R_k则是只与比值K/N有关的实常数，其值随K的增加而单调下降。将（8）式代入（4）式，最终得：

式中

$W_N\mathrm{=e}{}^{\mathrm{-J}\frac{\mathrm{2x}}{N}}$

F＝（X₀-X_N）/N （12）

（10）式就是本文高精度频谱分析的基本公式。

（10）式中上述的实常数U_k、V_k、R_k的数学式如下：

（10）式所示的具有抗混（Anti-Aliasing）效果的A_k精确表达式如不能实现快速计算，本文方法将难以推广应用;用按（10）式快速计算A_k时，如能尽可能利用现有FFT程序，则本文方法更易于推广普及。下面我们来寻找满足这种要求的算法，并将这种算法命名为FAFT（Fast Anti-Aliasing Fourier Transform）。

（10）式中的U_k、T_k、V_k都是与采样值无关的常数，可预先算好储存备用;F也只与首尾两点的采样值有关，容易计算。关键是如何实现

于是

_k为：

_k＝2U_kAeven（K）+R_kAodd（K）W_m-Z_kF （18）

这时如能同时实现Aeven（K）与Aodd（K）的快速计算就等于实现了A_k的快速计算。由（16），（17）两式可知，Aeven（K）与Aodd（K）分别为N/2个偶次采样数据与奇次采样数据对应的复频谱系数，为计算它们须先将采样数据按奇、偶顺序重新排列。这并非本文方法的额外要求，现有FFT也须如此重新排列数据：时域抽取法（Decimation in Time）开始时进行，频域抽取法（Decimation in Frequency）结尾时进行。所以这一步可直接利用现有FFT的有关程序实现。

计算

even（K）与

odd（K）也可借用现有FFT程序，比如时域抽取法的Cooley-Tukey算法。仔细分析这种算法可知，对于N＝2^S个被分析数据，该算法共循环S步，其中第（S-1）步变换出的数据，正好是这里的 even（k）与

odd（K）。所以利用现有FFT算法即可实现

even（K）与 odd（k）的快速计算。第（S-1）步后，将现有FFT程序稍加改动，即不继续计算Aeven（K）W^K _N与

odd（K）W_N，而是只计算后者，然后按（15）式要求合成A^A _k。这样，借用几乎没作什么改动的现有FFT程序，就可一举实现高精度A_k的快速计算。

下面将结合附图对本发明进行详细说明。

图1是实现本发明方法的频谱分析仪的方框示意图;

图2是本发明信号频谱分析方法的流程图;

图3是第一实验的信号波形;

图4是第二实验的信号波形。

为了便于理解，首先对实现本发明方法的频谱分析仪的方框图作一介绍（参见图1）。它主要包括数据采集电路（1）、存储器（2）、中央处理单元CPU（3）、显示器（4）、打印机（5）以及键盘（6）等。数据采集电路首先将要分析的信号（模拟形式）拾取，采样并转换成数字信号，该数字信号通过总线（7）送入存储器（2），中央处理单元CPU（3）将存储器（2）中的数据进行处理并计算，最终得出各次谐波分量的系数A_k值，并显示在显示器上或通过打印机打印出来。

下面将结合图1和图2对本发明的频谱分析方法作一详细介绍。

按照本发明的频谱分析方法，它包括如下步骤：

1、初始数据采集

数据采集电路（1）对输入的被测信号进行N个点的采样，并将采样点的值变成数字信号存入存储器（2）中。设该N个采样值为：

X₀，X₁，X₂……X_N

其中采样点数N可以根据测量的精度要求而下，比如：N为4、8、16……，实际操作时可以通过键盘来设定N的值。N值定了以后，按照 N＝2^S则S的值自然就定了，S为现有FFT中的计算步骤。

根据所存储的采样值计算吉卜斯现象因子F：

F＝（X_O-X_N）/N

2、利用现有的快速富里叶变换FFT程序作奇偶排列

即将存储在存储器（2）中的初始采样数据按奇数和偶数重新排列，以便于以后的计算，

X₀，X₂，X₄……X_N-2

X₁，X₃，X₅……X_N-1

3、执行现有快速富里叶变换FFT方法前（S-1）步，计算：

even（K）、

odd（K）分别为偶次采样数据与奇次数据对应的复频谱系数。

4、完成现有FFT方法第S步的一半计算内容，求出：

′odd（K）＝ dod（K）W_N

A′odd是旋转加数后的奇次采样数据的复频谱系数，W^K _N是旋转因了。

5、合成A_k

按照下式求出K次谐波的复振幅

_k：

_k＝2U_kAeven（K）+R_kA′odd（K）-Z_kF

其中实常数U_k、R_k、Z_k是事先分别按照公式（13）、（15）和（9）计算好并存在存储器中的。

6、显示或打印

将上面求出来的K次谐波系数A_k显示在显示器（4）上或通过打印机（5）打印出来。

本发明的快速高精度信号频谱分析方法（FAFT）用三项多项式逼近各区段的X₃（t），精确算出了A_kn对应的积分值，精度必然比FFT高;这样，采样间隔取大一些仍能比FFT精度高。对于同一个被采样信号，采样间隔大就意味着采样点数N少，而分析过程中所需的复数乘、加运算次数M与采样点数N之间有以下关系：M＝Nlog₂N，N小即意味着计算时间短，因此FAFT与FFT相比，不仅可大大提高精度，还可缩短计算时间。另外，FAFT不是用离散信号频谱代替连续信号频谱，不受采样定理限制。对带限信号，即使用低于采样定理要求的速率采样，仍能有很高的精度。在硬件采样速率已定的前题下，这等于扩大了可分析的频谱范围。而在精度要求一定的情况下，与FFT相比，FAFT可采用采样速率低的采样电路与位数低好A/D变换器，这就可降低硬件成本。

FFT的另一缺点，即采样N点只能确定N/2个频谱系数的缺点，是这一方法算出的

_k+N＝

_k这一周期性重复性质造成的。FAFT无此弊病。

由

even（K），

odd（K）与W^K _N的递推性质：

even（K+mN/2）＝ even（K） （19）

odd（K+mN/2）＝

even（K） （20）

W^N _N+mN＝W^N _K（21）

W^N _N+（m+1/2）N＝-W^K _N（22）

可得FAFT法算出的A_k的递推公式为：

由于U_k，R_k，Z_k无周期重复性质，故FAFT法算出的A_k+mN/2与A_k之间不存在周期性重复关系，即利用N个采样值可算出N/2个以上的频谱系数。

FAFT方法可直接利用现有信号分析仪实现，也可另行设计更能充分利用FAFT优点的硬件来实现。与现有信号分析仪相比，采用FAFT后不仅可降低对采样速率与A/D精度的要求，还可省去数据采集电路中的抗混滤波器及有关部分，使硬件成本降低，结构简化。

为了证实FAFT的上述优点，对已知频谱系的两种常见波形-锯齿波X₁（t）与余弦全波整流波形X₂（t），分别用FAFT与标准FFT（10）作了频谱分析。X₁（t），的波形分别见图3、4，它们在一周期内的表达式与理论频谱系数分别为：

式中，

b_k＝+ 1/(kπ) （26）

（27）

式中，

a₀＝ 2/(π) （28）

a_2k=(-1)^k-12/(π(4K²-1)) (29）

分析结果见表1-表4。表中f₁＝1/T为信号基波频率，f_s＝1/△为采样频率，K为谐波次数，f_k＝Kf₁为K次谐波频率，a_k，b_k为分析值，a_k，b_k为理论值。为节省篇幅，每种比较只列出了一种波形的结果。

表1是同样的采样速率与A/D精度时两种方法分析精度对比。可看出，不仅f_s相同时FAFT的精度远比FFT的精度高，即使FAFT的f_sFFT的f_s低时，FAFT的精度也远比FFT的精度高。

表2是同样的A/D精度，不同的采样速率时，两种方法的分析精度对比。可看出，对波形1即使FAFT的f_s有FFT的f_s的1/64，分析出的前31次谐波系数的精度仍比FFT的精度高。

表3是同样的分析精度时两种方法所需的采样速度与A/D精度对比。可看出，对波形2在分析出的前62次谐波系数精度基本相同的情况下，FFT需采用每周期采样256点的采样电路与18位的A/D时，FAFT只需采用每周期采样128点的采样电路与10位的A/D。

表4是不同的采样速度时，FAFT分析出的f_k＝50f_s、120f_s、250fs的K次谐波的谐波系数及精度。这不仅说明FAFT可扩大谐波分析范围，还说明此法的确不受采样定理限制。按采样定理，欲分析第3998次谐波，采样频率至少应为基频的7996倍，即f_s≥7996f₁;FAFT方法，f_s＝16f₁即可以相当高的精度分析出3998次谐波及更高次数谐波的谐波系数。本例还证明与FFT相比，FAFT可大大提高分析速度。比如欲分析0～3998次谐波，FFT需处理7996个以上的采样数据，FAFT只处理（16⁺1）个数据即可;这时两种方法所需计算时间长短不言自明。

Claims (1)

1、一种快速高精度信号频谱分析方法，它包括下列步骤：

(1)初始数据采集：即通过数据采集电路对输入的被测信号进行N个点的采样：X₀，X₁，X₂……X_N，并将它们变成数字信号存入存储器中，再计算吉卜斯现象因子F值，F=(X_O-X_x)/N，

(2)利用现有的快速富里叶变换FFT程序对上述N个采样值作奇偶排列，即：

X₀，X₂，X₄……X_N-2

X₁，X₃，X₅……X_N-1

(3)按照现有快速富里叶变换FFT方法前(S-1)步计算偶次采样数据的复频谱系数

even(k)：

(4)进行现有快速富里叶变换FFT方法第S步的一半计算内容，求出旋转加数后的奇次采样数据的复频谱系数

odd(k)：

odd(k)=Aodd(k)WW^K _N

其中W^K _N是旋转因子，

其特征在于该方法还包括下述步骤：

(5)利用事先存储的实常数U_K、R_K和Z_K按下式合成K次谐波的复振幅 _K：

=2U_K

even(k)+R_K

odd(k)-Z_kF

(6)显示或打印K次谐波的复振幅 数据。

Images