CN101388001A - 基于全相位fft的高精度瞬间相位估计方法 - Google Patents

Abstract

本发明属于数字信号处理技术领域，涉及一种基于全相位FFT的高精度瞬间相位估计方法，包括下列步骤：对输入信号x(t)直接进行采样而得到离散数据序列{x(n－N+1)...，x(n)，...，x(n+N－1)}；将采样所得的2N－1个数据存储在DSP的RAM中，利用长度为2N－1个的卷积窗对所述的序列进行加窗，再将距离为N个采样间隔的数据两两叠加形成N个数据，对这N个数据进行快速傅立叶变换，输出N个复数，即对应为N个谱值Y(k)；进行谱峰搜索，找出N个谱值Y(k)中的幅值最大的谱Y(k*)，再将其虚部和实部进行比值，对该比值作反正切计算即可得时刻n的瞬间相位估计结果。

Description

基于全相位FFT的高精度瞬间相位估计方法

技术领域

本发明属于数字信号处理技术领域，具体涉及一种对各类工程领域信号的瞬间相位进行估计的方法。

背景技术

在通信、仪表、电力、光学应用、故障诊断等工程领域，存在大量对信号的相位进行高精度、高效、快速的估计问题。例如在采用QAM、PSK调制的通信系统中，只有在接收端对各个码元周期内的载波初相位进行快速估计，才能正确恢复星座图，实现正确解码^[1]；在各类仪表中，相位计是常见的一种，相位计的基本功能就是准确测算出信号各个时刻的相位^[2]；而在电力系统的并网合闸场合，要求精确地测量两电网系统中工频信号间的相位差；另外在电力系统谐波分析中，受发电机绕组工艺、各种非线性电力设备负荷(如各种晶闸管整流装置、变频装置)等因素的影响，电网频率通常会发生偏移，这就要求实时测出50Hz基波和各次谐波的相位、频率、幅值^[3]；相位式高精度激光测距仪就是通过测量连续调幅信号在待测距离上往、返传播所产生的相位延迟，来间接地测定距离的^[4]；在旋转机械的故障诊断中，常常需通过对振动信号进行谱分析(包括振幅谱和相位谱)来判断故障类别，然而有时一个振幅谱往往对应着两三个解释(比如某台从旋转机械采集到的信号频谱中存在着较大的2倍频分量，则表明这台机器可能存在不对中或者轴弯曲或者机械松动的问题，到底是哪一种情况，很难说清楚)，这时如果结合相位分析，就很容易作出准确的判断^[5]。总之，信号相位的估计问题普遍存在于与国民经济紧密相关的生产环节中，因而研发出高性能的相位估计法具有很高的工程价值和经济价值。

相位测量与估计的方法很多，传统依靠模拟器件的方法，如：矢量法、二极管鉴相法、脉冲计数法等，其测量系统复杂，需专用器件，硬件成本高。近年来，相位估计逐渐向数字化方向发展，其优点在于硬件成本低、适应性强，只需单片机、DSP、FPGA等通用器件就可完成，对不同的测量对象只需改变程序算法即可，且其精度一般高于模拟式测量。因此，选定一套精确的相位估计算法是关键。

然而，为研发出高精度的相位测量算法(即参数估计算法)，仅在工程领域中去考虑问题是远远不够的。因为无论是在哪种应用场合，也无论是测量哪种物理量，经采样后信号的表现形式都是离散观测数据，若要精确、快速、有效地从观测数据中提取出参数信息，这就要涉及很多的理论问题，所涵盖的知识可延伸到数字信号处理、信号检测与估计理论、信息论、概率论与数理统计、随机过程等多个学科领域，只有加深对这些领域的基础知识的理解，从新的角度提出一些优化参数估计性能的措施，才有可能开发出性能更为优良的算法。本发明即是在由申请人深入把握以上各个学科知识的基础上，加以灵活应用而形成。

而现有的数字化相位估计主要包括以下几种方法：

(1)脉冲计数法

这是最常用的相位测量法，其测相原理如图1所示，主要过程为：产生与所测正弦信号同频率f、且已知初相的参考正弦波，然后分别对这两路正弦波进行脉冲计数(假设其基本计数周期为T₀)，根据其脉冲计数差值(假设n个脉冲差)而算出其延时nT₀，进而测算出其相位差2 π f n T₀，由参考相位加上此相位差，即可得到相位估计。

可看出，这种方法的缺陷在于：(1)需产生参考正弦波，且要求参考信号与所测正弦信号的频率完全一致，若稍有偏离，则会给测量结果带来很严重的偏差；(2)需专门产生脉冲计数信号，且测量精度很大程度上取决于基本脉冲宽度T₀，只有减小T₀才能提高测量精度，这对硬件设备的要求很高；(3)仅能对单频信号进行测量估计。

(2)希尔伯特变换法^[6]

对于余弦类的信号x(t)＝cos(2π f t+θ₀)，为估计t＝0时的相位θ₀，以采样频率f_s对x(t)进行采样后得到序列x(n)＝cos(2π f/f_sn+θ₀)，n＝0，1，…，N-1；若对x(n)作希尔伯特变换可得其解析信号 $z\left(n\right)=x\left(n\right)+j\hat{x}\left(n\right),$ n＝0，1，…，N-1，则根据下式可得到的θ₀估计式

${\widehat{\mathrm{\θ}}}_0=\arctan \frac{\hat{x}\left(0\right)}{x\left(0\right)}---\left(1\right)$

然而希尔伯特变换法对采样频率f_s要求很高，当信号频率f＝kf_s/N，k∈z⁺时，上式测得的相位是准确的，一旦采样频率f_s稍发生偏离，使得这种整数倍的关系不成立，式(1)的估计就不正确；此方法的抗噪性能差；另外，此方法仍只能对单频信号进行估计。

(3)正弦曲线拟合法^[7~9]

John Kuffel提出了正弦曲线拟合法^[7，8]，该方法可获得很高的估计精度。它分为四参数(频率、幅度、相位和直流分量均未知)法和三参数法(除频率已知外，其它3个均未知)两种情况，文献[9]指出，四参数正弦曲线拟合过程并不是闭合的线性过程，尚无确切的数学公式可直接计算拟合参数，若拟合初始条件选择不当，易使得迭代过程发散或收敛到局部最优，且拟合需大量运算时间。为提高效率，文献[9]提出先估算频率，再进行三参数的正弦曲线拟合(为闭合线性过程，且绝对收敛)。然而，此方法要求频率估计非常精确，否则一个直接后果就是会把频率估计误差带入到相位估计中，另外，此方法仍只能对单频信号进行相位估计。

(4)相关法^[10~11]

文献[10-11]采用了相关法进行测相，假设参考采样序列为x(n)＝cos(2πf/f_sn)，n＝0，1，…，N-1，待测序列为y(n)＝cos(2πf/f_sn+θ₀)，n＝0，1，…，N-1，为估测出y(n)的初相，需分别求出参考采样序列x(n)和y(n)的自相关函数R_xx(τ)、R_yy(τ)，及其两者的互相关函数R_xy(τ)，根据式(2)可得到θ₀的估计

${\widehat{\mathrm{\θ}}}_0=\arccos \frac{R_{\mathrm{xy}}\left(0\right)}{\sqrt{R_{\mathrm{xx}}\left(0\right)R_{\mathrm{yy}}\left(0\right)}}---\left(2\right)$

然而从式(2)可看出，相关法测相仍需一路与代测信号频率完全一样的参考信号，且仍只适合于单频信号的相位估计；另外，实验证明，此方法的抗噪性能很差。

(5)各类基于FFT频谱校正的相位估计法

基于FFT的相位估计法是近年来工程界和学术界研究得最热的、也是实际应用最广的参数估计法。原因是，其一，由于FFT(Fast Fourier Transform，快速傅立叶变换)是DFT的快速算法，具有其它算法无法比拟的方便、简洁的计算优势；其二，在理论上，当N足够大时，若对离散采样值进行DFT可获得最大似然解的近似逼近；其三，各种丰富的信息可在FFT的结果中直接得到反映，如FFT谱线位置的分布反映了信号频率大小，而FFT结果为复数，因而隐含了信号的相位信息等(这是其它正交变换所不具备的，如DCT、沃尔什变换、K-L变换、Wigner-ville变换等)。虽然FFT的栅栏效应限制了其频率分辨率，但其栅栏效应可通过对FFT谱线进行内插的方法加以解决，这样就可衍生出多种离散频谱校正法，文献[12]总结了4类频谱校正法，包括能量重心法^[13]、比值法^[14-15]、DFT+FFT谱细化法^[16]、相位差法^[17-19]，这些方法都可对存在多种频率成份的信号进行较精确的相位估计，因而广泛应用于振动分析及故障检测、电力系统谐波分析、电介质损耗角测量、雷达测速等领域。

然而可通过实验和理论证明，各类基于FFT频谱校正的相位估计法存在如下缺陷：(1)FFT存在谱泄漏现象，当信号包含多种频率成份时，各种频率成份间会产生谱间干扰，这会很大程度地影响相位的测量精度；(2)由于存在谱泄漏，这些校正法尤其不适合于存在密集频谱的相位估计场合；(3)经调研发现，这些FFT频谱校正法在估计相位时，都是按先估计频率，得到频率偏离值后再进行相位估计的步骤进行的，因而这样做的一个直接影响是会把频率估计的误差带入到相位估计中。

本发明提出基于全相位FFT方法进行相位估计，全相位FFT是专利申请人王兆华在专利文献[22]提出的一种新型谱分析方法，谱分析包括振幅谱分析和相位谱分析两个方面，专利文献[22]只研究了全相位FFT的振幅谱特性，没有涉及相位谱的研究，本发明基于全相位FFT的相位谱性能，提出一种新型相位估计方法。

本发明的课题背景是国家自然科学基金与民航联合基金项目(No.60776815)和教育部博士点新教师专项科研基金项目(No.20070056105)。

发明内容

本发明的目的在于，克服现有相位估计技术存在的各种缺陷，提供一种高精度、低复杂度、快速、高效的相位估计方法。在各个工程领域，经采样获得信号的离散数据后，借助于本发明所提出的新型瞬间相位估计法对有限的离散数据进行处理，可以实现高精度、低复杂度、快速、高效地估计信号各个时刻的相位值；并结合具体的工程应用，将得出高精度的相位测量转化为实际工程中的其他物理量的测量值，使之为提高生产效率和生产力而服务。

为实现高精度、高效率、快速的相位估计，本发明提出如下的技术方案：

一种基于全相位FFT的高精度瞬间相位估计方法，包括下列步骤：

(1)对输入信号x(t)直接进行采样而得到离散数据序列{x(n-N+1)…，x(n)，…，x(n+N-1)}；

(2)将采样所得的2N-1个数据存储在DSP的RAM中，利用长度为2N-1个的卷积窗对所述的序列进行加窗，再将距离为N个采样间隔的数据两两叠加形成N个数据，对这N个数据进行快速傅立叶变换，输出N个复数，即对应为N个谱值Y(k)；

(3)进行谱峰搜索，找出N个谱值Y(k)中的幅值最大的谱Y(k*)，再将其虚部和实部进行比值，对该比值作反正切计算即可得时刻n的瞬间相位估计结果。

本发明提出的基于全相位FFT的瞬间相位估计法，应用于各工程领域，可产生如下有益的效果：

第一.可大大节省相位估计的硬件成本。这是因为

(1)apFFT测相法无需计数脉冲，也无需参考信号。

apFFT相位测量法不需要脉冲计数法所需的计数脉冲。

而且，脉冲计数法、相关函数法都需要另外产生与所测信号同频的、初相已知的参考信号，而apFFT测出的相位是样点的瞬间相位，实际上在FFT的计算矩阵中，矩阵的直流行元素对于其它行的元素来说，可起到隐含的信号参考作用，故不需要参考信号，节省了相应的信号电路成本。

(2)apFFT的“相位不变性”可大大地降低对A/D采样器件的要求。

不论是希尔伯特变换法，还是各种基于FFT的频谱校正法，都对A/D采样的采样率提出很高的要求：即要求所测信号的频率f恰好为FFT频率分辨率Δf＝f_s/N的整数倍，一旦采样频率f_s稍发生偏离(即为“不同步采样情况”)，使得这种整数倍关系不成立，就会大大地增大希尔伯特变换测相法的误差。同样地，这也会影响用FFT直接测相的精度，需借助附加的校正措施才可完成较精确的相位估计。

而apFFT因具有“相位不变性”，即所测相位与信号的频偏大小无关，故不对A/D采样的器件作苛刻要求。

(3)apFFT结构简单，本身即可大大降低相位估计的硬件电路成本。

由于图6的apFFT处理过程与图4的处理完全等效，N次FFT仅需1次FFT即可实现，故相应计算成本大大降低。

第二.apFFT测相具有很高的精度，尤其适合于精密测量场合。

前面已经证明，apFFT具有很高的相位估计精度；在信噪比很高时，相比于传统FFT测相来说，还具有数量级上的优势。故apFFT非常适合于精密测量场合(如激光测距仪就是通过精确测量往返信号的相位差，来精确测距)。

第三.apFFT测相效率高、速度快，尤其适合于实时快速测量场合。

apFFT测相延用了FFT的快速算法，且无需任何附加的频谱校正措施。此外对单频信号进行相位估计时，甚至还有无需FFT、仅需借助几次简单的乘累加操作即可实现精确测相位的简化结构。故apFFT测相的效率高、速度快，非常适于进行实时快速估计的场合(如用于在通信系统中，需通过快速估计载波相位来在接收端和发送端快速建立同步等)。

第四.apFFT尤其适合于信号包括多种频率成份和存在密集频谱的相位估计场合。

这是由于apFFT具有优良的抑制谱泄漏的性能所决定的。因为apFFT能够很好地抑制谱泄漏，当信号包含多种频率成份时，各频率成份的apFFT谱间干扰比FFT情况要小的多，因而非常适合于多频成份的测量，而脉冲计数法、希尔伯特变换法、正弦曲线拟合法、相关函数法都只适于单频信号测量。

另外，apFFT的谱泄漏范围比FFT情况要窄得多，故而当信号包含密集频谱成份时，谱泄漏相互覆盖的程度相比于FFT情况要低得多，故apFFT非常适合于存在密集频谱的相位估计场合，而文献[21]指出：“对于包含密集频率成份的离散信号，其校正方法目前还不是很成熟，现有的大部分校正方法都存在算法复杂、计算量大的缺点，很多方法只能校正仅含有两个频率成份的信号，因此改进现有的细化分析方法或者研究新的适用于校正密集频谱的方法将是未来研究的方向之一”——显然，apFFT测相法可填补此空白。

附图说明

图1脉冲计数法示意图。

图2全数字“瞬间相位计”硬件示意图。

图3传统加窗FFT谱分析(N＝4)示意图。

图4全相位FFT谱分析过程示意图。

图5各子分段DFT和全相位DFT的振幅谱。

图6本发明采用的等效的全相位FFT谱分析过程示意图(N＝4)。

图7传统FFT和全相位FFT的振幅谱和相位谱对照图。

图8本发明提出的基于全相位FFT的高精度瞬间相位估计方法的流程图。

具体实施方式

下面首先对实施本发明的硬件予以简单说明。参见图2，为精确估计出信号x(t)的瞬间相位值，需借助信号调理电路对输入信号进行模拟预处理，以对信号幅度范围进行必要调整，并去除外干扰噪声等；再经过A/D(模数转化器)采样得到样本序列x(n)以并行数字输入的形式进入DSP器件，经过DSP器件的内部的算法处理，而得到信号相位值的估计，最后借助输出驱动及其显示模块显示出相位估计值，即图2整个系统构成一个“全数字的瞬间相位计”。

其中图2的DSP(Digital Signal Processor，数字信号处理器)为核心器件，在瞬间相位估计过程中，完成如下主要功能：

(1)调用核心算法，完成瞬间相位的估计处理。

(2)根据实际需要调整采样率f_s，使之在该采样率条件下，尽量高精度地估计出瞬间相位。

(3)内部RAM存储数据不足时，将处理数据与外部RAM进行数据交换，以配合核心算法处理。

(4)将相位估计结果实时输出至驱动和显示模块，

需指出，由于采用了数字化的相位估计方法，因而决定图2系统的复杂度、实时程度和稳定度的主要因素并不是图2的DSP器件的外围连接，而是DSP内部程序存储器所存储的核心相位估计算法。若核心算法选用得当，则可使得整个硬件设备成本和软件开销都得到最大程度的节省，表现在：

(1)核心算法选用得当，则可降低DSP对A/D器件的采样频率f_s的控制复杂度，从而节省软件开销。

(2)核心算法选用得当，则可尽量节省处理过程中的数据存储空间，图2的外部RAM甚至都可省去。

(3)核心算法选用得当，则可大幅度地简化处理过程中的数学计算过程，实现高速、高精度的相位估计。

为此，本发明提出将基于“全相位FFT”的核心估计算法植入DSP器件内，基于此完成高精度、高速、高效的瞬间相位估计。相比于传统FFT，本专利提出的基于全相位FFT的核心算法在以下4个方面改进了相位估计性能。

第一.考虑包含某样点所有分段的FFT谱，以改善传统谱分析的谱泄漏问题。

传统FFT因具有谱泄漏的缺陷而降低了其参数估计精度。经典数字信号处理理论认为，由于点数为N的FFT只能对长为N的有限长序列进行谱分析，因而为产生有限长序列，就不得不对输入序列{…，x_-1，x₀，x₁，…，x_N-1，x_N…，}进行截断而产生序列{x₀，x₁，…，x_N-1}，但这种强制截断的直接结果就是会引起谱泄漏问题。虽然可对截断序列进行加窗来改善谱泄漏问题，即改为对加窗后的序列{x₀ f₀，x₁f₁，…，x_N-1 f_N-1}进行FFT谱分析(如图3所示，其延时器z^-1的延时量即为图2中A/D转换器的时钟周期T＝1/f_s)，然而这又会使得频谱的主瓣加宽，使之无法适合于密集频谱场合的参数估计。因而需根本善传统FFT的谱泄漏的问题。

为改善FFT的谱泄漏缺陷，仍需从解决截断问题入手，但又必须跳出对输入序列进行简单截断的思维定势。我们可看到，传统FFT所涉及的是整体序列的截断问题，倘若我们改从观察序列中的某个样点(如x₀)的角度来研究问题，则可看出，传统截断序列{x₀，x₁，…，x_N-1}仅仅考虑到了包含样点x₀的一种截断情况，若将所有包含x₀的截断情况全部考虑，则可形成如图4左边的各序列x_i(i＝0，1，…，N-1)，根据FFT的循环移位性质，相应地将各子序列x_i分别循环左移i个位置，则可形成图4各序列

(i＝0，1，…，N-1)，再将各

分别按照图4进行传统FFT(可加窗进行)，并将所有FFT分析结果叠加起来(或用窗序列b的值进行加权叠加)即可生成一种新的谱分析结果Y(k)。从图4可看出，样点x₀遍历了输入分段x_i的所有可能的位置，而在数字信号处理中，模拟连续信号被采样后，其相位只能通过离散的时刻位置来表示，故将图4的谱分析称为“全相位”FFT(all-phase FFT，简称apFFT)谱分析。

从图4可看出，全相位FFT蕴含了N路的子FFT的过程，故谱分析性能大为改善，而传统FFT仅是其中的第一路x₀的分析结果。

令N＝6，研究序列{x(n)＝2cos(2.2×2nπ/6+π/3)，-5≤n≤5}＝{-1.0000，-0.6180，1.8271，-1.8271，0.6180，1.0000，-1.9563，1.6180，-0.2091，-1.3383，2.0000}，对这11个数据按照图4处理：考虑所有包括中心样点x₀＝1.0000的长度为6的截断情况，并对各截断序列分别进行循环移位再进行DFT，即可得图5各子谱X_i(k)(如图5a~5f所示)，再对X_i(k)进行求和平均后得到全相位DFT谱Y(k)(如图5g所示)。

从图5可看出，x(n)两个共轭的频率成份(即x(n)＝e^{j(2.2×2nπ/6+π/3)}+e^{-j(2.2×2nπ/6+π/3)})在各个谱图中均得到反映(即k＝2和k＝4附近的两簇谱线)。相比于传统FFT(即图5a的X₀(k))，由于各子谱X_i(k)的泄漏在叠加中具有相互抵消的作用，使得图5g的全相位DFT大大抑制了谱泄漏：两共轭频率成份间的谱间干扰明显比图5a小得多，主谱更突出。因而apFFT非常适合于对存在密集频谱成份的信号作谱分析。

第二.构造与图4完全等价的、简化的“全相位FFT谱分析流程”。

文献[20]证明，图6的谱分析过程与图4完全等价。

而对比图4和图6所示的全相位FFT的两种等价结构可看出，图6相比于图4节省了(N-1)次的FFT计算，谱分析处理流程大大简化。故为提高谱分析的效率和节省硬件开销，实际工程应用时apFFT应采用图6过程进行处理。

从图6可看出，只需用长为(2N-1)的卷积窗w_c对中心样点x(0)前后(2N-1)个数据进行加权，然后将间隔为N的数据两两进行重叠相加可得到N个数据y(n)，将这2N-1个数据映射为N个数据的过程称为“全相位数据预处理过程”。对预处理后的N个数据y(n)进行FFT即得全相位谱分析结果。

对比图3和图6可看出，apFFT仅比传统FFT多出(N-1)次的乘累加运算，与DFT的快速算法FFT所需的N/2×log₂N次复数乘法相比，这部分附加的计算开销很小，但这却会换来性能的很大改善。

文献[20]指出，图6的卷积窗w_c为前窗和翻转的后窗的卷积，即

w_c(n)＝f(n)*b(-n)     -N+1≤n≤N-1                   (3)

令R_N为矩形窗，若f＝b＝R_N时，则称为无窗全相位FFT；若f、b中只有一个为R_N时，则称为单窗全相位FFT；若f＝b≠R_N时，则称为双窗全相位FFT。

第三.分别从各子分段FFT谱中提取相位信息并加以综合，以提取精确的相位谱信息。

以序列

为例，对该序列按照图4的apFFT各个子过程进行处理，可得到如下子DFT向量X_l，即

X₀＝[0.867 e^-j3°，1.286e^j27°，5.172e^j57°，2.258e^-j93°，1.041e^-j63°，0.819e^-j33°]^T

X₁＝[0.867e^-j141°，1.286e^-j51°，5172e^j39°，2.258e^-j51°，1.041e^j39°，0.819e^j129°]^T

X₂＝[0.867e^j81°，1.286e^-j129°，5.172e^j21°，2.258e^-j9°，1.041e^j141°，0.819e^-j69°]^T

X₃＝[0.867e^-j57°，1.286e^j153°，5.172e^j3°，2.258e^j33°，1.041e^-j117°，0.819e^j93°]^T

X₄＝[0.867e^j165°，1.286e^j75°，5.172e^-j15°，2.258e^j75°，1.041e^-j15°，0.819e^-j105°]^T

X₅＝[0.867e^j27°，1.286e^-j3°，5.172e^-j33°，2.258e^j117°，1.041e^j87°，0.819e^j57°]^T

将这N个子谱求和平均即得无窗全相位DFT的结果，即Y＝[0.1252e^j12°，0.2754e^j12°，4.4576e^j12°，0.8494e^j12°，0.1806e^j12°，0.1118e^j12°]^T注意Y(k)的各谱线相位值均为12°，即为输入中心样点x(0)的初相位θ₀＝12°。

文献[20]证明，序列 $\{x\left(n\right)=e^{j({\mathrm{\ω}}_{0}n+\mathrm{\θ})},-N+1\≤n\≤N-1\}$ 归一化后的不加窗情况的FFT谱分析结果X(k)和无窗apFFT谱分析结果Y(k)可表示成如式(4)所示(令ω₀＝β2π/N，β可以是小数)

$\left\{\begin{array}{c}X\left(k\right)=e^{j[\mathrm{\θ}+(1-\frac{1}{N})(\mathrm{\β}-k)\mathrm{\π}]}\·\frac{\sin \left[(\mathrm{\β}-k)\mathrm{\π}\right]}{\sin [(\mathrm{\β}-k)\mathrm{\π}/N]}\\ Y\left(k\right)=\frac{e^{\mathrm{j\θ}}{\sin }^2\left[(\mathrm{\β}-k)\mathrm{\π}\right]}{{\sin }^2[(\mathrm{\β}-k)\mathrm{\π}/N]}\end{array}\right.k=0,..,N-1---\left(4\right)$

再取式(4)的相位值，则可分别得到传统FFT相位谱φ_X(k)和apFFT相位谱φ_Y(k)

式(5)表明，传统FFT的相位φ_X(k)与信号的频偏值(β-k)密切相关，而全相位FFT则不同，其相位谱φ_Y(k)恒等于θ，与频偏值无关，称apFFT的这种性质为“相位不变性”，这包括三层含义：(a)apFFT测出的相位与信号的频率无关；(b)apFFT测出的相位为输入2N-1个数据的中心样点x(0)的瞬间相位，该相位接近理论值；(c)apFFT测相时不需要附加其它校正或补偿的计算，直接取峰值谱线的相位值即得样点的相位估计。

图7给出x(n)的加汉宁窗FFT和加汉宁双窗的apFFT的振幅、相位谱(N＝256)的对照图，x(n)是以f_s＝256Hz的采样频率对x(t)＝2cos(50.3×2π t+30°)+2cos(85.3333×2π t+60°)+2cos(121.444×2π t+90°)进行采样而得。

图7表明，FFT谱X(k)泄漏较严重，而apFFT谱Y(k)则不然，由各条主谱线泄漏出的旁谱线根数很少；且在各主谱k＝50、85、121处，传统FFT相位谱φ_X(k)很乱，而φ_Y(k)几乎等于中心样点x₀的理论值30°、60°、90°，即验证了apFFT的“相位不变性”。

为从数值上说明apFFT的高精度测相特性，以频率f_s＝3000Hz对x(t)＝cos(759.47×2π t+π/3)进行采样，分别采用无窗apFFT法和第一类相位差法(需对存在延时关系的两序列作FFT求取主谱线相位差，其延时量本文设置为N个采样间隔，故所用数据2N个，而apFFT只需2N-1个数据，详见文献[17])进行5000次蒙特卡洛测相仿真，并考虑到了量化噪声(分别将量化步长设置为m＝8和m＝10)。表1给出了不同高斯信噪比SNR条件下的测相方差数据。

表1 测相方差数据(单位：rad²)

表1实验数据表明，在比传统FFT测相法少用1个实验数据的情况下，apFFT测相法的均方误差仅为传统FFT相位差法的1/3~1/2左右；且信噪比SNR越大，优势越明显，信噪比进一步增大时，其精度可获得数量级意义上的提高。

第四.找出不同加窗方式下，相位谱估计的精确程度的差别，以优化基于全相位FFT谱分析的结构参数。

根据图6的卷积窗w_c的不同，全相位FFT谱分析分为无窗、单窗和双窗三种模式。文献[20]指出，相比于同阶数的传统FFT，这三种不同加窗模式的全相位FFT的抑制谱泄漏性能均有很大提高。但若对这3种加窗模式的apFFT进行对比，则抑制谱泄漏的程度双窗好于单窗，单窗好于无窗；而谱泄漏的宽度范围则恰好相反：无窗小于单窗，单窗小于双窗；另外文献[20]还指出，apFFT具有抑制噪声的能力，其抑制噪声的能力则为无窗强于单窗、单窗强于双窗。

因而可根据具体应用需求来选择合适的apFFT加窗模式；若信噪比很大时，建议采用双窗模式，反之若噪声干扰很大时，则采用无窗模式；若存在信号频率成份不是呈密集分布时，则选用双窗模式，反之选择无窗模式(因为谱泄漏范围小，谱间干扰不严重)；从而达到优化图6的卷积窗结构参数而提高相位估计精度的目的。

图8中，对输入信号进行A/D转换后(A/D转换可内置于DSP器件中，如DSP器件TMS320C2812)，将2N-1个数据存储在RAM中，调用存储于DSP内部程序存储器中的apFFT程序，该程序流程非常简单，具体处理过程如下：对输入数据进行apFFT得到N个谱值Y(k)，对这些谱值进行谱峰搜索，找出其幅值最大的谱Y(k^*)，再将其虚部和实部进行比值，对该比值作反正切计算即可得相位估计结果，最后将其输出显示。

图8流程中各处理阶段所涉及的具体硬件资源耗费情况，分别作如下阐述：

(1)对输入信号x(t)进行采样而得到离散数据x(n)。

这里对采样频率f_s值并不需要作严格限制，不要求信号频率为频率分辨率(由采样值f_s除以FFT长度N而决定)的整数倍。另外，对采样器的转换精度也没有提出过高要求，选用常用的字长数大于8bit以上的采样器即可。事实上，一般DSP器件所配备的最低的A/D转换精度就很高(如TMS320LF2407与TMS320F2812内部自带的A/D转换器为10bit的转换精度)，因而这足以保证apFFT在硬件实现时具有足够高的测相精度。

(2)内部数据的存储

可把apFFT的程序存放在DSP器件的程序存储器内，从图5可看出，apFFT需要2N-1个输入数据进行存储处理。需指出的是，由于apFFT具有很优良的抑制谱泄漏的性能，这保证了apFFT具有很高的测相精度；反过来考虑，这意味着在同样的测相精度的前提下，apFFT测相所需的样点数比传统FFT测相情况要少得多。可证明，在信号频偏量为0.1个频率分辨率，且信噪比为10dB的噪声环境下，若要达到0.0003rad²的测相方差精度(对应于均1度的绝对误差)，若采用文献[20]的传统DFT的时移相位差法，至少需1700个样本数据，而apFFT测相则仅仅需不到470个样本数据。就算以16位的采样精度进行采样，在每个数据需2个字(2Bytes)的存储空间，这意味着仅需940个字节(不到1kBytes)的内部RAM空间就可获得高精度相位估计，这对通用DSP芯片的内部RAM来说(如TMS320c2407的内部RAM为5kBytes，TMS320c2812的内部RAM为20kBytes，TMS320c5407的内部RAM为40kBytes，TMS320c6x的内部RAM为256kBytes)，足以存下这些数据；另外，从图5可看出，apFFT的主要处理还是在FFT环节，而FFT算法的各级蝶形运算都可实现原位运算，即在处理过程中，无需开辟新的数据存储空间进行算法处理。

从而在整个apFFT的实现过程中，内部RAM的数据空间就足够了，无需进行外部存储器的扩展。

(3)apFFT程序调用

而就apFFT测相程序本身而言，包括以下部分：(1)全相位数据预处理过程；(2)FFT处理；(3)谱峰搜索；(4)反正切；4个步骤组成，程序简单。进行数据运算时，将2N-1个离散数据序列x(n)存储在DSP的RAM中，利用长度为2N-1个的卷积窗对所述的序列进行加窗，再将距离为N个采样间隔的数据两两叠加形成N个数据，对这N个数据进行快速傅立叶变换，输出N个复数，这N个复数的幅值形成N个谱值Y(k)，然后进行谱峰搜索，找出其中幅值最大的谱Y(k*)，再将其虚部和实部进行比值，对该比值作反正切计算即可得相位估计结果。其中全相位数据预处理过程仅需对输入数据进行简单2N-1次乘累加(对应于DSP的MAC指令)操作；而FFT即为DFT的快速算法，其计算量从耗费N(N-1)次的乘法计算量减少到N/2·log₂N次，故复杂度得以大大简化；另外，由于FFT模值具有对称性，谱峰搜索只需对前N/2的FFT模值找出最大值，并记下其谱序号即可，而反正切操作仅仅是单个指令即得到瞬间相位值的精确估计。故整个apFFT测相程序非常简单，无需外部扩展程序存储器进行程序调用。

(4)计算结果输出

在由DSP硬件计算得到θ₀的估计值后，直接可通过DSP的输出总线输出至外部显示驱动设备进行数码显示。

需指出，由于采用了DSP实现，故测相操作变得更为灵活，可根据信号所包含的各种频率成份的具体情况，通过编程灵活改变apFFT的内部参数设置(如谱分析的阶数N、加窗模式等)。

【参考文献】

[1]丁康，罗江凯，谢明.离散频谱时移相位差校正法[J].应用数学和力学，2002，23(7)：729-735

[2]齐国清.利用FFT相位差校正信号频率和初相估计的误差分析[J].数据采集与处理，2003，18(1)：7-11

Claims (1)

1.一种基于全相位FFT的高精度瞬间相位估计方法，包括下列步骤：

(1)对输入信号x(t)直接进行采样而得到离散数据序列{x(n-N+1)…，x(n)，…，x(n+N-1)}；

(2)将采样所得的2N-1个数据存储在DSP的RAM中，利用长度为2N-1个的卷积窗对所述的序列进行加窗，再将距离为N个采样间隔的数据两两叠加形成N个数据，对这N个数据进行快速傅立叶变换，输出N个复数，即对应为N个谱值Y(k)；

(3)进行谱峰搜索，找出N个谱值Y(k)中的幅值最大的谱Y(k*)，再将其虚部和实部进行比值，对该比值作反正切计算即可得时刻n的瞬间相位估计结果。

Images