当前位置:首页 >> 工学 >>

FFT的C语言编程


FFT 的 C 语言编程
1.程序:
#include <stdio.h> #include <math.h> #define re 0 /* re=0,用 re 表示实部*/ #define im 1 /* im=1,用 im 表示虚部*/ main() { float x[128][2],w[2],temp[2]; /* x[128][2]: 复数变量; x[i][re]: 第 i 个复数变量的实部; x[i][im]: 第 i 个复数变量的虚部; w[2]: 存储旋转因子 WN P,w[re]、w[im]分别代表旋转因子的实部和虚部; temp[2]: 蝶形计算中的临时变量,temp[re]、temp[im]分别代表其实部 和虚部; */ float arg,wreal,wimag; /* arg 存储旋转因子指数 p(数值上相差-2 π/N)。wreal 存储 cos(arg), wimag 存储 -sin(arg) */ float tem,tr,ti; int L,M,B,j,i,k, N,N2; char c='i'; scanf("%d %d",&N,&M); /* 输入复数信号长度 N,蝶形运算级数 M */ N2=N>>1; for(j=0;j<N;j++) /* 输入复数信号实部和虚部*/ { scanf("%f %f",&x[j][re],&x[j][im]);

} printf("\n"); /*输入倒序*/ for(j=0,i=1;i<N-1;i++) { k=N2; while(k<=j) { j=j-k; k=k>>1; } j=j+k; if(i<j) { tr=x[j][re]; ti=x[j][im]; x[j][re]=x[i][re]; x[j][im]=x[i][im]; x[i][re]=tr; x[i][im]=ti; } } /*FFT 三重循环模块*/ for(L=1; L<=M; L++) /* 逐级进行计算共 M 级*/ { B=1<<L-1; /* 第 L 级共有 B=2L-1 个不同的旋转因子*/ arg=-acos(-1)/B; /* 旋转因子初始化注释见结尾处*/ w[re]=cos(arg); w[im]=-sin(arg); for(j=0; j<B; j++) /* j 代表第 L 级不同旋转因子的个数*/ { /* 旋转因子*/ arg=acos(-1)/B; /* arg=π/Β */ wreal=cos(arg); wimag= -sin(arg); tem=w[re]*wreal-w[im]*wimag; w[im]=w[re]*wimag+w[im]*wreal;

w[re]=tem; for(k=j; k<N; k+=2*B) /*第 L 级具有相同旋转因子蝶形计算,每个蝶形相距 2L=2B 个点*/ { temp[re]=x[k+B][re]*w[re]-x[k+B][im]*w[im]; temp[im]=x[k+B][im]*w[re]+x[k+B][re]*w[im]; x[k+B][re]=x[k][re]-temp[re]; x[k+B][im]=x[k][im]-temp[im]; x[k][re]=x[k][re]+temp[re]; x[k][im]=x[k][im]+temp[im]; /* 编写蝶形运算程序*/ /* 第 L 级每个蝶形计算的输入节点距离为 B */ /* 蝶形运算 ? + + ] ] B B P P N ? ? + [ [ i i N */ ? = = ] ] x x + [ [ W W ] [ x[i] i x B i x i x

/* 利用临时存储变量 temp[2]计算 WN Px[i+B] */ /* 复数运算:(a+bj)(c+dj)=(ac-bd)+(bc+ad)j */ /* temp[re]= ac-bd,temp[im]= bc+ad */ } } } for(j=0;j<N;j++) /*输出*/ { printf("%f %c%f\n",x[j][re],c,x[j][im]); }

} /* 计算旋转因子说明 ) ( ^ 2 * L M j P ? = 、1 2 ? = L B 、M N 2 = 有θ θ θ π π sin cos ) / ( 2 i e e e W i B j i P N i P N ? = = = = ? ? ? 其中 B j / π θ= 下一个旋转因子旋转角度为 arg (arg) * i i e e ? ? θ cos(arg)) * sin sin(arg) * (cos sin(arg)) * sin cos(arg) * (cos θ θ θ θ + ? ? = i */

2. 3.

8 点复数信号的离散傅里叶变换: 16、32 点复数信号的 FFT 频谱图

1)16 点: function kuangxin1 x=[ 0.000000 0.000000 0.000000 0.000000 0.000000 -4.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000

0.000000 4.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.828400 -2.828400 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 2.828400 2.828400 0.000000 0.000000]; N=16; for i=1:16 y(i)=sqrt(x(i,1)^2+x(i,2)^2); end i=0:15; stem(i,y)

4 3.5 3 2.5 2 1.5 1 0.5 0

0

5

10

15

2)32 点: function kuangxin2 x=[0.000000 0.000000 0.000000 0.000000 0.000000 -4.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 4.000000 0.000000 0.000000 0.000000 0.000000

0.000000 0.000000 2.828400 -2.828400 0.000000 0.000000 8.000000 0.000000 0.000000 0.000000 2.828400 2.828400 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 1.530800 -3.695600 0.000000 0.000000 5.656800 0.000000 0.000000 0.000000 1.530800 3.695600 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 3.695600 -1.530800 0.000000 0.000000 5.656800 0.000000 0.000000 0.000000 3.695600 1.530800

0.000000 0.000000];

for i=1:32 y(i)=sqrt(x(i,1)^2+x(i,2)^2); end i=0:31; stem(i,y)

8 7 6 5 4 3 2 1 0

0

5

10

15

20

25

30

35


赞助商链接
相关文章:
更多相关标签: