/* 新手上路还望见谅。 *
- #include <stdio.h>
- #include <math.h>
- #include <stdlib.h>
- #define N 8 //64 输入样本总数
- #define M 3 //DFT运算层数 //2^m=N
- #define PI 3.1415926
- float twiddle[N/2] = {1.0, 0.707, 0.0, -0.707};
- float x_r[N] = {1, 1, 1, 1, 0, 0, 0, 0}; //输入数据，此处设为8个
- float x_i[N]; //N=8
- /**
- * 初始化输出虚部
- */
- static void fft_init( void )
- {
- int i;
- for(i=0; i<N; i++) x_i[i] = 0.0;
- }
-
- /**
- * 反转算法.将时域信号重新排序.
- * 这个算法有改进的空间
- */
- static void bitrev( void )
- {
- int p=1, q, i;
- int bit_rev[ N ];
- float xx_r[ N ];
-
- bit_rev[ 0 ] = 0;
- while( p < N )
- {
- for(q=0; q<p; q++)
- {
- bit_rev[ q ] = bit_rev[ q ] * 2;
- bit_rev[ q + p ] = bit_rev[ q ] + 1;
- }
- p *= 2;
- }
-
- for(i=0; i<N; i++) xx_r[ i ] = x_r[ i ];
-
- for(i=0; i<N; i++) x_r[i] = xx_r[ bit_rev[i] ];
- }
- void fft( void )
- { fp = fopen("log2.txt", "a+");//此处
- int cur_layer, gr_num, i, k, p; //cur_layer代表正要计算的当前层，gr_num代表当前层的颗粒数
- float tmp_real, tmp_imag, temp; // 临时变量, 记录实部
- float tw1, tw2;// 旋转因子,tw1为旋转因子的实部cos部分, tw2为旋转因子的虚部sin部分.
-
- int step; // 步进
- int sample_num; // 颗粒的样本总数(各层不同, 因为各层颗粒的输入不同)
-
- /* 对层循环 */
- for(cur_layer=1; cur_layer<=M; cur_layer++)
- {
- /* 求当前层拥有多少个颗粒(gr_num) */
- gr_num = 1;
- i = M - cur_layer;
- while(i > 0)
- {
- i--;
- gr_num *= 2;
- }
-
- /* 每个颗粒的输入样本数N' */
- sample_num = (int)pow(2, cur_layer);
- /* 步进. 步进是N'/2 */
- step = sample_num/2;
-
- /* */
- k = 0;
-
- /* 对颗粒进行循环 */
- for(i=0; i<gr_num; i++)
- {
- /*
- * 对样本点进行循环, 注意上限和步进
- */
- for(p=0; p<sample_num/2; p++)
- {
- // 旋转因子, 需要优化...
- tw1 = cos(2*PI*p/pow(2, cur_layer));
- tw2 = -sin(2*PI*p/pow(2, cur_layer));
-
- tmp_real = x_r[k+p];
- tmp_imag = x_i[k+p];
- temp = x_r[k+p+step];
-
- /* 蝶形算法 */
- x_r[k+p] = tmp_real + ( tw1*x_r[k+p+step] - tw2*x_i[k+p+step] );
- x_i[k+p] = tmp_imag + ( tw2*x_r[k+p+step] + tw1*x_i[k+p+step] );
- /* X[k] = A(k)+WB(k)
- * X[k+N/2] = A(k)-WB(k) 的性质可以优化这里*/
- /*旋转因子, 需要优化...
- tw1 = cos(2*PI*(p+step)/pow(2, cur_layer));
- tw2 = -sin(2*PI*(p+step)/pow(2, cur_layer));
- x_r[k+p+step] = tmp_real + ( tw1*temp - tw2*x_i[k+p+step] );
- x_i[k+p+step] = tmp_imag + ( tw2*temp + tw1*x_i[k+p+step] );*/
- x_r[k+p+step] = tmp_real - ( tw1* temp - tw2*x_i[k+p+step] );
- x_i[k+p+step] = tmp_imag - ( tw2* temp + tw1*x_i[k+p+step] );
-
- printf("k=%d, x_r[k]=%f, x_i[k]=%f\n", k+p, x_r[k+p], x_i[k+p]);
- printf("k=%d, x_r[k]=%f, x_i[k]=%f\n", k+p+step, x_r[k+p+step], x_i[k+p+step]);
- }
- /* 开跳!:) */
- k += 2*step;
- }
- }
- }
- void display( void )
- {
- printf("\n\n");
- int i;
- for(i=0; i<N; i++)
- printf("%f\t%f\n", x_r[i], x_i[i]);
- }
- int main( void )
- {
- fft_init( ); //初始化
- bitrev( ); //将输入直接按FFT计算要求排序，如8点FFT计算，排序为x[0]、x[4]、x[2]、x[6]、x[1]、x[5]、x[3]、x[7]
- fft( ); //进行FFT计算
- display( ); //显示计算结果
-
- system( "pause" );
- return 1;
- }
|
