/*
项目名称:基2时间快速傅里叶变换
内容:n为点数,输入值为2的整数幂,l为变换的级数,第一级为最底层。wnk(n,k)输出W以N为周期,k的复数。
使用方法:将Array内容改为x(t)时域采样值,更改n为
时间:2021-11-16
*/
#include <stdio.h>
#include<complex.h>
#include<math.h>
#include<stdlib.h>
#define PI 3.14159265
complex* fft(int* array,int N,int L,int s);
int reversal(int num,int m) //倒序数组,如1倒置后为4,x[1]对应的是采样数组中下标为4的元素
{
int answer=0x00,i,temp=1;
for(i=1;i<=m;i++)
{
temp=num%2;
temp=temp<<(m-i);
num=num/2;
answer|=temp;
}
return answer;
}
complex wnk(int n,int k)//求复数e-j2pik/N
{
return cos(2*PI*k/n)-sin(2*PI*k/n)*1i;
}
int main(void)
{
int i=0,n=16,l=log2(n),s=0;
l=((pow(2,l)+0.0001)<n)?(l+1):l;
int Array[]={1,-2,3,4,5,6,7,8,9,10,11,12,13,14,15,16};//时域采样数组
int *array=malloc(((pow(2,l)+1)*sizeof(int)));
for(i=0;i<n;i++)
array[i]=Array[i];
n=pow(2,l)+0.0001;//自动选取2的整数次幂作为采样点数
for(;i<=n;i++)//将数组补零
array[i]=0;
complex *X;
X=fft(array,n,l,s);
for(i=0;i<n;i++)
{
if(cimag(X[i])>=0)
printf("第%2d个点:%f+%fi\n",i,creal(X[i]),cimag(X[i]));
else
printf("第%2d个点:%f%fi\n",i,creal(X[i]),cimag(X[i]));
}
return 0;
}
complex* fft(int* array,int n,int l,int s)
{
int i=0,m=log2(n);
if(l==1)//最后一级
{
complex *array3=malloc(3*sizeof(complex));
int x1=array[reversal(s,m)],x2=array[reversal(s+1,m)];
array3[0]=x1+x2*wnk(n,0);
array3[1]=x1-x2*wnk(n,0);
return array3;
}
else
{
int interval=(int)(pow(2,l-1)+0.0001);
complex *array1 =fft(array,n,l-1,s);
complex *array2 =fft(array,n,l-1,s+interval);//X[k]=X1[k]+X2[k]*Enk
complex *array3 =malloc((n+1)*sizeof(complex));//X[k+n/2]=X1[k]-X2[k]*Enk
for(i=0;i<interval;i++)
{
complex x1=array1[i];
complex x2=(complex)(array2[i])*(complex)(wnk(n,i*(int)(pow(2,m-l)+0.0001)));
array3[i]=x1+x2;
array3[i+interval]=x1-x2;
}
return array3;
}
}