/*FFT.C*/
#include "Stdio.h"
#include "Conio.h"
#include "Math.h"
#define pi 3.141592f
/* 参量与变量 */
FILE *fp;
int sample_rate; /* 采样速率 */
double sample_x[256]; /* 采样信号*/
double signal_fft[256]; /* FFT结果 */
short int bit_reversing_flag[256];
double x_0_r[256];
double x_0_i[256];
double x_1_r[256];
double x_1_i[256];
int step; /* 数据传输时用的步进值 */
int fft_length; /* fft计算终点数值 */
int fft_times; /* fft计算次数 */
/* 函数 */
void Get_Data( void ); /* 获取数据 */
void Init( void ); /* 清空*/
void FFT( void ); /* 快速傅里叶变换FFT */
void FFT_Reverse_Data( void ); /* 倒序 */
void FFT_Algorithm( void ); /* 蝶形 */
int Reverse_Order( int num, int kind ); /* 求出一个下标的倒序下标 */
void Complex_Add( double *a, double *b, double c, double d ); /* 两个复数相加 */
void Complex_Rotate( double *a, double *b, double N, double n ); /* 一复数乘旋转因子 */
void Show_FFT_Result( void ); /* 显示快速傅里叶变换结果 */
void main(void)
{
Init(); /* 清空数据 */
Get_Data(); /* 获取数据 */
FFT_Reverse_Data(); /* 倒序 */
FFT_Algorithm(); /* 蝶形 */
Show_FFT_Result();
}
void Get_Data( void ) /* 获取数据 */
{
int i, j;
for( i = 0, j = 0; i < 256; i+=step, j++ ) /* 数据传输 */
x_0_r[j] = 0.0,
x_0_r[j] = sample_x[i];
x_0_r[0] = -20;
x_0_r[1] = 21;
x_0_r[2] = 10;
x_0_r[3] = 10;
x_0_r[4] = 9;
x_0_r[5] = -4;
x_0_r[6] = 32;
x_0_r[7] = 29;
x_0_r[8] = -27;
x_0_r[9] = 7;
x_0_r[10] = 43;
x_0_r[11] = -13;
x_0_r[12] = -17;
x_0_r[13] = 18;
x_0_r[14] = 2;
x_0_r[15] = 0;
}
void Init( void ) /* 清空 */
{
int i;
sample_rate = 0;
for( i = 0; i<256; i++ ) {
sample_x[i] = 0.0f;
signal_fft[i] = 0.0f;
bit_reversing_flag[i] = 0;
x_0_r[i] = 0.0f;
x_0_i[i] = 0.0f;
x_1_r[i] = 0.0f;
x_1_i[i] = 0.0f;
}
sample_rate = 0; /* 设置采样速率 */
if( sample_rate == 0 ) { step =16; fft_length = 16; fft_times = 4; } /* 16 */
if( sample_rate == 1 ) { step = 8; fft_length = 32; fft_times = 5; } /* 32 */
if( sample_rate == 2 ) { step = 4; fft_length = 64; fft_times = 6; } /* 64 */
if( sample_rate == 3 ) { step = 2; fft_length = 128; fft_times = 7; } /* 128 */
if( sample_rate == 4 ) { step = 1; fft_length = 256; fft_times = 8; } /* 256 */
}
int Reverse_Order( int num, int kind ) /* 求出一个下标的倒序下标 */
{
int ref_0[13], ref_1[13];
int ref_2[13] = { 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048, 4096 };
int i;
int change_data;
int end_order;
if( kind == 16 ) end_order = 4;
if( kind == 32 ) end_order = 5;
if( kind == 64 ) end_order = 6;
if( kind == 128 ) end_order = 7;
if( kind == 256 ) end_order = 8;
for( change_data = num, i = 0; i < end_order; i++ ) {
ref_0[i] = change_data % 2;
change_data = change_data / 2; }
for( i = 0 ; i < end_order; i++ )
ref_1[i] = ref_0[ end_order - i - 1 ];
for( change_data = 0, i = 0; i < end_order; i++ ) {
ref_1[i] = ref_1[i] * ref_2[i];
change_data += ref_1[i]; }
return change_data;
}
void Complex_Add( double *a, double *b, double c, double d ) /* 两个复数相加 */
{
*a = *a + c;
*b = *b + d;
}
void Complex_Multiply( double *a, double *b, double c, double d ) /* 两个复数相乘 */
{
double x_r = *a, x_i = *b;
double y_r = c, y_i = d;
*a = x_r * y_r - x_i * y_i;
*b = x_r * y_i + y_r * x_i;
}
void Complex_Rotate( double *a, double *b, double N, double n ) /* 一复数乘旋转因子 */
{
double x_r = *a, x_i = *b;
double y_r, y_i;
y_r = (double)cos( -2.0f * pi * n / N );
y_i = (double)sin( -2.0f * pi * n / N );
Complex_Multiply( &x_r, &x_i, y_r, y_i );
*a = x_r;
*b = x_i;
}
void FFT_Reverse_Data( void ) /* 倒序 */
{
int i, j;
double ref_data;
for( i = 0; i < fft_length; i++ ) {
if( bit_reversing_flag[i] == 0 ) {
j = Reverse_Order( i, fft_length );
ref_data = x_0_r[j];
x_0_r[j] = x_0_r[i];
x_0_r[i] = ref_data;
bit_reversing_flag[i] = 1;
bit_reversing_flag[j] = 1;
}
}
}
void FFT_Algorithm( void )
{
int group_list[14] = { 1, 2, 4, 8, 16, 32, 64, 128, /* 用于蝶形运算 */
256, 512, 1024, 2048, 4096, 8192 };
double twiddle_factor_list[13] = { 2.0f, 4.0f, 8.0f, 16.0f, /* 用于旋转因子系数 */
32.0f, 64.0f, 128.0f, 256.0f, 512.0f,
1024.0f, 2048.0f, 4096.0f, 8192.0f };
int i, j, k, l, m, n;
double ref_order;
double ref_f;
double x_r, x_i, y_r, y_i;
for( i = 0; i < fft_times; i++ ) { /* 列的计算总次数 */
m = group_list[i];
for( j = 0; j < group_list[ fft_times - i - 1 ]; j++ ) { /* 一列内的组数 */
for( k = 0; k < m; k++ ) { /* 一组内的一半成员的个数 */
l = j * m * 2;
x_r = x_0_r[ l + k + m ]; /* 乘旋转因子 */
x_i = x_0_i[ l + k + m ];
y_r = twiddle_factor_list[ i ];
y_i = (double)k;
Complex_Rotate( &x_r, &x_i, y_r, y_i );
y_r = x_0_r[ l + k ]; /* 相加 */
y_i = x_0_i[ l + k ];
Complex_Add( &x_r, &x_i, y_r, y_i );
x_1_r[ l + k ] = x_r;
x_1_i[ l + k ] = x_i;
x_r = x_0_r[ l + k + m ]; /* 乘旋转因子 */
x_i = x_0_i[ l + k + m ];
y_r = twiddle_factor_list[ i ];
y_i = (double)k;
Complex_Rotate( &x_r, &x_i, y_r, y_i );
x_r = x_r * ( -1.0f ); /* 相减 */
x_i = x_i * ( -1.0f );
y_r = x_0_r[ l + k ];
y_i = x_0_i[ l + k ];
Complex_Add( &x_r, &x_i, y_r, y_i );
x_1_r[ l + k + m ] = x_r;
x_1_i[ l + k + m ] = x_i;
}
}
for( n = 0; n < fft_length; n++ ) { /* 交换数据以循环计算 */
x_0_r[n] = x_1_r[n];
x_0_i[n] = x_1_i[n];
}
}
}
void Show_FFT_Result( void ) /* 显示快速傅里叶变换结果 */
{
int n;
if((fp=fopen("Result.txt","w"))==NULL) {
printf("Can't open file!");
exit(0);
}
for( n = 0; n < fft_length; n++ ) {
fprintf( fp, "X%2d : %f, %f ", n, sqrt( x_0_r[n] * x_0_r[n] + x_0_i[n] * x_0_i[n] ), x_1_r[n] );
if( x_1_i[n] >= 0.0 ) fprintf( fp, "+%fi\n", x_1_i[n] );
else fprintf( fp, "%fi\n", x_1_i[n] );
}
} |
楼主
标题置顶
标题高亮
