CN104268123B

CN104268123B - 一种离散数字信号任意步长的滑动离散傅里叶变换方法

Info

Publication number: CN104268123B
Application number: CN201410492336.4A
Authority: CN
Inventors: 阎啸; 王茜; 秦开宇
Original assignee: University of Electronic Science and Technology of China
Current assignee: University of Electronic Science and Technology of China
Priority date: 2014-09-23
Filing date: 2014-09-23
Publication date: 2017-05-17
Anticipated expiration: 2034-09-23
Also published as: CN104268123A

Abstract

本发明针对HDFT变换方法在实际工程应用中HDFT变换存在潜在不稳定和累计误差的问题，提出了一种数字离散信号任意步长的滑动离散傅里叶变换方法即调制任意步长SDFT(Modulated Hopping DFT，mHDFT)变换方法，利用DFT循环频域移位特性，采用修正的UVT变换，有效地去除了HDFT谐振器反馈回路中的旋转因子，将HDFT极点精确固定在单位圆上，然后通过相位修正得到n时刻频点k的DFT变换结果，从而在确保HDFT恒稳定的同时有效降低了累计误差。

Description

一种离散数字信号任意步长的滑动离散傅里叶变换方法

技术领域

本发明属于离散数字信号处理技术领域，更为具体地讲，涉及一种离散数字信号任意步长的滑动离散傅里叶变换(Sliding Discrete Fourier Transform，SDFT)方法。

背景技术

离散傅里叶变换(DFT)作为信号时域与频域转换的桥梁，是现代数字信号处理理论的基石之一。它被广泛应用于包括高精度实时频谱分析、认知无线电中的频谱感知、宽带无线通信中正交频分复用信号处理等在内的众多信号处理领域中。为了提高DFT计算效率，在实际工程应用中通常采用快速傅里叶变换(Fast Fourier Transform，FFT)算法或Goertzel算法来实现DFT变换。

FFT算法利用了DFT算法中旋转因子的对称性、周期性等特点，通过不断将长序列分解成更短的序列，有效减少乘加次数以及旋转因子的数目，降低算法运算量。在利用FFT算法进行M点DFT运算时，通常要求采用M个信号时域连续样点构成数据帧，利用FFT对其进行成块变换。这就造成FFT算法每隔M个信号等效采样周期输出一帧由M个频点构成的信号频谱。基于二阶无限冲激响应滤波器(IIR)结构的Goertzel算法虽然能够通过迭代运算获得低于FFT算法计算量，但该算法只能在间隔M个信号等效采样周期的特定时刻输出对应频点的频谱信息。在某些强实时信号处理应用中,要求DFT变换能够在每个时钟周期产生并输出部分给定频点频谱信息，传统的FFT算法和Goertzel算法都无法满足这一实时性要求。虽然可以利用帧重叠FFT(Overlap FFT)技术，在每个输入信号采样样点对其进行全部M点DFT变换，并同时提高时-频分析光谱图的时间分辨率，但这将使算法运算量大大提高导致工程上难以实现。

另一方面，在某些实时性相对要求不高的特殊信号处理应用中，例如空间电磁环境实时监测、协作通信中的信号同步处理和盲信号特征参数提取等，并不要求在每信号输入样点时刻都计算产生所有频点DFT输出，只需在每间隔L(0<L<M)个时钟周期更新输出部分选定频点的频域信息。传统的FFT算法虽然可以利用Overlap FFT技术实现每间隔L(0<L<M)个时钟周期输出全部M点频谱信息，但是无法单独计算输出部分选定频点的频域信息，计算效率极低；Goertzel算法能够通过单个频点的递归运算实现部分选定频点的频域信息计算输出，但其频谱信息更新时间间隔必须为M个时钟周期。

滑动离散傅里叶变换(Sliding Discrete Fourier Transform，SDFT)是一种基高效连续递归DFT算法，其输出频点数据速率能够接近或等于输入信号时域样点速率，被广泛应用于宽带无线通信、自动控制和航天航空等强实时数字信号处理领域中。

但是在逐样点(Sample-by-Sample)运算机制的任意步长SDFT(Hopping DFT，HDFT)中，每次迭代运算中都会引入的量化误差，在运算结果中引入随迭代次数增加的累积误差，因此，在实际工程应用中HDFT算法存在潜在的不稳定性。

发明内容

本发明的目的在于克服现有滑动离散傅里叶变换的不足，提供一种离散数字信号任意步长的滑动离散傅里叶变换方法，以提高离散数字信号傅里叶变换的精度和稳定性。

为实现以上目的，本发明离散数字信号任意步长的滑动离散傅里叶变换方法，其特征在于，包括以下步骤：

(1)、将输入信号序列的信号样点x(n)减去M个采样间隔之前的信号样点x(n-M)，得到差值信号d(n)，即：

d(n)＝x(n)-x(n-M)，

其中，M是DFT变换点数，n信号样点的时域索引值；

(2)、然后进行修正后的UVT变换：

其中，L为滑动步长，k为DFT变换的频域索引值，W_M为复旋转因子并且W_M＝e^j2π/M；

(3)、将信号样点x(n)乘以调制序列将频点k的DFT变换平移到k＝0处，根据n-L时刻的DFT变换结果X_n-L(0)计算n时刻的DFT变换输出：

其中，m为调制序列的索引值，每个采样时刻增加1，初始值为0，增加到M-1时，下一采样时刻恢复到初始值0，作为迭代的n-L时刻的DFT变换结果X_n-L(0)，其初始值可以采用传统DFT变换方法得到；

(4)、通过相位修正得到n时刻频点k的DFT变换结果即信号第k个频点的频谱信息：

本发明的目的是这样实现的。

附图说明

图1是连续两帧DFT变换所采用的时域信号样点示意图；

图2是连续两个等效采样时刻16点DFT变换信号样点选取示意图；

图3是SSDFT单个频点信号处理链路结构图；

图4是M＝16、k＝3时SSDFT变换单个频率单元的零极点图；

图5是式(13)对应的mSSDFT变换结构图；

图6是式(14)对应的mSSDFT变换结构图；

图7是带相位补偿的mSSDFT变换信号处理链路结构图；

图8是单频点HDFT变换的结构图；

图9是单个频点mHDFT变换的结构图；

图10是带相位补偿的mHDFT变换信号处理链路结构图；

图11是L＝4、M＝16时单频点mHDFT变换信号处理链路结构图；

图12是各种递归DFT变换的计算误差曲线图。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

为更好地理解本发明，在具体实施方式中，首先对基于逐样点(Sample-by-Sample)运算机制的单步长SDFT(Single-step SDFT，SSDFT)变换和调制单步长SDFT(Modulated Single-step SDFT，mSSDFT)变换的基本原理和实现方法进行分析。在此基础上，将SSDFT推广为更为一般的任意步长SDFT(Hopping DFT，HDFT)算法，分析了HDFT稳定性问题。针对HDFT变换在实际工程应用中存在潜在不稳定和累计误差的问题，研究并提出了一种高精度、恒稳定的调制任意步长SDFT(Modulated Hopping DFT，mHDFT)变换。该变换利用DFT循环频域移位特性有效去除了HDFT谐振器反馈回路中的旋转因子，将HDFT极点精确固定在单位圆上，从而在确保HDFT恒稳定的同时有效降低了累计误差。最后通过软件仿真对几种SDFT变换性能进行了评估和比较，证明了mHDFT变换的优越性。

1、单步长滑动DFT(SSDFT)变换

1.1、单步长滑动DFT变换的引入

对于一个离散时间信号序列x(n)，如果在每个信号样点时刻对其连续进行M点DFT变换，连续两帧DFT计算所使用的信号样点如图1所示。

从图1中可以看出，对于在连续两个信号样点输入时刻进行的两帧DFT变换，第N帧和第N+1帧DFT变换所使用的两个信号时域样点数据块具有很大的相似性，即后一个时刻的样本数据块只是将前一个时刻的样本块中的第一个输入信号样点舍弃，而在最后添加一个新的输入信号样点，而其余M-1个信号样点都是相同的。因此，我们从这个角度来寻找一种DFT快速变换方法。对于两个连续时刻的频谱，已知前一时刻DFT变换结果，如果可以通过简单的递归迭代运算，得到后一时刻DFT变换输出，这在实时信号处理中具有重要的实际意义。

1.2、单步长滑动DFT变换的基本原理

对于一个离散时间信号序列x(n)，其M点DFT变换定义为：

其中，M是DFT变换点数；变量q是一个虚拟变量，q＝n-M+1；W_M＝e^j2π/M是复旋转因子；k是DFT变换的频域索引值；n信号样点的时域索引值。X_n(k)表示的是给定时刻n时DFT变换的第k个频率点值。

通过分析，我们可以看到对于给定时刻n，用于计算M点DFT的信号样点序列为：

x_n＝[x(n-M+1),x(n-M+2),…,x(n-1),x(n)]

(2)

＝[x(q),x(q+1),…,x(q+M-2),x(q+M-1)]

对于下一时刻n+1，用于计算M点DFT的信号样点序列为：

x_n+1＝[x(n-M+2),x(n-M+3),…,x(n),x(n+1)]

(3)

＝[x(q+1),x(q+2),…,x(q+M-1),x(q+M)]

当M＝16时信号样点序列x_n和x_n+1的取值如图2所示。

图2中用实点表示用于DFT运算的信号时域样点。其中，(a)中用于计算n时刻DFT变换的信号样点序列为x_n＝[x(0),x(1),…x(15)]；(b)中用于计算下一时刻n+1DFT变换的信号样点序列为x_n＝[x(1),x(1),…x(16)]。可见在时刻n+1，用于运算的信号时域样点只是增加了当前样点x(16)，去掉了历史样点x(0)。因此，对于连续的两个时刻，用于DFT变换的信号时域样点序列具有很大的相似性，后一个时刻的样点序列只是将前一个时刻的样点序列的第一个输入样点舍弃，而在最后添加一个新的信号样点。

根据式(1)，对于时刻n+1，滑动窗中的M点DFT的第k个频率点，其频谱值为：

对式(4)进行如下变换：

将式(1)代入式(5)得：

其中，X_n(k)是n时刻频点k的DFT计算结果，X_n+1(k)是n+1时刻频点k的DFT计算结果。根据式(6)，n+1时刻的频谱即DFT计算结果X_n+1(k)只需要在前一时刻的频谱X_n(k)基础上加上当前的输入值x(n+1)，减去M个时刻之前的输入值x(n-M+1)，再进行相移即可计算得到。参考式(6)，通过移动索引值n，我们也可以得到单步长滑动DFT(Single-step SlidingDFT，SSDFT)算法X_n(k)与X_n-1(k)的关系表达式：

根据式(7)，SSDFT单个频率单元的实现结构如图3所示。

如图3所示，我们可以看出对于给定单个频点，SSDFT仅仅需要2次复数加法和1次复数乘法就可以完成该频点DFT变换。

同时，我们注意到SSDFT信号处理链路的后半部分由一个谐振器构成，在该谐振器的反馈回路上存在一个旋转因子该旋转因子将在Z平面单位圆上引入一个极点。对式(3)进行Z变换，可以得到SSDFT的传递函数为：

1.3、单步长滑动DFT变换稳定性分析

由式(8)SSDFT的传递函数H_SSDFT(z)可以看出，SSDFT变换具有一个位于Z平面单位圆上极点和M个等间距分布在单位圆上的零点，并且该单极点与一个零点相抵消。虽然SDFT变换不满足因果线性时不变系统的有界输入有界输出(Bounded-Input Bounded-Output，BIBO)稳定条件，即当且仅当系统传递函数的所有极点都位于单位圆内时，因果线性时不变系统是BIBO稳定的。但是，由于H_SDFT(z)中的这个极点被一个零点抵消，从理论分析上来说SSDFT算法是稳定的。图4给出了SSDFT算法单个频率单元的零极点图，其中k＝3，M＝16。

然而，在实际工程应用中必须用有限字长来表示谐振器反馈回路上的旋转因子的量化误差将造成SSDFT的极点位于Z平面单位圆内或位于单位圆外。当极点位于单位圆外时，造成该极点无法被系统零点精确抵消，使SSDFT不再BIBO稳定；当极点位于单位圆内时，虽然此时SSDFT是BIBO稳定的，但是在SDFT每次迭代运算中都会引入的量化误差，在运算结果中引入随迭代次数增加的累积误差。因此，在实际工程应用中SSDFT变换存在潜在的不稳定性。

为了保证SDFT变换在实际工程应用中BIBO稳定，可以在其迭代方程中引入一个阻尼因子r。通过选择合理的r，迫使极点位于单位圆内。修正后的SSDFT(rSSDFT)算法可表示为：

其传输函数为：

显然，rSSDFT变换是通过牺牲系统运算精度来保证系统的BIBO稳定。

2、调制单步长DFT(mSSDFT)变换

2.1、调制单步长DFT变换的基本原理

在第1节中可以看到在实际工程应用中，SSDFT变换中旋转因子的量化误差将造成变换存在潜在的BIBO非稳定性和累积误差的引入，造成DFT变换性能的严重恶化。为了解决这一问题，我们注意到当k＝0时，SSDFT的迭代方程退化为：

X_n(0)＝X_n-1(0)+x(n)-x(n-M) (11)

在式(11)中，X_n(0)的迭代运算中将不会出现旋转因子此时，SDFT变换不仅将无条件BIBO稳定，并且不会引入累积误差。可以利用DFT变换圆周频域移位特性，将频点k的DFT变换平移到k＝0处进行SDFT变换。DFT变换圆周频域移位特性如下：

可见，X(k)可以通过将输入信号样点x(n)乘以调制序列时移到索引k＝0。

根据式(7)，调制单步长DFT(mSSDFT)变换迭代表达式可以写为：

由式(13)，mSSDFT变换结构如图5所示。mSSDFT变换中谐振器的反馈回路无复数旋转因子使mSSDFT变换的极点将精确位于的Z平面单位圆z＝1上，并且不存在有限精度数值误差引入的累积误差。同时，我们注意到当多个DFT频点需要计算时，比如在频域应用时域窗函数，每一个频点需要一个长度为M的延迟缓冲器。

利用的周期性，式(13)中差分方程可以被重新写为：

mSSDFT变换结构可进一步简化，如图6所示。

2.2、mSSDFT变换的调制序列及相位修正

从式(14)中可以看出，调制序列的相位是随索引值m变化而变化的。对于当前采样时刻m，其对应的相位为而在下一个采样时刻m+1，其对应的相位就变为可见，调制序列的相位在m＝0时等于0，以后每增加一个样点，相位就增加

根据式(1)DFT变换定义中虚拟变量q＝n-M+1，当时域索引值n＝M时SSDFT运算是从q＝1开始，即此时被分析信号的相位对应于时刻q＝1。同时，我们注意到在整个迭代运算中，每步进一个样点进行一次迭代运算，调制序列的相位就增加当q＝1时，被分析信号的相位是常数而不是0。因此，调制信号序列和被分析的信号存在的相位差，X_n(k)与X_n(0)的相位关系为：

采用式(16)中的旋转因子对X_n(0)进行相位校正即可得到信号第k个频点的频谱信息。

2.3、单频点mSSDFT变换的信号处理链路结构

由图3和图6比较可知，调制序列是随m的变化而变化，必须在每个信号样点输入时重新计算为了解决这一问题，这里采用一个复数振荡器来实现图6中的调制序列将式(16)重写为：

调制序列是以M为周期的，每M个样点就自动从(对应于m＝0)开始，避免了计算过程中累积误差的产生。采用递归方法的mSSDFT信号处理链路结构如图7所示。

由于|X_n(k)|＝|X_n(0)|，如果只需要k频点DFT变换幅度输出，则图7中最后一级用于相位修正的复数乘法可以省略。

3、任意步长滑动DFT(HDFT)变换

3.1、任意步长滑动DFT变换的引入

从前面SDFT变换论述可以看出，SDFT和mSDFT变换采用的是逐样本(sample-by-sample)计算方式，即每输入一个信号采样样点立即进行一次DFT运算，更新输出对应频点上频谱信息。这使DFT变换产生的信号频点值输出速率等于信号输入样点速率，很好的满足了强实时信号处理的要求。

而在某些实时性要求不高的特殊信号处理应用中，并不要求在每信号输入样点时刻都计算产生对应频点DFT输出，只需每间隔L(0<L<M)个时钟周期更新一次对应频点DFT输出。例如：在可变时间分辨率的信号时频分析中，要求调节连续两次DFT输出的时间间隔L(0<L<M)是任意可调节的。通过调节时间间隔L控制每次DFT运算输入信号样点更新个数，即调节DFT变换窗口的滑动步进L，实现任意时间尺度的信号时频分析。虽然单步长滑动DFT变换能够根据滑动步进L在对应时刻更新信号频谱输出，但是单步长滑动DFT变换递归运算机制要求在每个信号输入样点时刻都必须进行DFT变换，而DFT变换结果仅在每L个时钟周期有效输出，其余时刻都被丢弃。因此单步长滑动DFT变换在这种实时性要求不高的信号处理应用中计算效率相当低。

针对这一问题，我们基于SSDFT变换递归迭代计算思想，对其滑动步进进行了一般化推广，推导出了一种任意步长滑动DFT(HDFT)变换。该变换能够根据实际应用需求任意设置DFT窗口的滑动步长L，每间隔L个时钟周期对信号进行一次DFT变换运算，大大提高了算法的计算效率。

3.2、HDFT变换的基本原理

设滑动步长为L(0<L<M)，它表示每次计算DFT输出时，信号输入样点更新个数。显然对于SDFT和mSDFT来说，其滑动步长L＝1。下面我们将滑动步长将单步长SDFT推广到任意步长SDFT。

3.2.1、滑动步长L＝2

首先考虑L＝2情况。由SSDFT变换

其中，d(n)＝x(n)-x(n-M)。

X_n-1(k)和X_n-2(k)之间的关系如下：

将式(19)代入式(18)得到：

这样，当滑动步长L＝2时，HDFT变换的递归方程为：

3.2.2、任意滑动步长L

对于任意滑动步长L，X_n(k)和X_n-L(k)之间的变换关系可以通过在式(8)中用X_n(k)连续代入X_n-1(k)L次得到。为了推导方便，我们这里主要考虑L＝2^a，a≥0的情况，其他情况可采用同样的推导方法很方便的得到。经过L次连续代入，任意步长滑动DFT(HDFT)表达式为：

这里定义为L点更新矢量变换(Updating Vector Transform，UVT)中的第k个频点信息。

将式(23)代入式(22)，HDFT递归计算公式为：

可见，n时刻的DFT输出可根据n-L时刻的DFT输出和的变换结果直接计算得到。

单频点HDFT变换结构图如图8所示。

3.2.3、UVT的高效计算方法

根据式(23)UVT的定义，我们发现UVT的计算公式与DFT的定义非常相似，可以采用FFT算法来实现的高效计算。这里我们采用按时间抽取(DIT)的基-2FFT算法来计算UVT变换。变换可表示为：

其中，

DIT基-2FFT算法按奇数点和偶数点将长度为L的序列d(n)分解两个为L/2的序列，分别对两个子序列进行FFT变换。而两个长度为L/2子序列又可以进一步分解为两个长度为L/4奇偶序列，分别进行FFT变换。以此类推，连续对序列进行奇偶抽取分解，直到将序列成为单点序列为止，即可计算得到3.2.4、HDFT传递函数

将式(24)进行Z变换，可以得到单频点HDFT的传递函数为：

从HDFT算法传递函数中可以看出系统存在一个位于Z平面单位圆上的极点，和M个均匀分布在单位圆上的零点。

3.3、HDFT变换稳定性分析

分析式(28)HDFT变换的传递函数H_HDFT(z)，单频点HDFT变换具有一个位于Z平面单位圆上极点和M个等间距分布在单位圆上的零点，并且该单极点与其中一个零点相抵消。因此，从理论分析上来说HDFT变换是BIBO稳定的。

但在实际工程应用中，HDFT变换存在与SSDFT变换相同的稳定性和累积误差问题。即在使用有限字长来表示谐振器反馈回路上的旋转因子时量化误差的引入而造成HDFT的极点位于Z平面单位圆内或位于单位圆外。当极点位于单位圆外时，造成该极点无法被系统零点精确抵消，使HDFT变换不再BIBO稳定；当极点位于单位圆内时，虽然此时HDFT变换是BIBO稳定的，但是在HDFT每次迭代运算中都会引入的量化误差，在运算结果中引入随迭代次数增加的累积误差。因此，在实际工程应用中，HDFT变换存在潜在的不稳定性和累积误差的缺陷，严重恶化算法性能。

4、调制任意步长DFT(mHDFT)算法

针对HDFT变换中旋转因子的量化误差将造成算法存在潜在的BIBO非稳定性和累积误差引入的问题。本发明借鉴mSDFT变换的思想，提出了一种全新的高精度、恒稳定调制任意步长DFT(Modulated Hopping DFT，mHDFT)变换方法。该方法利用DFT变换圆周频域移位特性，通过直接将第k个DFT频点平移到k＝0处，然后根据n-L时刻的DFT变换结果计算n时刻的DFT变换输出。通过mHDFT算法有效去除了HDFT谐振器反馈回路中的旋转因子，将HDFT变换极点精确固定在单位圆上，从而在确保HDFT变换恒稳定的同时有效提高了算法计算精度。

4.1、调制任意步长DFT变换的基本原理

根据式(24)HDFT变换递归计算公式，对于滑动步长为L，k＝0处频点的HDFT变换输出为：

由式(29)可以看出X_n(0)的迭代运算中将不会出现旋转因子此时，HDFT变换不仅是无条件BIBO稳定的，并且不会引入累积误差。因此，我们可以利用这一特性来消除HDFT存在的潜在BIBO非稳定性和累积误差引入的问题。

利用式(12)给出的DFT圆周频域移位特性，通过将输入信号样点x(n)乘以调制序列将频点k的DFT变换平移到k＝0处，根据n-L时刻的DFT变换结果计算n时刻的DFT变换输出。最后，通过相位修正得到n时刻频点k的DFT变换结果。

4.1.1、mHDFT变换的递归运算

在对信号输入样点进行圆周频域移位处理时，我们必须注意到HDFT中UVT移位处理的特殊性。将UVT变换式(23)重写为：

可以看出在中存在一个与UVT变换中虚拟分量t无关的固定相移因子虽然当k＝0时，使该相移因子没有明显出现在式(29)中，但是在输入信号样点整个圆周频域移位处理过程中，频点k处的UVT所包含的固定相位因子并不会受到影响。因此，圆周频域移位处理后的mHDFT算法可以表示为：

由于式(31)可以简化为：

其中，定义为修正后的UVT(Revised UVT，rUVT)变换：

由式(32)，单个频点mHDFT变换的结构如图9所示。

4.1.2、rUVT变换的高效计算

根据式(33)rUVT的定义，我们发现rUVT的计算公式与DFT的定义相同，可以采用FFT算法来实现的高效计算。这里我们同样采用按时间抽取(DIT)的基-2FFT算法来计算rUVT变换。变换可表示为：

其中，

由式(34)，L点信号序列d(n)被分成两个长度为L/2的子序列，分别对应d(n)中奇数索引和偶数索引的子序列。这种信号序列的抽取可以不断重复，直到导致的序列简化为1点序列。可以根据抽取得到的两个子序列的DFT变换直接合成得到。

4.2、mHDFT变换的调制序列及相位修正

从式中可以看出，调制序列的相位是随索引值m变化而变化的。同时，我们也注意到对于滑动步长为L的mHDFT变换是根据n-L时刻的k频点的DFT变换结果，递归计算得到n时刻的k频点DFT变换输出。因此，mHDFT变换连续两次输出之间间隔为L个采样周期。这将造成在每次迭代运算中，在调制序列中都会引入的固定相移。即：对于当前DFT变换输出时刻m，其对应的相位为而在下一个当前DFT变换输出时刻m+L，其对应的相位就变为可见，调制序列的相位在m＝0时等于0，以后每步进一次，相位就增加

根据式(1)DFT变换定义，我们同样可以得到与SSDFT相识的结论。当时域索引值n＝M时mHDFT运算是从q＝1开始，即此时被分析信号的相位对应于时刻q＝1。同时，我们注意到在整个迭代运算中，每步进一个样点进行一次迭代运算，调制序列的相位就增加当q＝1时，被分析信号的相位是常数而不是0。因此，调制信号序列和被分析的信号存在的相位差，X_n(k)与X_n(0)的相位关系为：

采用式(38)中的旋转因子对X_n(0)进行相位校正即可得到信号第k个频点的频谱信息。

4.3、单频点mHDFT变换的信号处理链路结构

由于调制序列是随m的变化而变化，必须在每个信号样点输入时重新计算为了解决这一问题，同样采用mSSDFT中的复数振荡器来取代图9中的调制序列将式(37)重写为：

因此，采用递归方法的mHDFT变换信号处理链路结构图如图10所示。

由式(32)和图10可以看出，当滑动步长L＝1时，mHDFT变换和mSSDT变换的递归方程以及信号处理结构图相同，mHDFT变换直接退化为mSSDFT变换。可见，mSSDFT变换是mHDFT变换的特殊情况。由于|X_n(k)|＝|X_n(0)|，如果只需要k频点DFT变换幅度输出，则图10中最后一级用于相位修正的复数乘法可以省略。图11给出了当L＝4和M＝16时单频点mHDFT变换信号处理链路结构。

5、mHDFT变换的计算量

M点mHDFT变换的计算量可以通过以下公式计算出来：

C_A＝L+M(1+log₂L) (41)

其中，C_M、C_A分别表示mHDFT变换中复数乘运算和复数加运算的数量。

如果一个复信号的M点DFT变换所有频点信息都需要计算产生，那么各种变换所需的计算量如表1所示。

表1

可以看到FFT算法的计算量由频点个数M直接决定，而HDFT和mHDFT的计算量不仅于M有关，也与滑动步进L有关。mHDFT的计算量略高于SDFT和HDFT，这是由于mHDFT的相位调制造成的。但是，当M或者L较小时，mHDFT相比FFT计算量优势明显。接下来，我们对各种变换计算单个DFT频点所需计算量进行了理论分析，结果如表2所示。

表2

mHDFT变换所需复数加运算个数与HDFT变换相同，而在计算下一DFT频点时所需复数乘法运算比HDFT多了一个。多出的这个乘法器是由相位调制所引入的。同时C_M会随DFT频点个数的增加而成比例增加。因此，对于只计算部分选定DFT频点时，mHDFT运算量明显优于DFT、SDFT和mSDFT。下一节我们将证明计算量的微小增加对于提高运算精度是值得的。

6、仿真结果

各种前面讨论的递归DFT变换的计算误差可由下式计算得到

其中表示对于长度为M的数据段标准DFT频点计算输出，表示SDFT、rSDFT、HDFT和mHDFT算法计算输出。零均值、标准方差为1的高斯白噪声信号用来对算法精度进行评估。其中DFT点数M，rSDFT中阻尼因子r和mHDFt中的滑动步长分别设定为16、0.9999999和4。在64位双精度MATLAB环境中采用各种算法对K＝3对应DFT频点信息进行了计算，其误差结果如图12所示。

mHDFT变换的计算误差是所有算法中最小的，在10^-15附近波动。mHDFT的计算精度相比原始的HDFT算法精度提高了10倍。因此mHDFT在保证算法稳定的同时提供了最高的计算精度。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims

1.一种离散数字信号任意步长的滑动离散傅里叶变换方法，其特征在于，包括以下步骤：

d(n)＝x(n)-x(n-M)，

其中，M是DFT变换点数，n是信号样点的时域索引；

(2)、然后对差值信号d(n)进行修正后的UVT变换：

G_{n}^{L} (k) = Σ_{t = 0}^{L - 1} d (n - t) W_{M}^{- (t + L - 1) k},

其中，是L点修正后的更新矢量变换即UVT变换中的第k个频点信息，L为滑动步长，k为DFT变换的频域索引值，W_M为复旋转因子并且W_M＝e^j2π/M；

X_{n} (0) = X_{n - L} (0) + W_{M}^{k m} G_{n}^{L} (k);

其中，所述调制序列采用一个复数振荡器来实现，其形式为：

W_{M}^{k m} = W_{M}^{k (m - 1)} W_{M}^{k};

X_{n} (k) = W_{M}^{- k (m + L)} X_{n} (0) .