# FFT算法研究及基2-FFT算法的C语言实现

I

Abstract

Abstract
Discrete Fourier Transform (DFT) is often used to calculate the signal processing to obtain frequency domain signals. DFT algorithm can get the frequency domain characteristics of the signal, because the algorithm is computationally intensive, long-time use, the computer is not conducive to real-time signal processing. So DFT since it was discovered in a relatively long period of time is not to be applied to the actual projects until a new fast discrete Fourier calculation method --FFT is found in discrete Fourier transform was able to actually project has been widely used. FFT is not a new way to get the frequency domain, but the discrete Fourier transform of a fast algorithm. Fast Fourier Transform (FFT) is a digital signal processing important tool that the time domain signal into a frequency-domain signal processing. matched filtering has important applications. FFT is a discrete Fourier transform (DFT) is a fast algorithm, it can be a signal from the time domain to the frequency domain. Some signals in the time domain is not easy to observe the characteristics of what is, but then if you change the signal from the time domain to the frequency domain, it is easy to see out of. This design is required to be familiar with the basic principles of FFT algorithm, based on the preparation of C language program to achieve FFT algorithm real number sequence. Keywords: FFT algorithm, C language compiler to achieve

II

Abstract

III

1 引言
1.1 课题背景

1.2 FFT 算法的现状

1 引言

1.3 研究内容

1.4 论文的研究成果

2

1 引言

2 数字信号处理综述
2.1 数字信号理论

2 数字信号处理综述

2.2 数字信号处理的实现

2.3 数字信号处理的应用及特点

4

2 数字信号处理综述

5

2 数字信号处理综述

3 基本理论
3.1 FFT 算法的基本概念

3.1.1 离散傅里叶变换（DFT）

2? nk N

X (k ) ? DFT [ x(n)] ? ? x(n)e
n ?0

N ?1

?j

nk ? ? x(n)WN n ?0

N ?1

6

4 基-2FFT 算法的 C 语言实现及仿真

x(n) ? IDFT [ X (k )] ?

2? N ?1 j nk 1 N ?1 ? nk N X ( k ) e ? x(k )WN ? ? N k ?0 k ?0

1 N

，所以有关

3.1.2 快速傅里叶变换（FFT）

kn X ( k ) = ? x(n)WN RN ( k ) n ?0 N ?1

(k ? N )n kn ( n? N ) k （1）周期性： WN ? WN ? WN

kn N

( k ? N / 2) k （2）对称性： WN ? ?WN

X (2) ? ? x(n)W42 n ? x(0)W40 ? x(1)W42 ? x(2)W44 ? x(3)W46
n ?0 3

? [ x(0) ? x(2)]W40 ? [ x(1) ? x(3)]W42 ＝[ x(0) ? x(2)]－[ x(1) ? x(3)]W
0 4

(周期性)

(对称性)

4 基-2FFT 算法的 C 语言实现及仿真

3.2 FFT 算法的分类
DFT 从基本上可以分为两类分解方法：一类是把时间序列 x(n) 进行连续的分 解，于是得到的 FFT 算法可以称之为按时间抽取的算法，另一类是把傅立叶变换 序列 X(k) (k 为频率基准)进行连续的分解 ，叫做按频率抽取的算法。对于每一种 算法，根据基本的蝶形运算的构成，又可把 FFT 算法分为基 2、基 4、基 8 以及 任意因子等的 FFT 算法，FFT 算法对于不同基所需要的计算量是有所不同的，之 所以说有差异是指的是并没有数量级的差别。表 3.1 列出了 N=4096 时各种基的 FFT 算法所需运算量次数的的比较。

3.2.1 按时间抽取（DIT）的 FTT

? x(2r ) ? x1 (r ) ? ? x(2r ? 1) ? x2 (r )

X (k ) ? DFT[ x(n)] ?

r ? 0,1,?,

N ?1 2

? x(n)W
n ?0

N ?1

kn N

8

4 基-2FFT 算法的 C 语言实现及仿真

?

n ?0 n为偶数
N / 2?1 r ?0

? x(n)WNkn＋

N ?1

n ?0 n 为奇数

? x(n)W
N / 2?1 r ?0

N ?1

kn N

2 rk ＝ ? x(2r )WN ?

? x(2r ? 1)W
k N N / 2 ?1 r ?0

( 2 r ?1) k N

＝ ? x1 (r )W
r ?0
N / 2 ?1 r ?0

N / 2 ?1

2 rk N

?W

?x
2

2

2 rk (r )WN

rk k ＝ ? x1 (r )WN / 2 ? WN

N / 2 ?1 r ?0

? x (r )W

rk N /2

2 rk rk （因为 WN ? WN /2 ）

k ? X 1 (k ) ? WN X 2 (k )

rk X 1 ( k ) ＝ ? x1 (r )WN /2 ? r ?0 N / 2 ?1 N / 2?1 r ?0

? x(2r )W
rk N /2

rk N /2

,0 ? k ?

N ?1 2

rk X 2 ( k ) ＝ ? x2 (r )WN /2 ? r ?0

N / 2 ?1

N / 2?1 r ?0

? x(2r ? 1)W

,0 ? k ?

N ?1 2

k X ( k ) ? X 1 (k ) ? WN X 1 ( k ) 和 X 2 (k ) 只有 N/2 个点 分析： （ k ? 0,1,?, N 2 -1 ） ， X 2 (k )

X 2 ( k ) 和 X (k ? N / 2) 的关系，其中 k ? 0,1,?, N 2 ? 1 。
k 是由式子 X ( k ) ? X 1 (k ) ? WN X 2 (k ) 可得

k?N / 2 X (k ? N / 2) ? X 1 (k ? N / 2) ? WN X 2 (k ? N / 2)

k X (k ? N / 2) ＝ ? X 1 (k ) ? WN X 2 (k ) ， k ? 0,1,?,

N ?1 2

k ? ? X (k ) ? X 1 (k ) ? W N X 2 (k ) ? k ? ? X (k ? N / 2) ? X 1 (k ) ? W N X 2 (k )

k ? 0,1,?,

N ?1 2

a
k a ? WN b

4 基-2FFT 算法的 C 语言实现及仿真

b

W

k N

-1

k a ? WN b

DFT 通过这样的分解以后，每一个 N 2 点的 DFT 只需要 ( N 2 ) 2 ? N 4 次复数乘法
2

2 2

2

2

2

? x1 (2l ) ? x3 (l ) ? ? x1 (2l ? 1) ? x4 (l )

l ? 0,1,?,

N ?1 4

2lk X 1 ( k ) ＝ ? x1 (2l )WN /2 ? l ?0

N / 4 ?1

N / 4?1 l ?0

? x (2l ? 1)W
1 N / 4?1 l ?0

( 2l ?1) k N /2

lk k ＝ ? x3 (l )WN / 4 ? WN / 2 l ?0

N / 4?1

?x

4

lk (l )WN /4

k ? X 3 (k ) ? WN / 2 X 4 (k )

k ? 0,1,?,

N ?1 4

k 2k 将系数统一为以Ｎ为周期，即 WN / 2 ? WN ，可得
2k ? ? X 1 (k ) ? X 3 (k ) ? W N X 4 (k ) ? 2k ? ? X 1 (k ? N / 4) ? X 3 (k ) ? W N X 4 (k )

k ? 0,1,?,

N ?1 4

10

4 基-2FFT 算法的 C 语言实现及仿真

“时间抽取法”就是在每一次分解时，其步骤是按有限长序列顺序在时域上的 奇偶进行抽取。 根据上面流图的分析， N ? 2 M ，总共需要进行 M 次分解，于是构成了从 x(n) 至 X(k)的 M 级运算过程。每一级的计算都是由 N/2 个蝶形运算构成的，因此每一 级计算都需要 N/2 次复数乘法和 N 次复数加法，因此按时间抽取的 M 级运算后总 共需要 复数乘法次数： m F ?
N N ? M ? log 2 N 2 2

11

4 基-2FFT 算法的 C 语言实现及仿真

3.2.2 按频率抽取（DIF）的 FTT
12

4 基-2FFT 算法的 C 语言实现及仿真

X (k ) ? DFT[ x(n)] ?
( N / 2 ) ?1 n ?0

? x(n)W
n ?0

N ?1

kn N

?

? x(n)WNkn＋ ? x(n)WNkn
n? N / 2
nk N

N ?1

＝ ? x(n)W
n ?0 N / 2 ?1 n ?0

N / 2 ?1

?

N / 2?1 n ?0

? x(n ? N / 2)W

( n ? N / 2) k N

( N / 2) k nk ＝ ? [ x ( n) ? W N x(n ? N / 2)] WN

N/2 ( N / 2) k 因为 WN ? ?1,WN ? (?1) k ，k 为偶数时 (?1) k ? 1，k 为奇数时 (?1) k ? ?1 ，

nk X (k )＝ ?[ x(n) ? (?1) k x(n ? N / 2)] WN n ?0 N / 2 ?1

2 nr X (2r )＝ ? [ x(n) ? x(n ? N / 2)] WN n ?0

N / 2 ?1

?

N / 2 ?1 n ?0

?[ x(n) ? x(n ? N / 2)]W
N / 2 ?1 n ?0

nr N /2

( 2 r ?1) n X (2r ? 1)＝ ? [ x(n) ? x(n ? N / 2)] WN

?

N / 2 ?1 n ?0

?[ x(n) ? x(n ? N / 2)]W

n N

W

nr N /2

? x1 (n) ? x(n) ? x(n ? N / 2) 令? n ? x 2 (n) ? [ x(n) ? x(n ? N / 2)]W N

n ? 0,1,?, N / 2 ? 1

nr X (2r )＝ ? x1 (n)] WN /2 n ?0 N / 2 ?1

rn X (2r ? 1)＝ ? x2 (n)WN /2 n ?0

N / 2 ?1

4 基-2FFT 算法的 C 语言实现及仿真

a

a?b
n (a ? b)WN

b -1

W

n N

14

4 基-2FFT 算法的 C 语言实现及仿真

x(n) ? IDFT[ X (k )] ?

1 N ?1 ? nk X (k )WN , n ? 0,1,?, N ? 1 ? N k ?0

nk X (k ) ? DFT[ x(n)] ? ? x(n)WN , k ? 0,1,?, N ? 1 n ?0

N ?1

3.2.3 快速傅里叶分裂基 FFT 算法

15

4 基-2FFT 算法的 C 语言实现及仿真

n ? pn1 ? n0 ?

N n1 ? n0 , 4

0 ? n1 ? 3, 0 ? n0 ?

N ?1 4

kn X (k ) ? ? x(n)WN n ?0 N ?1

?

N /4?1 3

??
n ?0 n1 ?0
N /4?1

x(

N k ( n1 ? n0 ) N n1 ? n0 )WN 4 4

?

?
n0

kn0 WN ? x( n1

3

N n1 ? n0 )W4kn1 ] 4 N k N )W4 ? x(n0 ? )W42 k 4 2

?

N/ 4 ? 1 n0 ?0

?

[ x(n)W40 ? x(n0 ?

? x(n0 ?

3N 3k 3k kn0 )W4 WN ]WN 4

X (k ) ? X (4k1 ? k0 ) ?
N /4 ?1 n0 ? 0

?

[ x(n0 ) ? x(n0 ?

N (4 k1 ? k0 ) )W4 4

? x(n0 ? ?
N /4 ?1 n0 ? 0

?

N 81 ? 2 k0 3N 3(4 k1 ? k0 ) (4 k1 ? k0 ) n0 )W4 ? x(n0 ? )W4 ]WN 2 4 N N [ x(n0 ) ? x((n0 ? )W4k0 ? x(n0 ? )W42 k0 4 2 3N 3k0 (4 k1 ? k0 ) n0 )W4 ]WN 4

? x(n0 ?

16

4 基-2FFT 算法的 C 语言实现及仿真

N /2 ?1 N N ? X (2 k ) ? [ x(n) ? x(n ? )]WN2 kn , 0 ? k ? ? 1 ? ? 2 2 n ?0 ? N /4 ?1 ? ? N N 3N ? 4 kn )]WNn ?WN ? X (4k ? 1) ? ? ?[ x(n) ? jx(n ? ) ? x(n ? ) ? jx(n ? 4 2 4 ? n ?0 ? ? ? N ? 0 ? k ? ?1 ? 4 ? N /4 ? 1 ? ? N N 3N 3n ? 4 kn )]WN ?WN ? X (4k ? 3) ? ? ?[ x(n) ? jx(n ? ) ? x(n ? ) ? jx(n ? 4 2 4 ? n?0 ? ? ? N 0 ? k ? ?1 ? 4 ? ?

N /2 ?1 N N ? X (2 k ) ? [ x(n) ? x(n ? )]WN2 kn , 0 ? k ? ? 1 ? ? 2 2 n ?0 ? N /4 ? 1 ? ? N N 3N ? 4 kn )]WNn ?WN ? X (4k ? 1) ? ? ?[ x(n) ? jx(n ? ) ? x(n ? ) ? jx(n ? 4 2 4 ? n ?0 ? ? ? N ? 0 ? k ? ?1 ? 4 ? N /4 ? 1 ? ? N N 3N 3n ? 4 kn )]WN ?WN ? X (4k ? 3) ? ? ?[ x(n) ? jx(n ? ) ? x(n ? ) ? jx(n ? 4 2 4 ? n?0 ? ? ? N 0 ? k ? ?1 ? 4 ? ? N N ? 0 ? n ? ?1 ? x2 (n) ? x(n) ? x(n ? 2 ), 2 ? N N ? x1 (n) ? [ x(n) ? jx(n ? ) ? x(n ? ) ? jx(n ? 3N )]W n N ? 4 4 2 4 ? ? x 2 (n) ? [ x(n) ? jx(n ? N ) ? x(n ? N ) ? jx(n ? 3N )]W 3n N ? 4 4 2 4 ? N ? 0 ? n ? ?1 ? 4

N /2 ?1 N ? X (2 k ) ? x2 (n)WN2 kn ? DFT [ x2 (n)], 0 ? k ? ? 1 ? ? 2 n ?0 ? N /4 ?1 ? N 1 4 kn 1 ? X (4k ? 1) ? ? x4 (n)WN ? DFT [ x4 (n)], 0 ? k ? ? 1 4 n ?0 ? N /4 ? 1 ? N 2 4 kn 2 ? X (4k ? 3) ? ? x4 (n)WN ? DFT [ x4 (n)], 0 ? k ? ? 1 4 n ?0 ?

17

4 基-2FFT 算法的 C 语言实现及仿真

3.2.4 N 为组合数的 FFT——混合基算法

?n ? n1Q ? n0 ? ?k ? k1P ? k0
n0 , k1 ，分别为0, 1,…，Q-1， n1 , k0 ，分别为0, 1,…，P-1

N 点 DFT 可以重新写成为

X (k ) ? X (k1P ? k0 ) ? X ? k1 , k0 ?
kn ? ? x(n)WN n ?0 N ?1

( k1P ? k0 )( n1Q ? n0 ) ? ? ? x ? n1Q ? n0 ?WN n0 ?0 n1 ?0

Q ?1 P ?1

X ? k1 , k0 ? ? ?

n0 ? 0 n1 ? 0 Q ?1 P ?1

? ? x ? n , n ?W
1 0 1 0

Q ?1 P ?1

k1n1PQ N

k0 n1Q k1n0 P k0 n0 WN WN WN

n0 ? 0 n1 ? 0

? ? x ? n , n ?W

k0 n1Q N

k1n0 P k1n0 WN WN

18

4 基-2FFT 算法的 C 语言实现及仿真
k0 n1Q N

W
Q ?1

?W

k0 n1 P

X ? k1 , k0 ? ?

? P ?1 ? kn ? k n k n ? ?[ ? x ? n1 , n0 ?WP 0 1 ]WN 0 0 ?WQ 1 0 ? n0 ? 0 ? ? n1 ?0 ? ?
P ?1

kn 令 X1 ? k0 , n0 ? ? ? x ? n1 , n0 ?WP 0 1 n1 ?0
Q ?1

X ? k1 , k0 ? ?

n0 ?0

? [X ? k , n ? ? W
1 0 0

k0 n0 N

k1n0 ]WQ

X1? ? k0 , n0 ? ? X1 ? k0 , n0 ? ?WNk0n0
X ? k1 , k0 ? ?
n0 ?0

? X ? ? k , n ? ?W
1 0 0

Q ?1

k1n0 Q

x(0,0) ? x(0), x(0,1) ? x(1), x(0,2) ? x(2), x(0,3) ? x(3)
x(1,0) ? x(4), x(1,1) ? x(5), x(1,2) ? x(6), x(1,3) ? x(7) x(2,0) ? x(8), x(2,1) ? x(9), x(2,2) ? x(10), x(2,3) ? x(11)

(2)求 Q 个 P 点的 DFT

X1 (k 0, n0 ) ? ? x(n1 , n0 )W3k0n1
n1 ?0

2

(3) X1 (k0 , n0 ) 乘以得到 X1' (k0 , n0 ) 。 (4)求 P 个 Q 点的 DFT，参变量是 k 0
3

X 2 ? k0 , k1 ? ? ? X1' ? k0 , n0 ?W4k1n0
n0 ?0

(5) 将 x2 (k0 , k1 ) 通过 X (k0 ? k1P) 恢复为 X ( k ) （1）求 Q 个 P 点 DFT 需要 QP2次复数乘法和 Q· P· (P-1)次复数加法； （2）乘 N 个 W 因子需要 N 次复数乘法； （3）求 P 个 Q 点 DFT 需要 PQ2 次复数乘法和 P· Q（Q-1)次复数加法。 总的复数乘法量:QP2+N+PQ2=N(P+Q+1)； 总的复数加法量:QP(P-1)+PQ(Q-1)=N(P+Q-2)； 上述分解原则可推广至任意基数的更加复杂的情况。 例如,如果 N 可分解为 m 个质数因子 p1 , p2 ,?, pm 即 19

4 基-2FFT 算法的 C 语言实现及仿真

N ? p1 , p2 ,? pm

p1 ? p2 ? ? ? pm ? 2 为基2FFT 算法。

X1 (k0 , n1 , n0 ) ? ? x(n2 , n1 , n0 )W4n2k0 ,
n2 ?0 3

0 ? k0 ? 3

? X 1 (0, n1 , n0 ) ? ?W40 W40 W40 W40 ? ? x(0, n1 , n0 ) ? ? X (1, n , n ) ? ? 0 ? 1 2 3?? ? 1 1 0 ? ? ?W4 W4 W4 W4 ? ? x(1, n1 , n0 ) ? ? X 1 (2, n1 , n0 ) ? ?W40 W42 W44 W46 ? ? x(2, n1 , n0 ) ? ? ? ? 0 ? 3 6 9?? ? X 1 (3, n1 , n0 ) ? ? ?W4 W4 W4 W4 ? ? ? x(3, n1 , n0 ) ? ?1 1 1 1 ? ? x(0, n1 , n2 ) ? ?1 ? j ?1 j ? ? x(1, n , n ) ? 1 2 ? ?? ?? ?1 ?1 1 ?1? ? x(2, n1 , n2 ) ? ? ? ?? ?1 j ?1 ? j ? ? x(3, n1 , n2 ) ?

3.3 傅里叶变换的应用
FFT 算法可以分析处理每一个可以使用傅里叶变换来进行分析、处理的领域， FFT 算法会对这些地方用数字计算处理然后利用软件实现其所表现出来的功能。 其表现形式为包括卷积分的积分，它是傅里叶变换实现近似处理的理论基础。 不仅在信号、通信里会用到它，还在需要仿真、分析一些系统，对语音功率谱的 估计时会用到。

20

4 基-2FFT 算法的 C 语言实现及仿真

4 基-2FFT 算法的 C 语言实现及仿真

4.1 码位倒序

21

i<j? Y x(i)和 x(j)交换

k?N/2

k<j+1? N j?j+k

Y

j?j-k k?k/2

4 基-2FFT 算法的 C 语言实现及仿真

} }

4.2 蝶形运算
FFT 的 C 语言编程只需用 3 层循环即可实现，这是由 N 点的 FFT 从左到右共 有 log2N 级蝶形，每级有 N/2L 组，每组有 L 个，具体的三层是：完成整个 FFT 共 log2N 级每一级的蝶形运算的最外层循环， 完成每一组的蝶形运算 （每一级有 N/2L 组）的中间层循环，完成每一组有 L 个的单独 1 个蝶形运算的最内层循环。 第一层循环：运算的级数进行控制是第一层循环完成的，需要 m 级计算 N=2^m。 第二层循环： 根据乘数进行控制第二层循环， 这是由有 2^(L-1)个蝶形因子 （乘 数）第 L 级决定的。第三层循环里面的每一个蝶形因子都必须要执行一次，每一 级要进行 2^(L-1)次循环计算，第三层循环此时是在第二层循环控制下。 第三层循环：第三层循环在第二层循环确定某一乘数后要将本级中每个群中 具有这一乘数的蝶形计算一次，因为同一级内不同群的乘数分布相同，第 L 级共 有 N/2^L 即 2^ （n-L） 个群， 即第三层循环每执行完一次要进行 N/2^L 个碟形计算。 （执行不同 group 中具有相同 W 的蝶形运算） 可以得出结论：第三层循环在每一级中需要完成 N/2^L 个碟形计算；第三层 循环在第二层循环下进行 2^(L-1)次， 碟形计算在第二层循环完成时共进行 2^(L-1) *N/2^L=N/2 个。第 L 级的计算在第二、第三层循环得以完成。 实序列 256 点 FFT 运算的主程序代码基 2FFT 算法的 C 语言程序编制，是在 利用 C 语言软件调用上述 FFT 算法实现的，这是在设计要求下完成的。 在这里 x（n） 是一个有限长序列，共有 256 点，由于 N ? 256 ? 28 ，M=8,所以 总共有 8 级蝶形运算，每级有 N 2 ? 128个蝶形，共需要 1024 次复数乘法和 2048 次 复数加法。 用 C 语言实现 256 点实序列基 2FFT 的源程序见附录 A，仿真见附录 B。

23

4 基-2FFT 算法的 C 语言实现及仿真

[1] 郑君里， 应启珩， 杨为理． 信号与系统． 第二版.北京： 高等教育出版社， 2000. [2] 管致中，夏恭恪．孟桥．信号与线性系统．第四版．北京：高等教育出版社， 2004. [3] 刘永健．信号与线性系统．修订版．北京：人民邮电出版社，1994. [4] 奥本海姆 A V 等．信号与系统．第二版．刘树棠译．西安：西安交通大学出 版社，2002. [5] 刘培森．应用傅里叶变换．北京：北京理工大学出版社，1990. [6] 奥本海姆 A V 等．离散时间信号处理．黄建国，刘树棠译．北京：科学出版 社，1998. [7] 帕普利斯 A ．信号分析．毛培法译．北京：科学出版社，1981. [8] 丁玉美．高西全.数字信号处理．第二版.西安：西安电子科技大学出版社， 2001. [9] 乌镇杨.数字信号处理.北京：高等教育出版社.2004. [10] 斯坦利 W D.数字信号处理.常迥，译.北京：科学出版社，1979. [11] 王华奎.数字信号处理及应用.北京：高等教育出版社，2009. [12] 姚天任.江太辉.数字信号处理.第三版.华中科技大学出版社，2007. [13] 周耀华，汪凯仁.数字信号处理[M].上海：复旦大学出版社，1993. [14] 赵红怡.DSP 技术[M].北京：电子工业出版社，2008. [15] 刘松强.数字信号处理系统及其应用[M].北京：清华大学出版社，1996. [16] 陈后金主编.数字信号处理[M].北京：高等教育出版社，2004. [17] 胡广书.数字信号处理—理论、 算法与实现[M].北京： 清华大学出版社,1997. [18] 陈佩青.数字信号处理教程[M].2 版.北京.清华大学出版社，2001. [19] 门爱东，杨波，全子一.数字信号处理[M].北京：人民邮电出版社，2003. [20] 冷建华，李萍，王良红.数字信号处理[M].北京：国防工业出版社，2002. [21] 王风文，舒冬梅，赵宏才 .数字信号处理[M].北京：北京邮电大学出版社， 2005. [22] 邹彦，唐冬，宁志刚.DSP 原理与应用[M].北京：邮电大学出版社，2004. [23] 吴湘淇.信号、系统与信号处理（上下册）[M].北京：电子工业出版社，199

25

#include <iostream> #include <iomanip> #include "math.h" using namespace std; double pi = 3.1415926535897932; class complex { public: //无参构造函数 complex() { re=0; im=0; } //有参构造函数 complex(double real,double imag) { re=real; im=imag; } //加法 complex operator + (complex& c) { return complex( re + c.re , im + c.im ); } //减法 complex operator - (complex& c) { return complex( re - c.re , im - c.im ); 26

} //乘法 complex operator * (complex& c) { return complex( (re * c.re)-(im * c.im) , (re * c.im)+(im * c.re) ); } //除法 complex operator / (complex& c) { return complex( ( re*c.re + im*c.im )/( c.re*c.re + c.im*c.im ), ((im * c.re)-(re * c.im))/((c.re*c.re)+(c.im*c.im)) ); } //除 2 void half() { re=re/2; im=im/2; } //输出到 ostream void insert(ostream& out) { if(re>=0) cout<<" "; out<<re<< ((im>=0)?'+':'-') <<'j'<< ((im>=0)?im:(0-im)); } //显示 void show() { cout<<re<< ((im>=0)?'+':'-') <<'j'<< ((im>=0)?im:(0-im)) ; } //设值 void setValue(double real,double imag) { if(real>0.00000001||real< -0.00000001) re=real; 27

else re=0; if(imag>0.00000001||imag< -0.00000001) im=imag; else im=0; } friend complex c_pow(complex c,int y); private: double re,im; }; //流输出 ostream& operator << (ostream& out, complex c) { c.insert(out); return out; } //乘方 complex c_pow(complex c,int y) { complex returnValue = c ; for(int i = 1 ; i < y ; i++ ) { returnValue = returnValue * c ; } return returnValue; } //获得 Wn void getW(complex w[],int len) { for(int i=0 ; i<len ; i++ ) { w[i].setValue( cos(0 - pi*2*(i)/len) , sin(0 - pi*2*(i)/len) ); 28

} } //获得 mWn void getMW(complex w[],int len) { for(int i=0 ; i<len ; i++ ) { w[i].setValue( cos(pi*2*(i)/len) , sin(pi*2*(i)/len) ); } } //倒序重排 void resort(complex x[],int len) { //初始化 int half=len/2; int p1=0,p2; int k; complex temp; //重排 for( p2=0 ; p2+1<len ; p2++ ) { if(p2<p1 && p1+1<len) { temp=x[p1]; x[p1]=x[p2]; x[p2]=temp; } k=half; while(k<p1+1 && p1+1<len) { p1=p1-k; k=k/2; } p1=p1+k; 29

} } //基 2 时间 FFT void FFT2t(complex x[],complex y[],int len) { //初始化 int m=0,n,k,j,q; int tlen=len/2; while(tlen!=0) { m++; tlen=tlen/2; } resort(x,len); //倒序重排 for(j=0;j<len;j++) { y[j]=x[j]; } complex *xt = new complex[len]; complex *w = new complex[len]; getW(w,len); //开始迭代运算 q=len; for(int i=1;i<=m;i++) { q=q/2; n=len/q/2; for(j=0;j<len;j++) { xt[j]=y[j]; } if(i!=1) { for(j=0;j<q;j++) 30

{ for(k=0;k<n;k++) { xt[n+j*2*n+k]=xt[n+j*2*n+k]*w[k*q]; } } } for(j=0;j<n;j++) { for(k=0;k<q;k++) { y[j+k*2*n]=xt[j+k*2*n]+xt[j+n+k*2*n]; } } for(;j<2*n;j++) { for(k=0;k<q;k++) { y[j+k*2*n]=xt[j-n+k*2*n]-xt[j+k*2*n]; } } } delete[]xt,w; } //基 2 频率 FFT 和 IFFT void FFT(complex x[],complex y[],int len,bool IFFT = false) { //初始化 int m=0,n,k,j,q; int tlen=len/2; while(tlen!=0) { m++; tlen=tlen/2; 31

} for(j=0;j<len;j++) { y[j]=x[j]; } complex *xt = new complex[len]; complex *w = new complex[len]; if(IFFT==true) //若是 IFFT 则获取 Wn 序列 { getMW(w,len); } else //若是 FFT 则获取 mWn 序列 { getW(w,len); } //开始迭代运算 n=len; for(int i=1;i<=m;i++) { n=n/2; q=len/n/2; for(j=0;j<len;j++) { xt[j]=y[j]; } for(j=0;j<n;j++) { for(k=0;k<q;k++) { y[j+k*2*n]=xt[j+k*2*n]+xt[j+n+k*2*n]; } } for(;j<2*n;j++) { 32

for(k=0;k<q;k++) { y[j+k*2*n]=xt[j-n+k*2*n]-xt[j+k*2*n]; } } if(i!=m) { for(j=0;j<q;j++) { for(k=0;k<n;k++) { y[n+j*2*n+k]=y[n+j*2*n+k]*w[k*q]; } } } if(IFFT==true)//若是 IFFT 则除以 2 { for(j=0;j<len;j++) { y[j].half(); } } } resort(y,len); //倒序重排 delete[]xt,w; } //主函数 void main() { //初始化 complex q(0.9,0.3); complex x[256],y[256],z[256],w[256],v[256],w1; w1.setValue(1,0); getW(w,256); //获取 Wn 序列 33

x[0].setValue(1,0); for(int i=1;i<256;i++)//定原序列的值 { x[i]=x[i-1]*q; } //输出原序列 cout<<"原序列："<<endl; for(i=0;i<256;i+=2) { cout<<setiosflags(ios::fixed)<<setprecision(6)<<x[i]<<"\t"<<x[i+1]<<endl; } for(i=0;i<256;i++) //进行公式计算 { complex a=w1-(x[255]*q),b=w1-q*w[i]; v[i]=a/b; } //输出公式计算值 cout<<"公式计算值："<<endl; for(i=0;i<256;i+=2) { cout<<v[i]<<"\t"<<v[i+1]<<endl; } FFT2t(x,y,256); //进行基 2 时间变换 //输出基 2 时间变换后序列 cout<<"基 2 时间变换后序列："<<endl; for(i=0;i<256;i+=2) { cout<<y[i]<<"\t"<<y[i+1]<<endl; } FFT(y,z,256,1); //进行反变换 34

//输出反变换后的序列 cout<<"反变换后的序列："<<endl; for(i=0;i<256;i+=2) { cout<<z[i]<<"\t"<<z[i+1]<<endl; } }

35

36

37

38

39

40

41

42

43

45

FFT算法分析
FFT 算法分析 FFT 算法的基本原理是把长序列的 DFT...在基 2 算法中,每个蝶形运算单元都包括 1 次复数...C、D、E 处的信号标号 如下图所示: 串并转换...

dsp基于matlab的fft算法实现

8. 实验结果 10.源代码(C语言) #include "myapp...kfft(pr,pi,n,k,fr,fi,l,il):基 2 快速傅...加深了我 们对 DSP 特别是 FFT 算法的原理的理解...

(2)利用 FFT 计算逆 DFT 四、算法具体流程 五、算法的实际应用:计算大整数乘法 六、参考文献 一、引言 基于 FFT 的离散傅里叶变换(DFT)技术,是当今信息传输...
FFT实时谱分析系统的FPGA设计和实现
4 原位算法坐标旋转数字式计算机(CORDIC)算法 ...2.3 复乘运算模块 复乘运算是 FFT 处理器中两种...1M点FFT算法的FPGA设计与... 4页 1下载券 基于...