快速傅里叶变换FFT的C语言实现及应用(2)

Wr 的确定N 蝶形运算两节点的第一个节点为k值，表示成L位二进制数，左移L – m 输入序列x(n) : N个存储单元 r的二进制数。 位，把右边空出的位置补零，结果为 系数 ：N / 2个存储单元

4）存储单元

快速傅立叶变换的C语言实现方法

rN 有了傅立叶变换，我们可以从信号的频域特征去分析信号。尤其在无线通信系统 中，傅里叶变换的重要性就更加明显了，无论是设计者还是测试工程师，在工作中都会和傅立叶变换打交道。

W我们要衡量一个处理器有没有足够的能力来运行FFT算法，根据以上的简单介绍可以得出以下两点：

1. 处理器要在一个指令周期能完成乘和累加的工作，因为复数运算要多次查表相乘

才能实现。

2. 间接寻址，可以实现增/减1个变址量，方便各种查表方法。FFT要对原始序列进

行反序排列，处理器要有反序间接寻址的能力。

下面为一份FFT（快速傅立叶变换）的源码（基于C） /************FFT***********/ #include #include #include

#define N 1000 typedef struct {

double real;/*实部*/ double img;/*虚部*/ }complex;

void fft(); /*快速傅里叶变换*/ void ifft(); /*快速傅里叶逆变换*/ 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 x[N], *W;/*输出序列的值*/

int size_x=0;/*输入序列的长度，只限2的N次方*/ double PI;

int main() {

int i,method;

system(\

PI=atan(1)*4;/*pi等于4乘以1.0的正切值*/

printf(\ /*输入序列的长度*/ scanf(\

printf(\ /*输入序列对应的值*/ for(i=0;i

scanf(\ initW();

/*选择FFT或逆FFT运算*/

printf(\ scanf(\ if(method==0) fft(); else ifft();

output(); system(\ return 0; }

/*进行基-2 FFT运算*/ void fft() {

int i=0,j=0,k=0,l=0; complex up,down,product; change();

for(i=0;i< log(size_x)/log(2) ;i++) /*一级蝶形运算*/ { l=1<

for(j=0;j

for(k=0;k

mul(x[j+k+l],W[size_x*k/2/l],&product); add(x[j+k],product,&up); sub(x[j+k],product,&down); x[j+k]=up; x[j+k+l]=down; } } } }

void ifft() {

int i=0,j=0,k=0,l=size_x; complex up,down;

for(i=0;i< (int)( log(size_x)/log(2) );i++) /*一级蝶形运算*/ { l/=2;

for(j=0;j

for(k=0;k

add(x[j+k],x[j+k+l],&up); up.real/=2;up.img/=2;

sub(x[j+k],x[j+k+l],&down); down.real/=2;down.img/=2;

divi(down,W[size_x*k/2/l],&down); x[j+k]=up; x[j+k+l]=down; } } }

change(); }

/*初始化变化核*/ void initW() {

int i;

W=(complex *)malloc(sizeof(complex) * size_x); for(i=0;i

W[i].real=cos(2*PI/size_x*i); W[i].img=-1*sin(2*PI/size_x*i); } }

/*变址计算，将x(n)码位倒置*/ void change() {

complex temp;

unsigned short i=0,j=0,k=0; double t;

for(i=0;i

k=i;j=0;

t=(log(size_x)/log(2)); while( (t--)>0 )

{

j=j<<1;

j|=(k & 1); k=k>>1; }

if(j>i) {

temp=x[i]; x[i]=x[j]; x[j]=temp; } } }

void output() /*输出结果*/ {

int i;

printf(\ for(i=0;i

printf(\ if(x[i].img>=0.0001)

printf(\ else if(fabs(x[i].img)<0.0001) printf(\ else

printf(\ } }

void add(complex a,complex b,complex *c) {

c->real=a.real+b.real; c->img=a.img+b.img; }

void mul(complex a,complex b,complex *c) {