滑动DFT算法研究

第一章绪论1滑动DFT算法研究黄寒华第一章绪论1.1引言离散傅立叶变换（DFT）是数字信号处理领域中遇到的最常用，功能最强大的处理方法之一，是数字信号处理的核心。利用DFT变换，可以实现信号与系统的分析与综合。同时，许多算法，如相关、滤波、谱估计等也可以转化为DFT来实现，因此DFT在许多工程领域都已获得了许多应用。离散傅立叶变换的目的是把信号由时域变换到频域，从而可以在频域分析处理信息[1]。同时得到的结果也可以再由离散傅立叶逆变换到时域，其变换对的算法表达式为：DFT：1010,)()(NnnkNNkWnxkX(1.1)IDFT：1010,)(1)(NknkNNnWkXNnx(1.2)其中，NjNeW2(1.3)虽然离散傅立叶变换有着明确的物理意义，是确定时域序列的频率成分的最直接的数学过程，也便于计算机进行处理，但它的运算效率非常低，对实时问题处理意义不大，因此长期以来并没有得到真正广泛的应用。直到库利（Cooley）和图基（Turkey）提出了一个计算DFT的高效算法--快速傅立叶变换（FFT），才使离散傅立叶变换由理论走向了实践，成为了数字信号分析的强有力工具，被广泛应用于各种数字信号处理系统中。继Cooley和Turkey之后，又有许多人致力于进一步减少DFT的运算量，先后提出了一些改进算法，如winograd算法[2]等。随着DFT的发展，DFT的应用领域还在扩大，因此研究DFT的快速算法有着极其重要的意义。1.2快速傅立叶变换算法介绍从离散傅立叶变换公式可以看出，直接计算DFT需要N2次复数乘法和N（N-1）复数加法，而1次复数乘法要做4次实数乘法和2次实数加法，1次复数加法要做2次实数加法。直接计算DFT的计算量与变换区间长度N的平方成正比，当N较大时，计算量剧烈增加。因此，虽然离散傅立叶变换第一章绪论2物理意义巨大，但由于计算机处理速率的限制，直接进行离散傅立叶变换不太可能。1965年，Cooley和Turkey发表了一篇关于计算DFT的高效算法的论文，FFT算法使得上千点的DFT运算由原来的只能在研究机构和大学的计算中心完成而变成了在家庭计算机上只用几秒钟就可以完成。FFT算法将离散傅立叶变换的运算量由N2降为Nlog2N,对于1024个采样，运算复杂度降低了两个数量级。使离散傅立叶变换在工程技术领域中真正得到了广泛的应用。FFT算法通常采用的类型有两种[3]。一种是针对N等于2的整数次幂的算法，如按时间抽取的基2、基4算法、按频率抽取的基2、基4算法和分裂基算法等；另一种是N不等于2的整数次幂的算法，如任意基N算法、winograd算法等。以时间抽取的基2算法为例，设序列)(nx的长度N为2的整数次幂，即lN2。将式1.1表示成两个部分之和的形式，一部分是n为偶数时的)(nx，另一部分是n为奇数时的)(nx，则10)()(NnnkNWnxkX120)12(1202)12()2(NnknNNnnkNWnxWnx12021202)12()2(NnnkNkNNnnkNWnxWWnx(1.3)因为N为2的整数次幂，为偶数，则nkNnkNjnkNjnkNWeeW222222(1.4)则式1.3可表示为12021202)12()2()(NnnkNkNNnnkNWnxWWnxkX(1.5)因此，式1.5中)(kX分成了不同的两个N/2点DFT。一个N点DFT分解成两个N/2点DFT，其中一个乘以因子kNW后，两者再相加。我们知道，N点DFT需要N2次复数乘法和N（N-1）复数加法，则计算式1.5中所需的复乘次数为：NNNNN2/2/2/222(1.6)计算式1.5中所需的复加次数为：2/2/2/222NNNN(1.7)对于2N，采用式1.5计算DFT的运算复杂度比直接计算要降低许多。特别要指出的是，若N为2的整数次幂，则可以对结果中的两个N/2点DFT再做如式1.5的处理，直至DFT的长度为2为止。第一章绪论3这样的递推算法，使得快速傅立叶变换复乘和复加次数变成了NN2log，与直接计算比较，运算次数降低了NN2log/倍。按时间抽取基4算法、按频率抽取的基2和基4算法与按时间抽取基2算法原理完全类似，基4算法可以不用乘法而只用加法来实现，其比基2算法更为有效。分裂基算法的原理同时采用了基2和基4算法，其基本思路是对偶序号输出使用基2算法，对奇序号输出使用基4算法，由于分裂基算法在目前已知的所有针对MN2的算法中具有最小的乘法次数和加法次数，并且具有和Cooley-Tukey算法同样好的结构，被认为是最好的快速傅里叶变换算法。研究表明，该算法最接近理论上所需乘法次数的最小值[4]。任意基N算法，也称素因子算法，因为其N不是2的整数次幂，无法基2或基4分解，而是将其N因式分解llNNNNNNN21321(1.7)将输入序列分成N1个长度为N2l的序列，再将N2l分成N2个长度为N3l的序列，以此类推，从而达到降低计算量的目的。winograd算法[5]是以数据寻址理论进行卷积运算为基础的一种算法，其复乘次数比FFT少，但控制复杂，对要求算法简单，模块性好的数字信号处理器来说，winograd算法并不一定比FFT快。快速傅立叶变换模块性好，在硬件实现上方便，目前普遍使用的仍然是FFT算法。上面提到的这些快速算法其着重点都是减少复乘、复加次数，降低计算的复杂度，数据的传输是全局性的。以往的算法设计思想主要强调的是尽量少的乘法次数，但随着现代超大规模技术的发展和并行计算技术的使用，不仅要求降低算法的计算复杂度，而且要提高算法的并行度和使之模块化，减少模块间的耦合度，同时还要满足实时处理的要求。从这方面来讲，经典的FFT算法就有一定的缺点。相比较而言，某些用其他方式实现的DFT算法(如本文推导的滑动DFT算法)具有数据传递简单，通信局部化，实时性强的优点。第一章绪论4第二章滑动DFT算法原理本章介绍了滑动DFT算法的算法原理。2.1滑动DFT算法概述在某些信号处理的应用领域，特别是在一些实时应用中，希望在频域上对离散信号进行连续分析。分析不同时刻其输入样本对应的频谱，这就相当于用一个固定长度随时间滑动的滑动窗口来选择样本，这种在一个滑动窗口内计算N点DFT的算法称为滑动DFT算法。逐点滑动DFT原理图如图2.1所示。1Z1Z1ZDFT)(nx)1(nx)2(nx)2(Nnx输入数据频率分量0k1k2k2Nk1Nk)(kX图2.1逐点滑动DFT原理图)1(Nnx注：这只是某时刻长度为N的窗里面的N点的DFT。按照传统的FFT算法，对于某一时刻，用FFT计算出其所有频谱，如果要计算下一时刻的频谱，则再进行一次FFT运算，这两次FFT运算是孤立进行的，它们之间没有任何关系。而实际上，对于连续的两个时刻，我们会发现，其窗口中的样本有着很大的相似性，后一个时刻的样本只是将前一个时刻的样本的第一个输入舍弃，而在最后添加一个新的样本。如果两个时刻相距不远，则后一个时刻的样本只是将前一个时刻的样本的前几个输入舍弃，而在最后添加几个新的样本。不同时刻的窗口中的样本只有一个或几个不同，其时域中的相似性必然会使其频谱有着一定的联系。滑动DFT算法正是基于这样的思想而进行的探讨，对于两个连续时刻的频谱，已知前一时刻频谱，则可以通过简单的递推运算，得到后一时刻的频谱，这在连续的实时谱分析中有着重要的实际意义。第一章绪论52.2滑动DFT算法推导2.2.1滑动DFT算法推导[6]对于一个离散时间信号，其DFT变换为：10/2)()(NnNnkjenxkX（2.1）给定离散时间信号)(nx，则经过DFT变换后可以得到此信号对应的N点频谱)(kX。设时刻q，滑动窗选中N个样本的序列为：x)(),1(,),2(),1()(qxqxNqxNqxq（2.2）根据傅立叶变换理论，时刻q，在滑动窗中的N点DFT的第k个频率单元，其频谱值为：10/2)()(NnNnkjqenxkX（2.3）通常，DFT结果的索引是频域索引k，但是在式2.3中，DFT结果的索引既有频域索引k，也有时域索引q，因而)(kXq所表示的为q时刻第k个频率单元的值。例如，图2.2为20点滑动DFT的滑动窗，图2.2（a）为时刻为0时的滑动窗口中的20采样点，即为计算)(0kX的窗口，图2.2（b）则为时刻为1时的滑动窗口中的20采样点，即为计算)(1kX的窗口，以此类推。图2.2（a）时刻为0的20点滑动窗第一章绪论6为了推导的简化与清晰，将时间索引修改(q取代1Nq推导结束后再还原），N个样本序列可表示为：x)1(),2(,),1(),()(NqxNqxqxqxq（2.4）根据傅立叶变换理论，时刻q在滑动窗中的N点DFT的第k个频率单元，其频谱值为：10/2)()(NnNnkjqeqnxkX（2.5）根据式2.5，时刻q+1，在滑动窗中的N点DFT的第k个频率单元，其频谱值为：10/21)1()(NnNnkjqeqnxkX（2.6）式（2.6）中，令p=n+1，则NpNkpjqeqpxkX1/)1(21)()(（2.7）改变求和号的上、下限，得NkNjNkjNpNkpjqeNqxeqxeqpxkX/)1(2/)1(210/)1(21)()()()(（2.8）则NNkjNpNpkjNkjqeNqxqxeqpxekX/210/2/21)()()()(（2.9）即可得到滑动N点DFT的递推表达式：图2.2（b）时刻为1的20点滑动窗图2.2滑动窗示意图第一章绪论7)()()()(/21qxNqxkXekXqNkjq（2.10）也即)1()1()()(1/2qxNqxkXekXqNkjq（2.11）其中，)(kXq是q时刻单个频率单元的结果，)(1kXq是1q时刻单个频率单元的结果。还原式2.11中的时间索引，将时域中的q还原为1Nq，并且输入采样点和输出采样点使用同样的索引n，则得到滑动DFT的时域差分方程：)()()()(1/2NnxnxkXekXnNkjn（2.12）从式2.12可以看出，计算)(kXn，只要通过前一个点)(1kXn减去)(Nnx，加上现在的)(nx，再进行相移来计算。因此，滑动DFT仅仅需要2次实加和一次复乘（4次实乘）就可以完成每个输出点的计算。第三章滑动DFT算法的改进3.1多点滑动DFT3.1.1多点滑动DFT算法推导前面讨论的滑动DFT算法是逐点滑动的DFT。处理的是连续时刻采样的情况。但如果为非连续采样，即多点滑动呢？当然我们可以逐点滑动，连续几次实现，但如果相隔的采样点数大于N2log（其中N为滑动窗中的总采样点数），当计算全部频谱时，滑动DFT还不如FFT计算简单[12]。因此，接下来我们讨论多点滑动的DFT实现。设时域信号波形的采样离散值为()xn，傅里叶变换为()iXk，现求以ip为起始位置的N点傅里叶变换()ipXk，其中p为两窗之间的间隔，2Lp，根据DFT变换公式，由于n只对旋转因子knNW起作用，而对输入()xn只是下标的延迟，故以ip为起点的DFT变换可写成下式：pNpnkpnNjpienxkX1)(2)()(第一章绪论8pNNnkpnNjNpnkpnNjenxenx1)(21)(2)()(10)(21)(210)(2)()()(pnkpnNjpNNnkpnNjNnkpnNjenxenxenxpkNjpmmkNjpmkpNmNjNnpkNjnkNjeemxeNmxee