1.2.1 数字信号处理部分知识重温(1)从DFT(离散傅里叶变换)到FFT(快速傅里叶变换)(2)

写在前面:本文不作为一篇严谨的理论推导文章,仅为博主学习时的随手小记,个人的主观认知为主,以求日后能在较快时间内理清思路,找到当下的学习状态。文章中会尽量避免出现过多公式及推导,如果能有幸给其他朋友带来一些帮助,不甚荣幸。

感谢以下文章及老师的指导:

【数字信号处理(北京交通大学 陈后金)】 2.1 离散傅里叶变换(DFT) - 4 DFT性质2_哔哩哔哩_bilibili

这篇文章我写的不是很满意,感觉自己还是没有完全讲明白,日后会再次修改,建议直接观看网课,参考我在文章中的注解即可

  • 一、从DFT到FFT的导入

承接上文,我们引入DFT,以此来计算离散周期序列的频谱,估算其他类型信号的频谱,计算线性卷积和循环卷积等。DFT十分有用,但就是有一个缺点,体现在计算量(见下图)上,我们以一道题为例,找出DFT运算量规律:

因此DFT总运算量=N^2+N(N-1),包含一项指数函数,当N达到一定规模时,这是一个惊人的数据!

因此我们希望设计一套新的算法提速,以下是可能的路径:

  • 二、基于时间抽取的基2-FFT算法   

将长序列按照下标的奇数项和偶数项(请注意,这里讨论奇偶,是存储序列的数组的下标,而不是序列本身数值的奇偶),分别归入两个子序列进行DFT运算(运用此方法有一个前提:原序列的项数必须是2的指数,如果不是需要补零来凑),降低DFT处理序列的长度,降低运算量。

此方法可递归使用,将子序列继续按奇偶项划分,进一步提速,直到产出的子序列仅剩2个元素

这里含有一个重点:原序列二分后进行DFT运算,每往下二分一次,得到序列长度更小的子序列时,原序列的元素就需要重新按照元素在原序列下标的奇数项和偶数项分别归入两个子序列,我们称呼为“奇偶重排”。

这个加速DFT运算的方法就是FFT(快速傅里叶变换),具体来说,叫做基于时间抽取的基2-FFT算法。

    下图是DFT算法,根据将时域序列二分的思想改写后的表达式。

注:在第三行第二项中,将Wn^(2k+1)m分成了两部分:Wn^m,(Wn^2km=Wn/2^km),因此在形式上比第三行第一项还要再多乘一个复数Wn^m,不要遗漏了

我们将长序列分成短序列运算DFT,必须基于一个前提:两个子序列的DFT运算结果,必须可以通过某种方式,运算得到原序列的DFT运算结果。因此,具体实现步骤如上图。

因为子序列的长度必然只有原序列的1/2,所以两个子序列运算得出的结果,也注定只能表示原序列的其中一半的情况(第一行反应了原序列的前半段频域序列),因此我们令m=m+N/2,表示原序列后半段的频域序列。

注:在第三行中,由于子序列的长度由原序列的N变为N/2,所以N/2便是子序列(内含复数Wn/2^km)X1[m+N/2],X2[m+N/2]的序列长度(周期),由复数的周期性,Wn/2^k(m+n/2)=Wn/2^km,由此X1[m+N/2]=X1[m](X2[m+N/2]同理)

从这我们还可以发现一件事:原序列不断的按二分的方式拆解成更多的短序列,除了初始级序列长度只有2的子序列,其它级(长度为4、8...的)的子序列进行DFT运算,实际上并不需要把这个序列直接进行DFT运算,他们的DFT结果(频域序列X[m])实际上是由更低一级的DFT结果(频域序列)通过复数加法和复数乘法(也就是上图的运算方法)组合得来的。

真正进行了DFT运算,将时域序列通过计算得出频域序列(DFT结果)的,完成了时域到频域的转化的,实际上只有初始级(长度只有2的子序列)。

由此我们导出了子序列DFT(除了初始级)得到的频域序列合成原序列的频域序列的计算表达式、矩阵表达式(太麻烦了看看就行)和运算流图(蝶形图)。(见上图)

我将计算过程中的规律总结如下:

    1.表达式中旋转因子中的N不由子序列长度决定,由将被合成的序列长度决定

    2.(可以观察蝶形图)两个子序列合成原序列的元素,第二个子序列的元素(自行决定哪个是第一个,哪个是第二个)在计算时要配上旋转因子,旋转因子的m取决于被合成的元素位于原序列前(后)半部分的第几个位置。

    3.(可以观察蝶形图)合成原序列的后半部分的元素时,第二个子序列的元素“*(+1)*旋转因子”;合成原序列的后半部分的元素时,第二个子序列的元素“*(—1)*旋转因子”。

    4.每一级频域序列结果的得出,只和更低一级的两个长度为其一半的两个子序列有关。那么在计算机实现时,我们短序列DFT合成长序列DFT的结果,可以直接覆盖在两个子序列的存储空间(数组)之上,而不需要建立大量空间,存储每一级的中间运算结果。

    5.还有一些规律,总结在例1之后,现在暂时先不用跳下去看

我们还得到了在初始级序列中真正出现的DFT运算表达式。(见上图)

需要说明的是:

  1. Wn^0=1(n可取任何正整数)
  1. 第二行的表达式原为X[1]=x[0]+W2^1*x[1],根据旋转因子的对称性,Wn^(mk+n/2)=—Wn^mk。所以X[1]=x[0]+W2^(0+2/2)*x[1]=x[0]—W2^0*x[1]。

3. 在初始级进行的是DFT运算而不是序列的二分,不需要进行奇偶重排。

在此再次列出DFT的运算表达式,以作参照。

 我们将上面两种适用于不同级的序列DFT再次列出,归纳如下:第一行适用于除初始级序列,第二行适用于初始级序列(长度=2)

我们注意到,无论是哪一级,当它的运算方式用矩阵表达时,它们的系数矩阵是一致的,这也就为后面将这两种在不同级的运算方法统一起来奠定基础。

例题1:用上述方法,完成八点时域序列的DFT

    我们采用自顶向下的解题方法,从原序列二分,往下得到子序列。值得注意的是,在二分的过程中,原序列开始进行了一次元素的奇偶重排(见下图最左边的子序列排布)

    关于四点DFT得到的两子序列的频域序列,组合得到原序列DFT后得到的频域序列的过程,可以参考上文总结的计算规律,此处不展开。

    完成运算后,我们将四点DFT展开,也就是对原来的两个子序列再次二分,得到四个子序列,这个时候,两个子序列的元素在分化为四个子序列时再次奇偶重排(见下图最左边的子序列排布)

最后,我们将初始级的2点DFT展开,得到的完整蝶形图如下图。

需要说明的是,为了保持计算过程的一致性和对称性,我们利用了旋转因子的可约性Wn^m=Wc*n^c*m,将初始级的W2^0变为W8^0,将倒数第二级的W4^0,W4^1分别变为W8^0,W8^2。

从此题中,我们还可以发现两个计算规律(主要体现在蝶形图中):

    1.初始级的运算,是间隔2^0=1个元素交叉相乘;第二级的运算,是间隔2^1=2个元素交叉相乘....以此类推。

    2.每一级DFT运算出现的旋转因子,都会间隔的出现在下一级的DFT运算中,中间新加入的旋转因子上标=(相邻旋转因子上标之差)/2。

    3.(重要!!)我们通过若干次二分后,原序列的元素经过多次奇偶重排得到了新的排布方式(上图第三列),用二进制表示下标的话(原序列长度是2的几次方,这里用二进制表达就统一用几位二进制)见上图第二列。而如果我们将对应的下标的二进制表达式颠倒顺序(上图第一列),正好符合顺序排列。

这就启示我们:只要将长度为N=2^c的原序列元素顺序下标按c位二进制表示,再颠倒顺序,就能得到经过若干次二分后,最终各元素的排布情况。我们称这种思想为:倒序运算。

例题2:用上述方法,完成四点时域序列的DFT

  本题不做讲解,可以用上一题的思路进行分析。

  • 三、频率抽取基2-FFT算法

    这个算法与时域抽取的基2-FFT算法相似,都是将原序列分为两部分以此降低运算量。时域的方法是,将序列各元素按照奇偶存储顺序(下标)进行二分,最终可以分别合成出频域序列的前半部分和后半部分。在此可以提前剧透一下,频率抽取的方法是将序列按照前半部分和后半部分进行二分,最终分别合成出频域序列的奇数下标部分和偶数下标部分。

所以我们可以做出总结:奇偶拆分得到频率上的前后,前后拆分得到频率上的奇偶。

    另:频率抽取基2-FFT算法不需要对时域序列奇偶重排

上图第二步,第二部分的旋转因子Wn^m(k+N/2)=Wn^mk*Wn^m(N/2),当m为奇数时,Wn^m(N/2)=—1;当m为偶数时,Wn^m(N/2)=1。故Wn^m(N/2)=(-1)^m。

于是,原序列的频域序列元素可按照下标的奇偶性,分别由时域序列前半部分和后半部分的元素运算得到:(见下图)

注:为了表示奇偶,频域序列元素X[m]偶数下标用2m表示,奇数下标用2m+1表示,同时表达式内的m也要一并换成2m或2m+1(见下图)

为了运算的方便,我们将第一行和第二行方框内的式子设为x1[k],x2[k](见下图)

由此,我们可以得出时域的前半部分和后半部分的序列得到x1[k]和x2[k]的运算公式,运算矩阵和蝶形图。(见下图)

  经过时域的两个序列运算得出的结果x1[k]和x2[k]依旧是时域序列!!而非已经能得到频域序列,x1[k]和x2[k]合成频域序列的方法见下图。

 在此我们可以窥见时域抽取FFT运算和频域抽取FFT运算的一些区别:

  1. 时域抽取FFT运算,自初始级往后,运算输入和输出的结果都是频域上的数值而非时域上的数值,并且运算的输入是两个短序列,得出的结果是一个长序列的前半部分和后半部分。也就是说,我们运算次数越多,序列会越来越长

2.频域抽取FFT运算,运算输入和输出的结果都是时域上的数值而非频域上的数值,并且运算的输入是一个长序列的前半部分和后半部分,得出的结果是两个短序列。也就是说,我们运算次数越多,序列会越来越短;且在最后一级运算时,序列被划分到长度为2的序列,这时的情况我们放在下文讨论。

和时域抽取基2-FFT算法一样,频域抽取基2-FFT算法也是到了长度为2的序列时才真正的实现时域到频域的转换,我们再次给出DFT公式以作参考(下图),得到运算矩阵如下图,可以发现,与时域抽取基2-FFT算法初始级的算法完全相同,在此便不再分析。

 由此可说明时域抽取FFT和频域抽取FFT的另一个区别:

 时域抽取FFT将时域序列转换为频域序列放在初始级,也就是一开始就将时域短序列转换到了频域,然后频域序列不断运算,长度开始“滚雪球”,最后得出了原时域序列对应的频域序列

频域抽取FFT将时域序列转换为频域序列放在最后一级。首先它是“逆着滚雪球”,将时域序列越分越小,直到序列长度=2时,也就是直到最后才将时域短序列转换到了频域,得出了原时域序列对应的频域序列。

我们以8点序列的频域抽取基2-FFT算法为例,来演示一下完整的信号流图(蝶形图)

借助这幅图我们还发现:时域抽取FFT运算时,需要先对时域序列“倒序运算”以完成奇偶重排,再输入进行运算,而频率抽取FFT运算,直接将时域序列顺序输入即可,但输出的结果需要完成“倒序运算”的逆过程才能得到频域序列的顺序排列

  • 四、其他基时间抽取FFT算法(大概看看就行)

  ①基4

    序列长度必须满足4的幂次,否则原序列需要补零。

我们在此处再次列出时域的基2-FFT算法中,DFT算法根据将时域序列二分的思想改写后的表达式:

以此为参照,我们可列出时域序列四分后,X[m]的四部分的表达式:

根据旋转因子的,进一步化简可得到:

  ②基3