叶东明:二维FFT的程序实现及其应用
F(u,v)?????????f(x,y)exp[?i2?(ux?vy)]dxdy (1-6-1)
为f(x,y)的二维Fourier变换2D-FT(2 dimensional Fourier transform或称为二维谱函数。
对应地称
f(x,y)?????????F(u,v)exp[i2?(ux?vy)]dudv (1-6-2)
为(1-6-1)的二维Fourier变换的逆变换,又称为反演公式。 1.7 二维离散Fourier变换
定义1.7.1
对N1?N2二维复序列A(k1,k2),称
x(j1,j2)?N1?1N2?1k1?0k2?0??A(k,k12?k1j1?k2j2 )WNW N12 (1-7-1)
j1?0,1,?,N1-1, j2?0,1,?,N2-1为二维离散Fourier变换2D-DFT(2 dimensional discrete Fourier transform)。对应
地称
1A(k1,k2)x(j1,j2)?N1N2N1?1N2?1j1?0j2?0??x(j,j12k1j1k2j2 )WNW N12 (1-7-2)
k1?0,1,?,N1-1, k2?0,1,?,N2-12?i2?i为二维离散Fourier变换(1-7-1)的逆变换。这里WN1?exp(),WN2?exp()。
N1N22 DFT快速计算(FFT)算法分析及其程序实例
2.1 DFT的快速算法
本文将长度为N的序列{Ak}变为长度为N的序列{xj}的DFT及其逆变换形式
取为
?jkxj??AkWN j?0,1,?,N-1 (2-1-1)
k?0N?11N?1Ak??xjWNjk k?0,1,?,N-1 (2-1-2)
Nj?0为了可以用统一的程序处理DFT及其逆变换,我们将(2-1-2)改写为
1N?1?jkAk??xjWN k?0,1,?,N-1 (2-1-3)
Nj?0两个只差一个常数因子,求xj和Ak可以使用相同的程序。
Danielson和Lanczos两人提出一种将长度为N的序列{Ak}的DFT分为两个长度
11
集美大学学士学位论文
为
N的序列的DFT。其中一个序列由原序列{Ak}的偶数项构成,记其DFT为x0j,2而另一个序列由{Ak}的奇数项构成,记其DFT为xej。这时原序列的DFT为
?jkxj??AkWN k?0N?12k?0N?12k?0N?1?2kj?(2k?1)j ??A2kWN ? ?A2k?1WNN?12k?0N?12k?0 (2-1-4)
?kj?j?kj ??A2kWN ? WN ?A2k?1WN22?j0 ?xej?WNxj , j?0,1,?,N-1 2xN2?j??AkWNk?0N?1k?0N?1?(N?j)k2 ?jk ? ?(?1)kAkWN N?12k?0N?12k?0 ??A2k2k?2kj?(2k?1)j (2-1-5) WN ??A2k?1WNN?12k?0 ??Ak?0N?12?kj?j?kjWN ? WN ?A2k?1WN22?j0 ?xej?WNxj , j?0,1,?,N-1 2如此将长度为N的序列{Ak}的DFT{xj}用它的偶序列{A2k}和奇序列{A2k?1}的DFT,
e即{x0由于{xe{x0j}和{xj}表示。j},j}长度为
N可以分别计算。而长度为N的序列{xj}2可按(2-1-4),(2-1-5)求得。 为
只要我们预先取N?2r,当我们将长度为N的序列{Ak}的DFT转变为求两个长
N的序列{A2k},我们可再对{A2k},{A2k?1}的DFT后,{A2k?1}分别作同样的转变,2N即转变为求长度的序列的DFT,这样该算法便可以递归地使用。
42.2 基2FFT算法分析
现在我们讨论N?2m,m为整数的FFT算法。将j,k写成m位2进制数,即
j?jm?12m?1?jm?22m?2???j121?j0?(jm?1,jm?2,?,j0)k?km?12
m?1?km?22m?2???k12?k0?(km?1,km?2,?,k0)12
1 (2-2-1)
叶东明:二维FFT的程序实现及其应用
其中jp,kp只取0或1,p?0,1,?,m?1。记
xj?x(jm?1,jm?2,?,j0)Ak?A(km?1,km?2,?,k0) j,k?0,1,?,N?1 (2-2-2)
j?k?(jm?1,jm?2,?,j0)(km?1,km?2,?,k0) ?j0km?12由{Ak}的DFT定义得
m?1?(j1,j0)km?22m?2???(jm?1,jm?2,?,j0)k0(modN) (2-2-3)
xj?x(jm?1,jm?2,?,j0) ? ?其中
A1(j0,km?2,?,k0) ?km?1?0k0?0k1?011???11km?1?01?j(km?1,km?2,?,k0)A(k,k,?,k)W?m?1m?20N1 (2-2-4)
k0?0k1?0????A(j,k10km?2?0m?2?j(0,km?2,?,k0),?,k0)WN j0,j1,?,jm?1?0,1 ?A(k1m?1?j(km?1,0,?,0),km?2,?,k0)WN (2-2-5)
?(j0,0,?,0) ?A(0,km?2,?,k0)?A(1,km?2,?,k0)WN j0,km?2,?,k0?0,1
显然A1(j0,km?2,?,k0)的第一项的因子为1,只有第二项需要作一次复数乘法。所
以计算A1(j0,km?2,?,k0)需要2m次乘法。而A1的序号(j0,km?2,?,k0)的左第一位为
j0,将j0左移m-1位和WN的指数相符。
继续对xj分解得
xj?x(jm?1,jm?2,?,j0) ? ?其中
A2(j0,j1,km?3,?,k0) ?km?2?0k0?0k1?011???11km?1?01?j(0,km?2,?,k0)A(j,k,?,k)W j0,j1,?,jm?1?0,1 (2-2-6) ?10m?20N1k0?0k1?0????A(j,j,k201km?3?0m?3?j(0,0,km?3,?,k0),?,k0)WN ?A(j,k101m?2?j(0,km?2,0,?,0),?,k0)WN (2-2-7)
?(j1,j0,0,?,0) ?A1(j0,0,km?3,?,k0)?A1(j0,1,km?3,?,k0)WN j0,j1,km?3,?,k0?0,1
13
集美大学学士学位论文
计算A2(j0,j1,km?3,?,k0)需要2m次复数乘法,同时A2的序号(j0,j1,km?3,?,k0)的左两位为(j0,j1),将其逆排后左移m-2位和WN的指数相符。继续分解则有
Al ( l?1,2,?,m) 。若记
A0(km?1,km?2,?,k0)?A(km?1,km?2,?,k0) (2-2-8)
则有
Al(j0,j1,?,jl?1,km?l?1,?,k0) ?km?l?0?A1l?1?j(0,?,km?l,0,?,0)(j0,?,jl?2,km?l,?,k0)WN ?Al?1(j0,?,jl?2,0,km?l?1,?,k0)?(jl?1,?,j0,0,?,0) ?Al?1(j0,?,jl?2,1,km?l?1,?,k0)WN (2-2-9)
j0,?,jl?1,km?l?1,?,k0?0,1 l?1,2,?,m
对某一确定的l,计算Al(j0,j1,?,jl?1,km?l?1,?,k0)需要2m次复数乘法,同时Al的
序号(j0,j1,?,jl?1,km?l?1,?,k0)的左l位为(j0,j1,?,jl?1),将其逆排后左移m?l位和
WN的指数相符。
当完成l?m步计算后得
xj?x(jm?1,jm?2,?,j0)?Am(j0,j1,?,jm?1) j0,j1,?,jm?1?0,1 (2-2-10)
可以看到,xj的序号(jm?1,jm?2,?,j0)逆排后得(j0,j1,?,jm?1)恰好为Am的序号。计算每一Al需要2m次复数乘法,所以由{Ak}求得{xj}共需要m?2m次复数运算。即当
N?2m时,求长度为N复序列{Ak}的DFT基2FFT算法运算总次数为
N (2-2-11) m2m?Nlog2
从以上分析可得: 基2FFT算法实现的关键在于下标的组织;
基2FFT算法的必要条件为N?2m,m为整数。
对不满足N?2m条件的序列{Ak},我们可用零扩充的方法,即取最小的m,满
足2m?1?N?2m,则可令:
14
叶东明:二维FFT的程序实现及其应用
?Ak 0?k?NAk?? (2-2-12) m?0 N?k?2当回到时间域后,我们可以将扩充的零元素删除。 2.3 复序列基2FFT算法
设{Ak}是长度为N的复序列,其中N?2m。 算法2.3.1 记A0?A: (1) 计算WNj?cos2?j2?j?isin, j?0,1,?,N-1 ; NN基2FFT正变换
(2) 对l?1,2,?,m计算Al:
a.记k?(km?l?1,?,k0) ,j?(j0,j1,?,jl?1); b.逆排(j0,j1,?,jl?2)得(jl?2,?,j0);
c.左移(jl?2,?,j0)m?l位得(0,jl?2,?,j0,0,?,0),由此查得WNd.对j0,j1,?,jl?2,km?l?1,?,k0?0,1计算
?(0,jl?2,?,j0,0,?,0);
Al(j0,j1,?,jl?2,0,km?l?1,?,k0) ?Al?1(j0,?,jl?2,0,km?l?1,?,k0)?(0,jl?2,?,j0,0,?,0) ?Al?1(j0,?,jl?2,1,km?l?1,?,k0)WN
Al(j0,j1,?,jl?2,1,km?l?1,?,k0) ?Al?1(j0,?,jl?2,0,km?l?1,?,k0)?(0,jl?2,?,j0,0,?,0) ?Al?1(j0,?,jl?2,1,km?l?1,?,k0)WN
(3) 对j0,j1,?,jm?1?0,1完成: a. 逆排(jm?1,?,j0)得(j0,j1,?,jm?1) b. 令
x(jm?1,jm?2,?,j0)?Am(j0,j1,?,jm?1)
正变换完成,由{Ak}得其DFT{xj}。
算法2.3.2
基2FFT逆变换
置Ak?xk,k?0,1,?,N?1: (1) 计算WNj?cos2?j2?j?isin, j?0,1,?,N-1 NN15
(2) 置ImAk??ImAk,k?0,1,?,N?1