信号处理 FFT算法

实验2 基2时域抽选的FFT 程序设计与调试

一、实验目的

掌握信号处理，尤其是数字信号处理的基本原理和方法。要求能通过实验熟练掌握基2时域抽选的快速傅立叶变换算法（FFT ）的基本原理，了解二维及多维快速傅立叶变换算法。

二、实验原理

1．复数类型

对于FFT 算法涉及的复数运算，使用自定义的COMPLEX 来定义复数类型，其使用方法与常规类型（如int,float,double ）相似。 typedef struct { float real, imag; } COMPLEX; 2．FFT 基本原理

FFT 改进了DFT 的算法，减少了运算量，主要是利用了旋转因子W 的两个性质：

（a ）W 的周期性：W = W （b) W 的对称性：W =-W

FFT 把N 点DFT 运算分解为两组N/2点的DFT 运算，然后求和：

)()()(21k X W k X k X k N +=

1,,1,0 ),()()2

(2

21-=-=+

N k

N k k X W k X N k X 其中，

∑∑∑∑-=-=-=-=+==

=

=

1

1

2

21

1

112

2

2

2

2

2

2

2

)12()()()2()()(N N N

N N N N N r rk

r rk r rk

r rk

W r x W r x

k X W

r x W

r x k X

在计算X 1(k)与X 2(k)时，仍利用上述公式，把它们看成是新的X(k)。如此递归下去，便是FFT 算法。 3．蝶形运算

从基2时域抽选FFT 运算流图可知：

① 蝶形两节点的距离为2m-1，其中，m 表示第m 列，且m =1,… ,L 。 例如N=8=23， 第一级(列)距离为21-1=1， 第二级(列)距离为22-1=2， 第三级(列)距离为23-1=4。

② 考虑蝶形运算两节点的距离为2m-1，蝶形运算可表为： X m (k)=X m-1(k)+X m-1(k+2m-1) W N r

X m (k+2m-1)= X m-1(k)-X m-1(k+2m-1) W N r

由于N 为已知，所以将r 的值确定即可确定W N r 。为此，令k=(n 2n 1n 0)2 ,再将k 左移(L-m)位，右边位置补零，就可得到(r)2 的值，即(r)2 =(k)22L-m 。

例如 N=8=23

(1) k=2 , m=3 的r 值,∵k=2=(010)2 左移L-m=3-3=0 ,∴ r=(010)2=2; (2) k=3 ,m=3 的r 值;k= 3=(011)2,左移0位,r=3；

(3) k=5 ,m=2的值;k=5=(101)2 左移L-m=1位 r=(010)2 =2。 4．存储单元

存输入序列 x(n)，n=0,1, ,N-1，计N 个单元；

存放系数W N r ， r=0,1, ,（N/2）-1；需N/2个存储单元； 共计(N+N/2)个存储单元。 5．程序框图

开 始

送入x (n )，M

N ＝2 M

倒 序

L ＝1 , M

Ｊ＝0 , B £ 1

P ＝2 M £L J

k ＝ J , N £1 , 2L

p

N

p

N

W B k X k X B k X W B k X k X k X )()()()()()(+-⇐+++⇐输出

结 束

B 2 L £

1

三、实验内容

1．编写一个实现FFT 算法的函数。

void fft(COMPLEX A[ ] , int m) 其中，COMPLEX 是复数类型的名称，A 是要变换的数据，m 是FFT 点数N 的以2为底的幂次(即N=2m )。

2．编写测试用的主程序对fft 函数进行测试，并给出测试结果。

四、实验代码

#include #include

#include

#define N 1000

#define PI atan(1)*4

typedef struct

{

double real;

double imag;

}complex;

void fft(complex A[],int m); //快速傅里叶变换

void initW();

void change();

void add(complex,complex,complex *); //复数加法void mul(complex,complex,complex *); //复数乘法void sub(complex,complex,complex *); //复数减法void divi(complex,complex,complex *); //复数除法void output();

complex A[N], *W; //输出序列的值

complex up,down,product;

int size=0; //序列的长度

int m;

void main()

{

fft(A,m);

output();

}

//进行基-2 FFT运算

void fft(complex x[],int m)

{ int i=0,j=0,k=0,l=0;

printf("Please input m:");

scanf("%d",&m);

size=pow(2,m);

printf("Please input %d datas in A[%d]:\n",size,size);

for(i=0;i

{

printf("A[%d]:",i);

scanf("%lf %lf",&A[i].real,&A[i].imag);

}

initW();

change();