离散傅里叶变换的矩阵表示及其运算量

(k N )
WN 2
N
WNk WN2
WNk
j 2 N
e N 2
WNk
WN2k
j 2 2k
e N
j 2 k N
e 2
W
k N
2
4.2 基2时间抽选的FFT算法
4.2.1 算法推导
已经知道：X (k) DFT [x(n)] N 1 x(n)WNnk n0
令DFT的长度N=2M，M为正整数。
WNn
WNrN
N
WN2
N 1
N 1
N 1
2 n0
x(n)WN
n /
r 2
WNn
2 n0
x(n
N 2
)WN
n /
r 2
WNn
(1)
2
[ p(n)WNn
n0
]WN
nr /2
其中， p(n) x(n) x(n N ) 2
n 0,1,, N 1 2
显然，X(2r) 为g(n) 的N/2点DFT，X(2r+1) 为p(n)WNn 的
显然，直接计算N点的IDFT所需的复乘和复加的次数也是这 么多。当N足够大时，N2 ≈ N(N-1), 因此，DFT与IDFT的运 算次数与N2 成正比，随着N的增加，运算量将急剧增加，而 在实际问题中，N往往是较大的，因此有必要对DFT与IDFT 的计算方法予以改进。
4.1.2
W
nk因子的特性
n0
IDFT :
x(n)
1 N
N 1
X (k )WN nk
k 0
,
n 0,1, N 1(4.2)
式中
j2 。 下面要用矩阵来表示DFT关系。
WN e N
wenku.baidu.com
一般情况下，信号序列x(n) 及其频谱序列X(k) 都是用复数 来表示的，WN当然也是复数。因此，计算DFT的一个值X(k) 需要进行N次复数乘法（与1相乘也包括在内）和N-1次复数 加法；所以，直接计算N点的DFT需要进行N2 次复数乘法和 N(N-1) 复数加法。
的长度都为 N。现将 (4.37) 式中各序列作周期延拓，有
Y~(k ) Y~r (k ) jY~i (k )
(4.38)
由周期性有：
Y~r (k) Y~r (N k), Y~i (k ) Y~i (N k ),
Y~r (k) Y~r (N k) Y~i (k ) Y~i (N k )
则由DFT的线性有： Y (k) P(k) jG(k)
(4.36)
P(k) 和G(k) 这两个复序列的实部都是周期性的偶序列，而 其虚部都是周期性的奇序列。
对复序列Y(k) 又有 Y (k ) Yr (k ) jYi (k )
(4.37)
这里下标 r、i 分别表示实部和虚部，Y(k) 与其实部、虚部
r 2
WN
rN
2
g
(n)W
n N/
r 2
n0
其中，
g(n)
x(n)
x(n
N
)
2
n 0,1,, N 1 2
X
(2r
1)
N 1 2
x(n)WN2nrWN n
n0
N 1 2
x(n
n0
N 2
(2r1)(n N )
)WN
2
N
N
1 2
x(n)WNn/r2
n0
WNn
1 2
x(n
n0
N 2
)WN2nr
P(k
)
k 0,1,, N -1
其中，
N 1
N 1
2
2
G(k)
g
(r
)WN
r /
k 2
x(2r
)WN
r /
k 2
r 0
r 0
N 1
N 1
2
2
P(k)
p(n)WN
r /
k 2
x(2r
1)WN
r /
k 2
r 0
r 0
是由x(n)的偶数抽样点形成的DFT；而
是由x(n)的奇数抽样点形成的DFT。但是这两个式子并不完 全是N/2点的DFT，因为k的范围仍然是由0到N-1，因此，还 应该进一步考虑k由N/2到N-1范围的情况。
k 0,1, , N -1
令：
g(r) x(2r)
p(r)
x(2r
1)
r 0,1,, N 1 2
于是有：
N
N
N
N
1
1
1
1
2
2
2
2
X (k)
x(2r)WN2rk
x(2r
1)W
(2r N
1)
k
g
(r
)WN
rk /2
WNk
p(r
)WN
rk /2
r 0
r0
r 0
r 0
G
(k
)
W
k N
(4.39) (4.40)
现在将序列 Y~r (k )与 Y~i (k ) 作如下分解：
Y~r (k)
1 2
[Y~r
(k
)
Y~r (N
k )]
图 4.3 2 点 DFT 信号流图（蝶形图）
对于2点DFT，有：
T2
1 1
1 1
W21
1
1 1
所以2点DFT的运算只需一次加法和一次减法，这样的运算
叫做蝶形运算，这样的信号流图叫做蝶形图。
该算法每次分解都是将时域序列按奇偶分为两组，因此要求 N等于2的正整数幂，故将这种FFT算法叫做基2时间抽选法。
M c2
N 2
(log 2
N
1)
N 2
log 2
N 2
复数加法的次数为： Ac2
N 2 2 log 2 N N log 2 N
因此，FFT算法比直接计算DFT大大减少了运算量，尤其是
当N较大时，计算量的减少更为显著。比如，当N=1024时，
采用FFT算法时复数乘法的次数，不超过直接计算DFT时复
N
DFT和IDFT的快速算法的导出主要是根据
W
n因k 子的特
N
性。
W n(kN) N
WN nk
WN nN
WN nk
1．周期性：
对离散变量n有同样的周期性。
2．对称性：
[WNnk ]*
W nk N
WN (n)k
WN (Nn)k
或 3. 其它:
[WNnk ]*
W kn N
WN (k)n
WN (Nk)n
4.5 快速傅里叶反变换（IFFT）
IFFT是IDFT的快速算法。由于DFT的正变换和反变换的表 达式相似，因此IDFT也有相似的快速算法。可以用三种不 同的方法来导出IFFT算法。
方法1
设 x(n) DFT X (k )
n=0,1,┅,N-1
，
则有：
x(n)
1 N
N 1 X (k )WN nk ，
2. 同址计算 从图4.4可以看出，整个算法流图可以分为四段，（0）段
为倒序重排，后面3段为3(log28=3)次迭代运算：首先计算2 点DFT，然后将2点DFT的结果组合成4点DFT，最后将4点 DFT组合为8点DFT。因此，对于N点FFT，只需要一列存储 N个复数的存储器。
3. 运算量
观察图4.4可知，图4.3所示的蝶形图实际上代表了FFT的基 本运算，它实际上只包含了两次复数加法运算。一个 N(N=2M) 点的FFT，需要进行M=log2N次迭代运算。每次迭 代运算包含了N/2个蝶型，因此共有N次复数加法；此外，除 了第一次的2点DFT之外，每次迭代还包括了N/2次复数乘法 （即乘WN的幂）。因此，一个N点的FFT共有复数乘法的次 数为：
方法3 对DFT的反变换式取共轭，有：
x* (n)
1 N
N 1
[ X (k )WN
k 0
nk
]*
1 N
N 1
[X
k 0
* (k )WN nk
],
n 0,1,, N 1
与DFT的正变换式比较，可知完全可以利用FFT 的计算程序，只需要将X*(k) 作为输入序列，并将 最后结果取共轭，再除以N就得到x(n)。
现在令 k 0,1,, N 1，故对于后半段有： 2
G(k
N) 2
N 1 2
r(k N )
g(r)WN / 2 2
r 0
N 1
2
N
g(r)W N r 2
r 0
2
W N
2
rk
N 1
2
g
(r
)W
r N
k
G(k)
r 0
2
同理： 又知：
P(k N ) P(k) 2
kN
WN 2
WNk
图 4.1 N点DFT分解为两个N/2点的DFT (N=8)
乘次数的千分之五。
4.3 基2频率抽选的FFT算法 时间抽选法是在时域内将输入序列x(n)逐次分解为偶数点
子序列和奇数点子序列，通过求子序列的DFT而实现整个序 列的DFT。而频率抽选法则是在频域内将X(k)逐次分解成偶 数点子序列和奇数点子序列，然后对这些分解得越来越短的 子序列进行DFT运算，从而求得整个的DFT。当然，同样要 求N为2的正整数幂。
第4章 FFT
4.1 引言 4.1.1 离散傅里叶变换的矩阵表示及其运算量 DFT在数字信号处理中起着非常重要的作用， 这是与DFT
存在着高效算法， 即快速傅里叶变换(FFT) 分不开的。快速 运算的关键是减少运算量。
离散傅里叶变换对为：
DFT :
N 1
X (k ) x(n)WN nk , k 0,1, N 1 (4.1)
N/2点DFT。这样，一个N点的DFT分解为两个N/2点的DFT。
将分解继续下去，直到分解为2点的DFT为止。当N=8时，
基2频率抽选的FFT算法的整个信号流图如图4.6所示。
将图4.6与图4.4比较，可知频率抽选法的计算量与时间抽选 法相同，而且都能够同址计算。时间抽选法是输入序列按奇 偶分组，故x(n)的顺序要按倒序重排，而输出序列按前后分 半，故X(k) 的顺序不需要重排；频率抽选法则是输出序列按 奇偶分组，故X(k) 的顺序要按倒序重排，而输入序列按前后 分半，故x(n) 不需要重排。
4.2.2 算法特点 1. 倒序重排 这种FFT算法的每次分解都是将输入序列按照奇偶分为两组，
故要不断地将每组输入数据按奇偶重排，直到最后分解为2 点的DFT，输入数据才不再改变顺序。这样做的结果，使得 作FFT运算时，输入序列的次序要按其序号的倒序进行重新 排列。
现在将图4.4中输入序号以及重排后的序号按二进制写出如 下（注：下标“2”表示二进制数）。可以看出，将输入序 号的二进制表示（n2n1n0）位置颠倒，得到（n0n1n2），就 是相应的倒序的二进制序号。因此，输入序列按倒序重排， 实际上就是将序号为（n2n1n0）的元素与序号为（n0n1n2） 的元素的位置相互交换。
图 4.2 N/2 点的DFT 分解为两个N/4点的DFT (N=8)
综上所述，可以得到：
X (k) G(k) WN k P(k)
X
(k
N 2
)
G(k)
WN k P(k)
k 0,1,, N 1 2
其中G(k)、P(k) 分别是x(n)的偶数点和奇数点的N/2点DFT。
这样，我们就将一个N点的DFT分解成了两个N/2点的DFT， 由于DFT的运算量与其点数的平方成正比，因此使运算量减 少了。但是，还应该将每一个N/2点的DFT再分解为两个N/4 点的DFT，如此下去，直到分解为2点的DFT为止，总共需 要进行log2N-1=log2(N/2)次分解。
方法2
将DFT的正变换式：
N 1
X (k ) x(n)WN nk
与其反变换式：
x(n)
1 N
N 1 n0
X (k )WN nk
k 0
比较，很容易知道，可以利用FFT算法的程序来 计算IFFT，只需要将X(k) 作为输入序列，x(n) 则 是输出序列，另外将因子WNk 变为WN-k， 当然， 最后还必须将输出序列的每个元素除以N。
设 r 0,1,, N 1 ，则可以分别表示出k为偶数和奇数时 2
的X(k)。
N 1
N 1
X (2r)
N 1
x(n)WN 2rn
n0
2
x(n)WN 2nr
n0
2
x(n
n0
N 2
)WN
2(n
N 2
)r
N 1
N 1
N 1
2
x(n)WN
n /
r 2
n0
2
x(n
n0
N 2
)W
N
n /
k 0
根据基2 FFT的时间抽选法的第一次分解的结果：
X (k) G(k) WN k P(k)
X
(k
N 2
)
G(k)
WN
k P(k)
k 0,1,, N 1 2
可以得到：
G(k)
P(k)
1[X
2
1 2
WNk
[
(k) X (k)
X
(k X (k
N )] 2 N )]
2
4.7.1 两个长度相同的实序列 可以将两个长度相同的实序列组合成一个复序列来进行
FFT运算，从而一次完成这两个实序列的FFT，减少了总 的计算量。
设p(n) 和g(n) 是两个长度均为N的实序列，并设：
又设
y(n) p(n) jg(n)
p(n) DFT P(k) ， g(n) DFT G(k) ，y(n) DFT Y (k )
k 0,1,, N 1 2
图 4.8 由 X(k)、X(k+N/2) 得到 G(k)、P(k)
再由N/2点的DFT求得N/4点的DFT，依次类推下去，就可推 到求出x(n)的各点，如图4.9所示。将此流图与图4.4比较，相 当于整个流向反过来，此外，因子WNk成为WN-k，还增加了 因子1/2。但是，图4.9的IFFT算法不能直接利用按照图4.4编 写的FFT算法程序，却可以利用图4.6的频率抽选FFT算法的 程序，只需要将X(k)作为输入序列，因子WNk变为WN-k，并 且将最后所得的输出序列的每个元素都除以N。