#include <stdio.h>
#include <math.h>
#define PI 3.1415926
#define FFT_LENGTH 8
double input[FFT_LENGTH]={1,1,1,1,1,1,1,1};
struct complex1{ //定义一个复数结构体
double real; //实部
double image; //虚部
};
//将input的实数结果存放为复数
struct complex1 result_dat[8];
/*
虚数的乘法
*/
struct complex1 con_complex(struct complex1 a,struct complex1 b){
struct complex1 temp;
temp.real=(a.real*b.real)-(a.image*b.image);
temp.image=(a.image*b.real)+(a.real*b.image);
return temp;
}
/*
简单的a的b次方
*/
int mypow(int a,int b){
int i,sum=a;
if(b==0)return 1;
for(i=1;i<b;i++){
sum*=a;
}
return sum;
}
/*
简单的求以2为底的正整数
*/
int log2(int n){
unsigned i=1;
int sum=1;
for(i;;i++){
sum*=2;
if(sum>=n)break;
}
return i;
}
/*
简单的交换数据的函数
*/
void swap(struct complex1 *a,struct complex1 *b){
struct complex1 temp;
temp=*a;
*a=*b;
*b=temp;
}
/*
dat为输入数据的数组
N为抽样次数 也代表周期 必须是2^N次方
*/
void fft(struct complex1 dat[],unsigned char N){
/*最终 dat_buf计算出 当前蝶形运算奇数项与W 乘积
dat_org存放上一个偶数项的值
*/
struct complex1 dat_buf,dat_org;
/* L为几级蝶形运算 也代表了2进制的位数
n为当前级蝶形的需要次数 n最初为N/2 每级蝶形运算后都要/2
i j为倒位时要用到的自增符号 同时 i也用到了L碟级数 j是计算当前碟级的计算次数
re_i i_copy均是倒位时用到的变量
k为当前碟级 cos(2*pi/N*k)的 k 也是e^(-j2*pi/N)*k 的 k
*/
unsigned char L,i,j,re_i=0,i_copy=0,k=0,fft_flag=1;
//经过观察，发现每级蝶形运算需要N/2次运算，共运算N/2*log2N 次
unsigned char fft_counter=0;
//在此要进行补2 N必须是2^n 在此略
//蝶形级数 (L级)
L=log2(N);
//计算每级蝶形计算的次数(这里只是一个初始值) 之后每次要/2
//n=N/2;
//对dat的顺序进行倒位
for(i=1;i<N/2;i++){
i_copy=i;
re_i=0;
for(j=L-1;j>0;j--){
//判断i的副本最低位的数字 并且移动到最高位 次高位 ..
//re_i为交换的数 每次它的数字是不能移动的 并且循环之后要清0
re_i|=((i_copy&0x01)<<j);
i_copy>>=1;
}
swap(&dat[i],&dat[re_i]);
}
//进行fft计算
for(i=0;i<L;i++){
fft_flag=1;
fft_counter=0;
for(j=0;j<N;j++){
if(fft_counter==mypow(2,i)){ //控制隔几次，运算几次,
fft_flag=0;
}else if(fft_counter==0){ //休止结束，继续运算
fft_flag=1;
}
//当不判断这个语句的时候 fft_flag保持 这样就可以持续运算了
if(fft_flag){
dat_buf.real=cos((2*PI*k)/(N/mypow(2,L-i-1)));
dat_buf.image=-sin((2*PI*k)/(N/mypow(2,L-i-1)));
dat_buf=con_complex(dat[j+mypow(2,i)],dat_buf);
//计算 当前蝶形运算奇数项与W 乘积
dat_org.real=dat[j].real;
dat_org.image=dat[j].image; //暂存
dat[j].real=dat_org.real+dat_buf.real;
dat[j].image=dat_org.image+dat_buf.image;
//实部加实部 虚部加虚部
dat[j+mypow(2,i)].real=dat_org.real-dat_buf.real;
dat[j+mypow(2,i)].image=dat_org.image-dat_buf.image;
//实部减实部 虚部减虚部
k++;
fft_counter++;
}else{
fft_counter--; //运算几次，就休止几次
k=0;
}
}
}
}
void main(){
int i;
//先将输入信号转换成复数
for(i=0;i<FFT_LENGTH;i++){
result_dat[i].image=0;
//输入信号是二维的，暂时不存在复数
result_dat[i].real=input[i];
//result_dat[i].real=10;
//输入信号都为实数
}
fft(result_dat,FFT_LENGTH);
for(i=0;i<FFT_LENGTH;i++){
input[i]=sqrt(result_dat[i].real*result_dat[i].real+result_dat[i].image*result_dat[i].image);
//取模
printf("%lf\n",input[i]);
}
while(1);
}
| 
楼主
标题置顶
标题高亮
