C语言快速傅里叶变化FFT

快速傅里叶变化：主要有几个难点
1.如果使用DIT首先要倒位序，网上的方法值得参考
思路是和二进制进位的思路相同只是是反着进位
2.要自己封装复数类型
3.难点是公式里的k不好求，因为k不是连续取得的，我用了一个数j来存储，这样可以把k计算出来
4.wn^r中r的取值，直接移位即可，不需要转化为二进制
5.输入不是2的倍数时自动补零
其余我都写在注释里了，感觉还是有很多地方要改进

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "复数类型.c"
//可以用M_E表示自然对数，M_PI表示圆周率

/*首先要导位序，也可以只要一个存储空间A[i],只需要改变A[i]中i
重点1.只改变当i小于j的时候的值，而i大于等于j时因为已经换过，不能再换回来了
重点2.思路是与二进制算法思路，小于N/2时即，最高位为0，此时的结果需要加上N/2
当大于N/2时，最高位为1，此时要向低位进位，意味着要先-N/2代表进位
再判断是否大于N/2，大于再向次低位进位，否则直接加N/4*/

void daoxu(double *x,int n)
{
    int i=0,j=0,k=0;//i表示原始值，而j表示倒位序后的值
    double temp;
    for(;i < n;i++)
    {
        if(i<j)
        {
            temp = x[j];
            x[j] = x[i];
            x[i] = temp;
        }

        k = n/2;
        while(j >= k && k != 0)
        {
            j = j - k;
            k = k/2;
        }
        j = j + k;     
    } 
    /*for(i=0;i<n;i++)
    {
        printf("x[%d]=%f\n",i,x[i]);
        //printf("result: new_k[%d]=%d\n",i,new_k[i]);
    }*/
}

int process(int* ori_len)//对输入的数据进行处理，还要对数据进行补零，问题在于如何分配内存，使用链表？
{
    int temp = 1,i;
    for(i=0;;i++)
    {
        if(temp >= *ori_len)
            break;
        temp =temp * 2;
    }
    *ori_len = temp;
    return i;//返回幂次方
}

complex w_n(int k,int m,int L,int length)
//如果这个值是复数意味着要包装复数类型，并且和这个数做四则运算的其他数也必须是复数
//需要新建复数这个数据结构，但是其实最后计算的是幅值
{
    int r=k<<(L-m);//不应该取这个k，而应该取倒位序的k
    r=r&(length-1);//去除高位
    //printf("得到的r为=%d\n",r);
    complex result;
    result.real=cos(2*M_PI/length*r);//用欧拉公式展开
    //printf("result.real=cos(2*pi/%d*%d)=%f\n",length,r,result.real);
    result.imag=-sin(2*M_PI/length*r);
    //printf("result.imag=-sin(2*pi/%d*%d)=%f\n",length,r,result.imag);
    return result;
}

void butterfly_operation(double *x,int powe,int length)
{
    int k_num = length/2,k_cal,k,l;
    complex r;
    complex *back_up=(complex*)calloc(length,sizeof(complex));//需要备份值来保留处理前的数据,定义为虚数类型
    complex *y=(complex*)calloc(length,sizeof(complex));
    
    for(int m=1;m<=powe;m++)//pow表示蝶形的次数
    {
        //printf("第%d次蝶形运算\n",m);
        //printf("要做%d次蝶形运算\n",k_num);
        k = 0;//运算之前将k清零,这里的k是实际是倒位序的
        l = 0;
        for(int i=0;i<length;i++)//运算完一次蝶形更新备份值
        {
            if(m==1)
                back_up[i] = cpx_make(x[i],0);
            else
                back_up[i] = y[i];
            //printf("back_up[i]=%f+%fi\n",back_up[i].real,back_up[i].imag);
        }
        k_cal=pow(2,m-1);//k_cal实际代表的是2^(m-1)的值
        //printf("计算2^(m-1)的值=%d\n",k_cal);
        for(int j=0;j<k_num;j++)//在做第m次蝶形的时候，要做k_num次蝶形运算
        {
            //printf("k=%d\n",k);
            r = w_n(k,m,powe,length);
            //printf("r=%f+%fi\n",r.real,r.imag);
            y[k] =cpx_add(back_up[k],cpx_mul(back_up[k+k_cal],r));//x[k]取的是模值
            //printf("y[%d]=%f+%fi\n",k,y[k].real,y[k].imag);
            y[k+k_cal] = cpx_sub(back_up[k],cpx_mul(back_up[k+k_cal],r)); 
            //printf("cpx_mul(back_up[k+k_cal],r)=%f+%fi\n",cpx_mul(back_up[k+k_cal],r).real,cpx_mul(back_up[k+k_cal],r).imag);
            //printf("y[%d]=%f+%fi\n",k+k_cal,y[k+k_cal].real,y[k_cal+k].imag);
            if(k_cal>=2&&k<l+k_cal-1)
            {
                k+=1;               
            }
            else
            {
                k=k+k_cal+1;
                l=k;
            }      
        }
        //printf("-------------------------------------------------------------------\n");
    }
    for(int i=0;i<length;i++)
    {
        x[i]=sqrt(y[i].real*y[i].real+y[i].imag*y[i].imag);
        //printf("y[i].real=%f,y[i].imag=%f\n",y[i].real,y[i].imag);
    }
    free(back_up);
    free(y);
}

int main()
{
    int length,power,old_length;//以后要更正
    double t;
    printf("请输入采样的长度length=");
    scanf("%d",&length);//输入采样的长度*/
    double *x,*y;
    y=(double*)calloc(length,sizeof(double));//可以动态分配内存
    for(int i=0;i<length;i++)
    {
        //y[i] = i+1;//在输入之前就会分配好内存空间因此不能临时输入
        t=i*0.001;
        y[i] = 1.2+2.7*cos(2*M_PI*33*t)+5*cos(2*M_PI*200*t+M_PI/2);
        //y[i]=0.6*sin(2*M_PI*500*i)+0.6*sin(2*M_PI*50*i);
        //y[i] = sin(M_PI*i/8);
        //printf("%f\n",y[i]);
    }
    old_length = length;
    power = process(&length);//改变了length    
    x=(double*)calloc(length,sizeof(double));   
    for(int i=0;i<old_length;i++)
    {
        x[i] = y[i];
        printf("x[i]=%f\n",x[i]);
    }
    for(int i=old_length;i<length;i++)
    {
        x[i]=0;
        printf("x[i]=%f\n",x[i]);
    }
    free(y);
    daoxu(x,length);
    butterfly_operation(x,power,length);
    for(int i=0;i<length;i++)
        printf("%d\t%f\n",i,x[i]);//在输入之前就会分配好内存空间因此不能临时输入
    free(x);
    return 0;
}

我自己封装的复数头文件

#include <stdio.h>
#define Real(c) (c).real
#define Imag(c) (c).imag

typedef struct
{
    double real;
    double imag;
    
}complex;

complex cpx_make(double real, double imag)
{
    complex ret;
    ret.real = real;
    ret.imag = imag;
    return ret;
}

complex cpx_add(complex a,complex b)
{
    return cpx_make(Real(a) + Real(b), Imag(a) + Imag(b));
}

complex cpx_sub(complex a,complex b)
{
    return cpx_make(Real(a) - Real(b), Imag(a) - Imag(b));
}

complex cpx_mul(complex a,complex b)
{
    return cpx_make(Real(a) * Real(b)-Imag(a) * Imag(b),Real(a) * Imag(b) + Imag(a) * Real(b));
}

complex cpx_div(complex a,complex b)
{
    Real(b) = 1/Real(b);
    Imag(b) = 1/Imag(b);
    return cpx_mul(a,b);
}