求基2、基4、基8FFT(快速傅里叶变换)的c语言程序,要能运行得出来的

谢谢!来个基8fft就行
2025-03-18 01:04:09
推荐回答(2个)
回答1:

1.FFT:
// data为输入和输出的数据,N为长度
bool CFFT::Forward(complex *const Data, const unsigned int N)
{
   if (!Data || N < 1 || N & (N - 1))
      return false;
   //   排序
   Rearrange(Data, N);
   //   FFT计算:const bool Inverse = false
   Perform(Data, N);

   return true;
}
 2.IFFT:
// Scale 为是否缩放
bool CFFT::Inverse(complex *const Data, const unsigned int N,
   const bool Scale /* = true */)
{
   if (!Data || N < 1 || N & (N - 1))
      return false;
   //   排序
   Rearrange(Data, N);
   //   FFT计算,ture表示是逆运算
   Perform(Data, N, true);
   //   对结果进行缩放
   if (Scale)
      CFFT::Scale(Data, N);

   return true;
}
3.排序:
void CFFT::Rearrange(complex *const Data, const unsigned int N)
{
   //   Swap position
   unsigned int Target = 0;
   //   Process all positions of input signal
   for (unsigned int Position = 0; Position < N; ++Position)
   {
      //   Only for not yet swapped entries
      if (Target > Position)
      {
         //   Swap entries
         const complex Temp(Data[Target]);
         Data[Target] = Data[Position];
         Data[Position] = Temp;
      }
      //   Bit mask
      unsigned int Mask = N;
      //   While bit is set
      while (Target & (Mask >>= 1))
         //   Drop bit
         Target &= ~Mask;
      //   The current bit is 0 - set it
      Target |= Mask;
   }
}
 4.FFT计算:
void CFFT::Perform(complex *const Data, const unsigned int N,
   const bool Inverse /* = false */)
{
   const double pi = Inverse ? 3.14159265358979323846 : -3.14159265358979323846;
   //   Iteration through dyads, quadruples, octads and so on...
   for (unsigned int Step = 1; Step < N; Step <<= 1)
   {
      //   Jump to the next entry of the same transform factor
      const unsigned int Jump = Step << 1;
      //   Angle increment
      const double delta = pi / double(Step);
      //   Auxiliary sin(delta / 2)
      const double Sine = sin(delta * .5);
      //   Multiplier for trigonometric recurrence
      const complex Multiplier(-2. * Sine * Sine, sin(delta));
      //   Start value for transform factor, fi = 0
      complex Factor(1.);
      //   Iteration through groups of different transform factor
      for (unsigned int Group = 0; Group < Step; ++Group)
      {
         //   Iteration within group 
         for (unsigned int Pair = Group; Pair < N; Pair += Jump)
         {
            //   Match position
            const unsigned int Match = Pair + Step;
            //   Second term of two-point transform
            const complex Product(Factor * Data[Match]);
            //   Transform for fi + pi
            Data[Match] = Data[Pair] - Product;
            //   Transform for fi
            Data[Pair] += Product;
         }
         //   Successive transform factor via trigonometric recurrence
         Factor = Multiplier * Factor + Factor;
      }
   }
}
 5.缩放:
void CFFT::Scale(complex *const Data, const unsigned int N)
{
   const double Factor = 1. / double(N);
   //   Scale all data entries
   for (unsigned int Position = 0; Position < N; ++Position)
      Data[Position] *= Factor;
}

回答2:

#include
#include
#include

typedef struct{
double r;
double i;
}my_complex
;

//检查a是否为2的整数次方数
#define NOT2POW(a) (((a)-1)&(a)||(a)<=0)
//pi
#define MYPI 3.14159265358979323846

my_complex* fft(const my_complex* x, unsigned int len){
unsigned int ex=0,t=len;
unsigned int i,j,k;
my_complex *y;
double tr,ti,rr,ri,yr,yi;

if(NOT2POW(len)) return NULL; //如果失败,返回空指针
for(;!(t&1);t>>=1) ex++; //len应该等于2的ex次方

y=(my_complex*)malloc(len*sizeof(my_complex));
if(!y) return NULL;

//变址计算,库里-图基算法
for(i=0;i k=i;
j=0;
t=ex;
while((t--)>0){
j<<=1;
j|=k&1;
k>>=1;
}
if(j>=i){
y[i]=x[j];
y[j]=x[i];
}
}

//用变址后的y向量进行计算
for(i=0;i t=1< for(j=0;j for(k=0;k ti=-MYPI*k/t;
rr=cos(ti);
ri=sin(ti);

tr=y[j+k+t].r;
ti=y[j+k+t].i;

yr=rr*tr-ri*ti;
yi=rr*ti+ri*tr;

tr=y[j+k].r;
ti=y[j+k].i;

y[j+k].r=tr+yr;
y[j+k].i=ti+yi;
y[j+k+t].r=tr-yr;
y[j+k+t].i=ti-yi;
}
}
}

return y;
}

//以下为测试
int main()
{
int i,DATA_LEN;
my_complex *x,*y;

printf("基二FFT测试\n输入生成序列长度:");
scanf("%d",&DATA_LEN);
x=(my_complex*)malloc(DATA_LEN*sizeof(my_complex));

for(i=0;i x[i].r=i;
x[i].i=i-1;
}

printf("处理前...\n实部\t\t虚部\n");
for(i=0;i printf("%lf\t%lf\n",x[i].r,x[i].i);

y=fft(x,DATA_LEN);
if(!y){
printf("序列长度不为2的整数次方!\n");
return 0;
}

printf("处理后...\n实部\t\t虚部\n");
for(i=0;i printf("%lf\t%lf\n",y[i].r,y[i].i);

free(y);
free(x);

return 0;
}