fft代码

fft代码

#include <stdio.h>
#define N 128;
#define PI 3.1416;
typedef struct
{
    double realpart;
    double imaginaryPart;
}complexNumber;
complexNumber W[N>>1];
complexNumber A[N];
int Array[N];

// 输出数组
void OutDisplayArray()
{
    printf("输出数组\n")
    int i = 0;
    for(i=0;i<N;i++)
    {
        printf("%f %f\n",A[i].realpart,A[i].imaginaryPart);
    }
}
//创造N个采样值数组
void MakeArray()
{
    printf("输入%d位采样数组:",N);
    int i =0;
    for (i=0;i<N;i++)
    {
        print("输入第%d位采样数组:\n",i);
        scanf("%f",&A[i].realpart,&A[i].imaginaryPart);
    }
    OutDisplayArray();

}

//对数组进行倒序置换
void ChangeArray()
{
    int k=0,j=0;
    complexNumber tmp;
    for(int i = 0;i<N;i++)
    {
        if(j>i)
        {
            tmp = A[i]
            A[i] = A[j];
            A[j] = temp;
        }
        k = N >> 1;
        while(j >= k)
        {
            j = j-k;
            k = k >> 1;
        }
        j= j + k;
    }
    OutDisplayArray();
}

//快速傅里叶变换 对采样数组进行蝶形运算 
void fft()
{
    int s = lg(N);
    int i = 0;
    int ss = 0;
    complexNumber up,down,product,w,t,u;
    for(i=0;i<s;i++)
    {
        m= 1 << (i+1);                                //蝶形运算基本单元 
        ss= (N >> (i+1));                             // 蝶形运算的组数                                   
        for(k=0; k<ss; k++)                           // 蝶形运算的组数循环
        {
            w->realpart = 1;
            w->imaginaryPart = 0;

            for(j=0;j < (m>>1);j++)                   //一组蝶形运算基本单元循环
            {
                mul(w,A[k*m +j + m>>1],&product);
                t=product;
                u=A[k*m +j];
                add(u,t,&up);
                sub(u,t,&down);
                A[k*m +j] = up;
                A[k*m +j + m>>1] = down;

                         
                w->realpart = cos(2*pi*(j+1)/m);        //旋转因子     
                w->imaginaryPart = -sin(2*pi*(j+1)/m);
            }

        }
    }

}
/*复数加法的定义*/ 
void add(complexNumber a,complexNumber b,complexNumber *c)
{
	c->realPart=a.realPart+b.realPart;
	c->imaginaryPart=a.imaginaryPart+b.imaginaryPart;
}
/*复数减法的定义*/ 
void sub(complexNumber a,complexNumber b,complexNumber *c)
{
	c->realPart=a.realPart-b.realPart;
	c->imaginaryPart=a.imaginaryPart-b.imaginaryPart;
}
/*复数乘法的定义*/ 
void mul(complexNumber a,complexNumber b,complexNumber *c)
{
	c->realPart=a.realPart*b.realPart - a.imaginaryPart*b.imaginaryPart;
	c->imaginaryPart=a.realPart*b.imaginaryPart + a.imaginaryPart*b.realPart;
}