#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N 1024
double sqrt(double n);
/*定义复数类型*/
typedef struct {
    double real;
    double img;
    double Amp;
    double Ang;
}complex;
float amp[N];
float ang[N];
float real[N];
complex x[N], *W; /*输入序列,变换核*/
int size_x = 0;      /*输入序列的大小，在本程序中仅限2的次幂*/
char unuseful;
double PI;        /*圆周率*/
void fft();    /*快速傅里叶变换*/
void initW();  /*初始化变换核*/
void change(); /*变址*/
void add(complex, complex, complex *); /*复数加法*/
void mul(complex, complex, complex *); /*复数乘法*/
void sub(complex, complex, complex *); /*复数减法*/
void output();/*输出快速傅里叶变换的结果*/
int main()
{
    int i=0;                            /*输出结果*/
    int tmp;
    PI = atan(1) * 4;
    printf("输出DIT方法实现的FFT结果\n");
    /*printf("Please input the size of x:\n");//输入序列的大小
    scanf("%d", &size_x);
    printf("Please input the data in x[N]:\n");//输入序列的实部和虚部
    for (i = 0;i<size_x;i++)
    {
        printf("请输入第%d个序列：", i);*/
    FILE *fp;
    if ((fp = fopen("C:\\", "r")) == NULL)
    {
        printf("Cannot open file, press any key to exit!\n");
        exit(1);
    }
    while (!feof(fp))
    {
        fscanf(fp,"%d,",&tmp);
        x[i].real = tmp;
        real[i] = (float)x[i].real;
        i++;
        if (i>1024) break;
    }
    fclose(fp);
    //for (i = 0;i < 1024;i++) {
        //printf("x[%d].real=%d ", i, (int)x[i].real);
        //printf("\n");
    //}
   
    FILE *fo;
    if ((fo = fopen("C:\\", "r")) == NULL)
    {
        printf("Cannot open file, press any key to exit!\n");
        exit(1);
    }
    i = 0;
    while (!feof(fo))
    {
        fscanf(fo,"%d,",&tmp);
        x[i].img = tmp;
        i++;
        if (i>1024) break;
    }
    fclose(fo);
    //for (i= 0;i < 1024;i++) {
        //printf("x[%d].img=%d ", i, (int)x[i].img);
        //printf("\n");
    printf怎么实现的//}
    for (i = 0;i < 1024;i++) {
        //real_temp = x[i].real;
        //img_temp = x[i].img;
        printf("%lf + %lf\n",x[i].real,x[i].img);
        //printf("%d %d +%d %d\n",x[i].real, x[i].img,&x[i].real,&x[i].img);
        //printf("%d + %d\n", *(int *)ptest++, *(int *)ptest++);
        //printf("%f,%f\n", x[i].real,x[i].img);
        //printf(buff);
        //printf("x[%d].real+x[%d].img=%d+%d", i, i, x[i].img, x[i].real);
        //printf("%d", x[i].img);
        //printf("\n");
    }
    printf("输出倒序后的序列\n");
    initW();//调用变换核
    fft();//调用快速傅里叶变换
    printf("输出FFT后的结果\n");
    output();//调用输出傅里叶变换结果函数
    for (i = 0;i < 1024;i++)
    {
        x[i].Amp =sqrt( pow(x[i].real, 2) + pow(x[i].img, 2));
        x[i].Ang = atan(x[i].img / x[i].real);
        amp[i] = (float)x[i].Amp;
        ang[i] = (float)x[i].Ang;
        //printf("%lf %lf", j[i],m[i]);
        //printf("\n");
    }
    return 0;
}
/*快速傅里叶变换*/
void fft()
{
    int i = 0, j = 0, k = 0, l = 0;
    complex up, down, product;
    change();  //调用变址函数
    for (i = 0;i< log(1024) / log(2);i++)        /*一级蝶形运算 stage */
    {
        l = 1 << i;
        for (j = 0;j<1024;j += 2 * l)    /*一组蝶形运算 group,每个group的蝶形因子乘数不同*/
        {
            for (k = 0;k<l;k++)        /*一个蝶形运算 每个group内的蝶形运算*/
            {
                mul(x[j + k + l], W[1024*k / 2 / l], &product);
                add(x[j + k], product, &up);
                sub(x[j + k], product, &down);
                x[j + k] = up;
                x[j + k + l] = down;
            }
        }
    }
}
/*初始化变换核，定义一个变换核，相当于旋转因子WAP*/
void initW()
{
    int i;
    W = (complex *)malloc(sizeof(complex) * 1024);  //生成变换核
    for (i = 0;i<1024;i++)
    {
        W[i].real = cos(2 * PI / 1024*i);  //用欧拉公式计算旋转因子
        W[i].img = -1 * sin(2 * PI / 1024*i);
    }
}
/*变址计算，将x(n)码位倒置*/
void change()
{
    complex temp;
    unsigned short i = 0, j = 0, k = 0;
    double t;
    for (i = 0;i<1024;i++)
    {
        k = i;j = 0;
        t = (log(1024) / log(2));
        while ((t--)>0)    //利用按位与以及循环实现码位颠倒
        {
            j = j << 1;
            j |= (k & 1);
            k = k >> 1;
        }
        if (j>i)    //将x(n)的码位互换
        {
            temp = x[i];
            x[i] = x[j];
            x[j] = temp;
        }
    }
    output();
}
/*输出傅里叶变换的结果*/
void output()
{
    int i;
    printf("The result are as follows：\n");
    for (i = 0;i<1024;i++)
    {
        printf("%.4f+%.4fj\n", x[i].real, x[i].img);
    }
}
void add(complex a, complex b, complex *c)  //复数加法的定义
{
    c->real = a.real + b.real;
    c->img = a.img + b.img;
}
void mul(complex a, complex b, complex *c)  //复数乘法的定义
{
    c->real = a.real*b.real - a.img*b.img;
    c->img = a.real*b.img + a.img*b.real;
}
void sub(complex a, complex b, complex *c)  //复数减法的定义
{
    c->real = a.real - b.real;
    c->img = a.img - b.img;
}

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。