打印

快速傅里叶变换(FFT)c语言实现

[复制链接]
654|1
手机看帖
扫描二维码
随时随地手机跟帖
跳转到指定楼层
楼主
linfelix|  楼主 | 2023-9-28 10:00 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>


#define PI 3.1415926
struct complex
{
    double real;
    double image;
};
struct complex complex_add(struct complex c1,struct complex c2);
struct complex complex_sub(struct complex c1,struct complex c2);
struct complex complex_multi(struct complex c1,struct complex c2);
struct complex rotation_factor(int N,int n,int k);
double mold_length(struct complex c);
void fft(int len, struct complex in_x[],struct complex out_y[]);




int main()
{
    int sam[8] = {1,-1,1,-1,1,-1,1,-1};
    int jhg[8] = {0,0,0,0,0,0,0,0};
    struct complex x[8];
    struct complex y[8];

    for(int i=0;i<8;i++)
    {
        x[i].real = sam[i];
        x[i].image = jhg[i];
    }
    printf("时域序列\n");

    for(int i=0;i<8;i++)
    {
        printf("(%.2f, %.2fi) \n",x[i].real,x[i].image);
    }

    fft(8,x,y);

    printf("频域序列\n");

    for(int i=0;i<8;i++)
    {
        printf("(%.2f, %.2fi)\n",y[i].real,y[i].image);
    }

   

    return 0;
}

struct complex complex_add(struct complex c1,struct complex c2)     //复数加法
{
    struct complex p;

    p.real = c1.real + c2.real;
    p.image = c1.image + c2.image;

    return p;

}

struct complex complex_sub(struct complex c1,struct complex c2)     //复数减
{
    struct complex p;

    p.real = c1.real - c2.real;
    p.image = c1.image - c2.image;

    return p;

}

struct complex complex_multi(struct complex c1,struct complex c2)  //复数乘法
{
    struct complex c3;
    c3.real=c1.real*c2.real - c1.image *c2.image;
    c3.image = c2.real* c1.image + c1.real*c2.image;

    return c3;
};



struct complex rotation_factor(int N,int n,int k)       //旋转因子
{
    struct complex w;

    w.real = cos(2* PI * n * k /N);
    w.image = - sin(2* PI * n * k /N);

    return w;
}

double mold_length(struct complex c)        //幅度
{
    return sqrt(c.real * c.real + c.image * c.image);
};


int reverse_num(int l,int oringin_num)          //反位序
{
    int q=0,m=0;

    for(int k=l-1;k>=0;k--)
    {
        q = oringin_num % 2;
        m += q*(1<<k);
        oringin_num = oringin_num / 2;
    }

    return m;
}


void fft(int len, struct complex in_x[],struct complex out_y[])
{
    /*
        param len 序列长度,目前只能是2的指数
        param in_x输入的序列
        param out_y输出的序列
    */

    int l,k,r,z,dist,q,j;       //l是级
    struct complex w,tmp;
    struct complex in_x_mem[len];

    l = log2(len);

    for(k=0;k<len;k++)
    {
        in_x_mem[k] = in_x[reverse_num(l,k)];       //反位序号操作
    }

    for(r = 0;r<l;r++)      //遍历每一级蝶形运算
    {

        dist = 1<<r;        //提前计算每一级的间隔距离

        for( j=0;j<dist;j++ )      //计算策略是拆成上下两组,先上计算,后下计算,j是计算的起始序号
        {
            for(k=j;k<len;k+=(dist<<1))     //不好解释,得画图理解
            {
                q = k+dist; //q同一组蝶形运算第二个序号
                z = k << (l - r -1);    //确定旋转因子的指数

                w = rotation_factor(len,1,z);
                //由于不是并行计算,必须先缓存
                tmp = in_x_mem[k];

                in_x_mem[k] = complex_add( in_x_mem[k] , complex_multi(w, in_x_mem[q]) );
                in_x_mem[q] = complex_sub(tmp , complex_multi(w, in_x_mem[q]) );

            }
        }
    }
    memcpy(out_y,in_x_mem,len*sizeof(struct complex));
}




使用特权

评论回复
沙发
AloneKaven| | 2023-9-29 22:13 | 只看该作者
还是没看懂是怎么实现的

使用特权

评论回复
发新帖 我要提问
您需要登录后才可以回帖 登录 | 注册

本版积分规则

38

主题

1274

帖子

1

粉丝