#include "myfft.h"
#include "data.h"
#define PI 3.1415926
//正弦函数（查表）
float mysin(float x){
if(x<0)
return -mysin(-x);
unsigned long deg=x/PI*180;
deg%=360;
return sin_tab[deg];
}
//余弦函数
float mycos(float x){
return mysin(PI/2-x);
}
//绝对值函数
float myfabs(float x)
{
return x>=0?x:-x;
}
//开方函数（牛顿迭代法）
float mysqrt(float x)
{
float val = x;//最终
float last;//保存上一个计算的值
do
{
last = val;
val =(val + x/val) / 2;
}
while(myfabs(val-last) > 0.01); //精度控制
return val;
}
//幂函数
int mypow(int i,int j){
int k=1;
for(;j>0;j--){
k*=i;
}
return k;
}
//对数函数
int mylog(int n,int a){//x=loga(n)
int x=0;
while(n>>=1){
x++;
}
return x;
}
//离散时间傅里叶变换
void dft(const float x[],float real[],float imag[],const int n,int isign){
if(isign!=1 && isign!=-1){//isign=1,正变换;isign=-1,逆变换
return;
}
int i,j;
float theta;
for(i=0;i<n;i++){
real[i]=0;
imag[i]=0;
for(j=0;j<n;j++){
theta=-2*PI/n*j*i*isign;//逆变换角度相反
real[i]+=x[j]*mycos(theta);
imag[i]+=x[j]*mysin(theta);
}
if(isign==-1){//逆变换需除以N
real[i]/=n;
imag[i]/=n;
}
}
}
//数位倒读
int rev(int i,int m){//i=0~2^m,m为二进制位数
int j=0;
while(m>0){
j+=(i & 0x01)*(0x01 << (m-1));//j+=(i%2)*mypow(2,m-1);
i>>=1;//i/=2
m-=1;
}
return j;
}
//快速傅里叶变换
void fft(const float real_in[],const float imag_in[],float real_out[],float imag_out[],const int n,int isign){
if(isign!=1 && isign!=-1){//isign=1,正变换;isign=-1,逆变换
return;
}
const int Lv=mylog(n,2);//蝶形级数
int L;//蝶形运算级数，用于循环
int N;//蝶形运算数据量，用于循环
int distance;//蝶形运算两节点间的距离，用于循环（distance=N/2）
int group;//蝶形运算的组数
float tmpr1,tmpi1,tmpr2,tmpi2;//临时变量
int i,j,k;
for(i=0;i<n;i++){//数位倒读
j=rev(i,Lv);
real_out[j]=real_in[i];
imag_out[j]=imag_in[i];
}
L=1;
distance=1;
N=2;
group=n>>1;
for(;L<=Lv;L++){//蝶形循环
for(i=0;i<group;i++){//每级蝶形各组循环
for(k=0;k<distance;k++){//每组蝶形运算
float theta=-2*PI*k/N*isign;//旋转因子,逆变换的角度与正变换相反
tmpr1=real_out[N*i+k];
tmpi1=imag_out[N*i+k];//X(k)
tmpr2=mycos(theta)*real_out[N*i+k+distance]-mysin(theta)*imag_out[N*i+k+distance];
tmpi2=mycos(theta)*imag_out[N*i+k+distance]+mysin(theta)*real_out[N*i+k+distance];//WN(k)*X(k+N/2)
real_out[N*i+k]=tmpr1+tmpr2;
imag_out[N*i+k]=tmpi1+tmpi2;//X(k)=X(k)+WN(K)*X(k+N/2)
real_out[N*i+k+distance]=tmpr1-tmpr2;
imag_out[N*i+k+distance]=tmpi1-tmpi2;//X(k+N/2)=X(k)-WN(K)*X(k+N/2)
if(isign==-1){//逆变换结果需除以N，即除以Lv次2
real_out[N*i+k]*=0.5;
imag_out[N*i+k]*=0.5;
real_out[N*i+k+distance]*=0.5;
imag_out[N*i+k+distance]*=0.5;
}
}
}
N<<=1;
distance<<=1;
group>>=1;
}
}
