#include <stdio.h>
#include <math.h>
#include "fft_lib.h"
//复数结构体定义
struct complex{
double real;
double imag;
};
//复数加法
complex complex_add(complex a1,complex a2)
{
complex out;
out.real=a1.real+a2.real;
out.imag=a1.imag+a2.imag;
return out;
}
//复数减法
complex complex_sub(complex a1,complex a2)
{
complex out;
out.real=a1.real-a2.real;
out.imag=a1.imag-a2.imag;
return out;
}
//复数乘法
complex complex_multiply(complex a1, complex a2)
{
complex out;
out.real = a1.real * a2.real - a1.imag * a2.imag;
out.imag = a1.real * a2.imag + a1.imag * a2.real;
return out;
}
//欧拉公式 :fft里面的旋转因子
complex euler_formula(double P,int N)
{
double x;
complex out;
x=2*PI/N*P;
out.real=cos(x);
out.imag=-sin(x);
return out;
}
/*
说明: m!=0:返回数值为该数值是2的m次方
m=0:输入数值有问题
*/
int charge_number(int N)
{
int m=0;
for(;N%2==0;N=N/2){
if(N%2==0)
m++;
}
if(N!=1)
m=0;
return m;
}
//倒序
void reverse_order(complex A[],int N)
{
int i,j; //正序和倒序位置
int lh,n1;
complex index;
int k;
lh=N/2;
n1=N-2;
j=lh;
//重新排序部分
for(i=1;i<n1;i++){
if(i<j){
index=A[j];
A[j]=A[i];
A[i]=index;
}
k=lh;
while(j>=k){
j=j-k;
k=k/2;
}
j=j+k;
}
}
//蝶形运算
void count_fft(complex A[],int M,int N)
{
int L; //当前在那一级蝶形运算
int B; //跨度
int J,K; //旋转因子以及蝶形运算
int P; //旋转因子系数
complex W; //旋转因子
complex index; //计算过程中的中间值
complex index1; //用来暂存A[K]某一时刻的值
for(L=1;L<=M;L++){
B=(int)pow(2,L-1);
for(J=0;J<=B-1;J++){
P=(int)(J*pow(2,M-L));
for(K=J;K<=N-1;K=K+2*B){
W=euler_formula(P,N);
index=complex_multiply(A[K+B],W);
index1=A[K]; //这个中间变量十分重要,如果不使用中间变量,下一行代码执行后,A[K]值就已经改变 后半段一定就是错的了
A[K]=complex_add(index1,index);
A[K+B]=complex_sub(index1,index);
}
}
}
}
void fft(complex A[],int M,int N)
{
int i=0;
reverse_order(A,N); //倒序
count_fft(A,M,N); //蝶形运算
printf("FFT的结果为:\n");
for(i=0;i<N;i++)
printf("第%d组的值为: %-f\t+ %fj \n",i+1,A[i].real,A[i].imag);
}