/////////////////////////////////////////////////////
//FFT算法
/////////////////////////////////////////////////////
typedef struct Complex//复数结构体
{
double Re;
double Im;
}Complex;
Complex Mul(Complex a,Complex b)//复数相乘
{
Complex c;
c.Re=a.Re*b.Re-a.Im*b.Im;
c.Im=a.Re*b.Im+a.Im*b.Re;
return c;
}
Complex Plus(Complex a,Complex b)//复数相加
{
Complex c;
c.Re=a.Re+b.Re;
c.Im=a.Im+b.Im;
return c;
}
Complex Sub(Complex a,Complex b)//复数相减
{
Complex c;
c.Re=a.Re-b.Re;
c.Im=a.Im-b.Im;
return c;
}
/******************************************************************************************
函数说明：倒序运算
参数说明：N点数，等于2的M次方，*x序列的首地址
返回值：无
******************************************************************************************/
/*void Invert_order
(uint N,Complex *x,uchar M)
{
uint power,i,j,m,tempp;//m是i的倒位序数
Complex temp;//是临时变量
for(i=0;i<N;i++)
{
power=1;
m=0;
for(j=0;j<M;j++)
{
tempp=i;
tempp>>=(M-1-j);//C语言 n>>=1 中的>>=意思是先将变量n的各个二进制位顺序右移1位，最高位补二进制0，然后将这个结果再复制给n。
if(tempp&0x01)//判断tempp最后一位是不是1
{
m=m+power;
}
power=power*2;
}
if(i<m)
{
temp.Im=x[i].Im;
temp.Re=x[i].Re;
x[i].Re=x[m].Re;
x[i].Im=x[m].Im;
x[m].Im=temp.Im;
x[m].Re=temp.Re;
}
for(int n=0;n<N;n++)
{
if(x[n].Re<-10000000)
{
x[n].Re=0;
}
if(x[n].Im<-10000000)
{
x[n].Im=0;
}
}
}*/
void Invert_order(uint N,Complex *x)
{
int LH,j,N1,i,k,n;
Complex T;
LH=N/2;
j=LH;
N1=N-2;
for(i=1;i<=N1;i++)
{
if(i<j)
{
T.Im=x[i].Im;
T.Re=x[i].Re;
x[i].Re=x[j].Re;
x[i].Im=x[j].Im;
x[j].Im=T.Im;
x[j].Re=T.Re;
}
k=LH;
while(j>=k)
{
j=j-k;
k=k/2;
}
j=j+k;
}
for(n=0;n<N;n++)
{
if(x[n].Re<-10000000)
{
x[n].Re=0;
}
if(x[n].Im<-10000000)
{
x[n].Im=0;
}
}
}
/******************************************************************************************
Invert_order函数结束：
******************************************************************************************/
/******************************************************************************************
函数说明：指数运算
参数说明：BaseNumber底数,IndexNumber指数
返回值：无符号整数
******************************************************************************************/
uint Pow(uchar BaseNumber,uint IndexNumber)
{
uint m=1;//指数函数值
uchar i;
for(i=0;i<IndexNumber;i++)
{
m=BaseNumber*m;
}
return m;
}
/******************************************************************************************
Pow函数结束：
******************************************************************************************/
/******************************************************************************************
函数说明：对数运算
参数说明：BaseNumber底数,AntiNumber真数
返回值：无符号整数
******************************************************************************************/
uint Log(uchar BaseNumber,uint AntiNumber)
{
uint m=0;//对数函数值
while(1)
{
AntiNumber=AntiNumber/BaseNumber;
if(AntiNumber)m++;
else break;
}
return m;
}
/******************************************************************************************
log函数结束：
******************************************************************************************/
/******************************************************************************************
函数说明：按时间抽取基二快速傅里叶算法
参数说明：N点数，等于2的M次方，*x序列的首地址
返回值：无
******************************************************************************************/
void Dit2Fft(uint N,Complex *x)
{
int i;
Complex temp;//临时复数变量
uint L,M=Log(2,N);//M级蝶形图
uint k,j;
uint StepLength;//k值步长
uint Bank ;//两点间的距离
const double pai=3.1415926;
double ps;
uint r;//旋转因子的幂
Complex W;//旋转因子
for(L=1;L<=M;L++)
{
StepLength=Pow(2,L);
Bank=StepLength/2;
for(j=0;j<=Bank-1;j++)
{
r=Pow(2,M-L)*j;
ps=2*pai/N*r;
W.Re=cos(ps);
W.Im=-sin(ps);
for(k=j;k<=N-1;k=k+StepLength)
{
Complex x_temp;
x_temp=Mul(W,x[k+Bank]);
temp=Plus(x[k],x_temp);
x[k+Bank]=Sub(x[k],x_temp);
x[k].Re=temp.Re;
x[k].Im=temp.Im;
}
}
}
}
/******************************************************************************************
Dit2Fft函数结束：
******************************************************************************************/
/******************************************************************************************
函数说明：傅里叶反变换算法
参数说明：N点数，等于2的M次方，*x序列的首地址
返回值：无
1、求旋转因子的共轭，再乘1/N
2、x(n)=IFFT[x(k)]=1/N*{FFT[x*(k)]}*
******************************************************************************************/
void Ifft1(uint N,Complex *x)
{
int i;
Complex temp;//临时复数变量
uint L,M=Log(2,N);//M级蝶形图
uint k,j;
uint StepLength;//k值步长
uint Bank ;//两点间的距离
const double pai=3.1415926;
double ps;
uint r;//旋转因子的幂
Complex W;//旋转因子
Invert_order(N,x);//位倒序输入
for(L=1;L<=M;L++)
{
StepLength=Pow(2,L);
Bank=StepLength/2;
for(j=0;j<=Bank-1;j++)
{
r=Pow(2,M-L)*j;
ps=2*pai/N*r;
W.Re=cos(ps);
W.Im=sin(ps);
for(k=j;k<=N-1;k=k+StepLength)
{
Complex x_temp;
x_temp=Mul(W,x[k+Bank]);
temp=Plus(x[k],x_temp);
x[k+Bank]=Sub(x[k],x_temp);
x[k].Re=temp.Re;
x[k].Im=temp.Im;
}
}
}
for(k=0;k<N;k++)
{
x[k].Re=x[k].Re/N;
x[k].Im=x[k].Im/N;
}
}
void Ifft2(uint N,Complex *x)
{
uint k;
Invert_order(N,x);//位倒序输入
for(k=0;k<N;k++)
{
x[k].Re=x[k].Re;
x[k].Im=-x[k].Im;
}
Dit2Fft(N,x);
for(k=0;k<N;k++)
{
x[k].Re=x[k].Re;
x[k].Im=-x[k].Im;
}
for(k=0;k<N;k++)
{
x[k].Re=x[k].Re/N;
x[k].Im=x[k].Im/N;
}
}
/******************************************************************************************
Ifft函数结束：
******************************************************************************************/
评论0