第六讲离散Fourier变换续与快速Fourier变换

第六讲离散Fourier变换续与快速Fourier变换目录一、离散Fourier变换（DFT）...............................................................................................1二、线性卷积与圆周卷积.......................................................................................................22.1线性卷积....................................................................................................................22.2圆周卷积....................................................................................................................32.3时域卷积定理............................................................................................................3三、离散Fourier变换的应用.................................................................................................5四、考察离散Fourier变换的计算量.....................................................................................5五、快速Fourier变换.............................................................................................................5六、FBP的频域滤波实现问题..............................................................................................11一、离散Fourier变换（DFT）对序列xk，其离散Fourier变换定义为(1)120,0,11,nkNiNkXnxkNen离散Fourier逆变换由下式给出(2)1201,0,11,nkNiNnxkXnekNN令(3)12iNWe(1)式写成矩阵形式(4)111110000012112101001111NNNNNXx注：有教材将离散Fouerier变换定义在周期序列上，N为序列xk和Xk的周期。这样定义对后面圆周卷积的处理是方便的。二、线性卷积与圆周卷积2.1线性卷积离散时间信号xn与hn的卷积定义为(5)mynxmhnm常记为(6)ynxnhn注1：卷积就是加权求和注2：公式(5)求和的范围是到，若序列xn和hn均为有限长序列，则当指标使得xm和hnm没有定义时，可令其值为0；注3：看个视频的例子线性卷积.swf注4：序列xn和hn均为有限长序列，长度分别为N和M，请考虑两序列卷积的结果yn的长度。注5：一个练习。1,2,3,4,5xn，1,2,3hn2.2圆周卷积对周期为N的信号xn和hn，其圆周卷积定义为(7)10Nmynxmhnm常记为(8)ynxnhn注1：上述圆周卷积定义在周期信号上注2：两信号的周期要相同注3：若考虑两长度为N的信号xn和hn，可将其圆周卷积定义如下(9)01NNmynxmhnm上式中Nr表示对r关于N求余所得的余数，0,1,2,,1NrN注4：视频圆周卷积.swf注5：一个练习。1,2,3,4xn，1,2,2,1hn注6：可以将(9)式写成矩阵形式2.3时域卷积定理以,,YXnHnn表示,h,xnnyn的离散Fourier变换，(10)ynxnhn则(11)YnXnHn证明：(12)11112200111222000011122200001111NmrnmkmNNNiiNNmkrrnmkmNNNiiNNkrmkrmrnNNNiiNNkrmynxmhnmXkeHreNNXkHreeNXkHreeN考虑102?krmNiNme(13)102,0,kkrmNiNmNkrer(14)01112220012011krmrnNNNiiNNkrmrnNiNrynXkHreeNXrHreN再看(14)右端恰为XrHr的IDFT。两端做离散Fourier变换即得(11)式。三、离散Fourier变换的应用离散Fourier变换的一个重要应用就是计算线性卷积。然而离散Fourier变换的卷积性质是跟圆周卷积相关的。考虑，如何利用离散Fourier变换计算线性卷积。四、考察离散Fourier变换的计算量考虑离散Fourier变换(15)111110000012112101001111NNNNNXx注意，上式是N个方程的计算，每计算一个Xn需要进行N次复数乘法和1N次复数加法。因此要完成全部运算共需2N次复数乘法和1NN复数加法。随着N增大，运算工作量将迅速增加，例如10N需要100次复数乘法运算，而当1024N时，就需要一百多万次复数乘法运算。按照这种规律，如果N较大的情况下，要求对信号进行实时处理，所需的运算速度就难以实现。幸运的是我们有FFT，它能够将DFT的计算复杂度将为logNN。五、快速Fourier变换快速Fourier（FastFourierTransform，FFT）是计算离散Fourier变换的一种特殊方法。FFT通过充分挖掘(15)中矩阵的性质，来简化计算。以下仅考虑2MN的情形。矩阵(16)111110000011201121NNNNN其中，(17)12iNWe下面我们就考察W和F的性质。(18)2kikNWek取某些值时，kW的表达是简单的242322024:1:43:42:1NiNNikNikNNikNNNkNkWeiNkWWNWekiee对称性(19)2NnknkWW周期性(20)NnknkWW矩阵F中，关于W的最高次幂是11NN。利用上述两个性质，F中最多有多少个不同的数？注意到在矩阵W中的某些数是非常简单的。例如，01W，21NW，实际上无需做乘法运算，在N较大的情况下，这一因素可使运算工作略有减少；再考虑到系数nkW的周期性和对称性，合理安排重复出现的相乘运算，将使计算工作量显著减少。（3）把N点DFT运算分解为两组2N点的DFT运算，然后求和下面证明这种分解是正确的，而且可以减少运算工作量对序列xn取2MN点DFT。把运算分为奇数和偶数两部分(21)101122212001122220011220202222222221221221NnkNnNNrkrkNNrrNNrkrkkNNNrrNNrkkrkNNNrrkNiNiNNNXkxnWxrWxrWxrWWxrWxrWWxrWGkWeWHekW其中(22)12021202221NrkNrNrkNrGkxrWHkxrW分别为xn的奇数部分和偶数部分构成的子列的DFT。也就是说，一个N点的DFT被分解为两个2N点的DFT。由周期性可知，(23)22NGkGkNHkHk又(24)22NNkkkNNNN将(23)，(24)代入(21)有(25)2,0,1,2,22122kNNkNkNNXkGkWHkNNNXkGkWHkkGkWHk或者(26)/2,//22/2/2,kkNkNNGkWHkkNXkGkNWHkNN伪代码,101000,/21040,1,,,1/2,4ditfft2,,#DFTof,,if1then#trivalsize1DFT,ditfft2,/basecaseelse#DFTof,2,2ditfft2,/2,2for,,#DFTof,,,1s2sNsNsNsNs2ssNssXxNsxxNXxXxxNsxsNskxxxXxxx/2/2/20to/21#combineDTFsoftwohalvesintofullDFTendforendifkkNkkNkkkNtWXXtNtXXWX考察算法的计算量(27),/2,k2/2kNkNXkGkWHkkNNXkkWkNGH例，8N蝶形运算单元每个蝶形运算单元需要一次复数乘法和两次复数加法。由Gk和Hk获得Xk的过程中，共包含2N个蝶形运算。因此需2N次复数乘法和N次复数加法。对于2MN的情况，需要把这种奇偶分解逐级进行下去。当2MN时，全部DFT运算可分解为M级流程图，其中每级都包含2N次复数乘法和N次复数加法，FFT的全部计算量为复数乘法运算：2log22NNMN复数加法运算：2logNMNN而直接DFT的计算量为复数乘法运算：2N复数加法运算：1NN下图为直接DFT和FFT算法中复数乘法次数与序列长度的关系曲线。以上是对序列长度是2N时的快速Fourier变换。当序列不是2MN时，也有相应的快速Fourier算法。请有兴趣同学自己查阅相关文献。FFTW是一个开源的FFT库，号称是世界上最快的FFT库。网址是六、FBP的频域滤波实现问题考虑滤波反投影重建问题I）若投影数据长度为L，卷积核应该多长？II）直接利用卷积公式的计算量是多少？III）利用FFT计算可以减少计算量？IV）补零的话应该补多少呢？V）如何得到频域的滤波器？卷积核做DFT，还是直接对加窗的斜坡滤波器采样，还是。。。