1. #include <stdio.h>
    2. #include <math.h>  
    3. #include <stdlib.h>  
    5. #define N 512                       //信号长度
    6. #define M 23                        //单位脉冲响应的长度  
    7. #define size_x 512                  //fft点数， 2的次幂  
    8. #define PI  3.1415926535            //圆周率 
    9. double HI[512] = { -0.00491000678664519, -0.0110940581309966,    0.00410194747948791,    0.0123959008767979,    0.0260846944086577,    0.0280068857178369,    0.0117031929574722, -0.0282213909190478, -0.0862193226866780, -0.148431590883831, -0.196291852743501,    0.785720061596288, -0.196291852743501, -0.148431590883831, -0.0862193226866780, -0.0282213909190478,    0.0117031929574722,    0.0280068857178369,    0.0260846944086577,    0.0123959008767979,    0.00410194747948791, -0.0110940581309966, -0.00491000678664519 };
    10. /*定义复数类型*/
    11. typedef struct {
    12.     double real;
    13.     double img;
    14. }complex;
    16. int Lenx = 4096;
    17. int M1 = M - 1;                     // 重叠部分的长度
    18. int L = N - 22;                     //不重叠部分的长度
    20. complex x[4096], * W;               //输入序列,变换核
    21. complex y[4410];                    //输出序列
    22. complex h[size_x], * W1;            //输入滤波器,变换核 
    23. complex x1[9][size_x], * W2;        //输入分段,变换核 
    24. complex y2[9][size_x], * W5;        //输出结果,变换核
    27. /*快速傅里叶变换*/
    28. void fft(complex h[512]);
    29. void ifft_y1(complex y_1[512]);
    31. /*初始化变换核*/
    32. void initW1();
    33. void initW2();
    35. /*变址*/
    36. void change1(complex H[512]);
    37. void change5(complex y_1[512]);
    40. void add(complex, complex, complex*); //复数加法  
    41. void mul(complex, complex, complex*); //复数乘法  
    42. void sub(complex, complex, complex*); //复数减法  
    43. void output();                        //输出快速傅里叶变换的结果  
    46. int main()
    47. {
    48.     int i;
    50.     //输出结果  
    51.     system("cls");
    52.     printf("输出滤波后的结果\n");
    54.     //输入序列
    55.     for (i = 0; i < 4096; i++)
    56.     {
    57.         x[i].real = sin(2 * PI * 150 * i / 8000);       //150Hz    
    58.         x[i].real += sin(2 * PI * 200 * i / 8000);      //200Hz
    59.         x[i].real += sin(2 * PI * 2000 * i / 8000);     //2000Hz
    60.     }
    62.     //滤波器    
    63.     for (i = 0; i < size_x; i++)
    64.     {
    65.         if (i < M)
    66.         {
    67.             h[i].real = HI[i];
    68.         }
    69.         if (i >= M)
    70.         {
    71.             h[i].real = 0;
    72.         }
    73.         h[i].img = 0;
    74.     }
    75.     int K = floor((Lenx + M1 - 1) / L) + 1;  //分的段数
    76.     int p = (K)*L - Lenx;                    // 最后一段结尾补零数
    77.     double zong[4432];
    78.     for (i = 0; i < 4432; i++)
    79.     {
    80.         if (i < M1)
    81.         {
    82.             zong[i] = 0;
    83.         }
    84.         if (i >= M1 && i < 4096)
    85.         {
    86.             zong[i] = x[i + M1].real;
    87.         }
    88.         else
    89.         {
    90.             zong[i] = 0;
    91.         }
    92.         h[i].img = 0;
    93.     }
    94.     double Y[9][512] = { 0 };
    95.     //输入函数的分段(重叠保留法)
    96.     for (int j = 0; j < 512; j++)
    97.     {
    98.         for (int k = 0; k < 9; k++)
    99.         {
    100.             x1[k][j].real = zong[k * L + j];
    101.         }
    104.     }
    105.     initW1();//调用变换核 
    106.     initW2();//调用变换核 
    107.     //对h进行DFT
    108.     fft(h); //调用快速傅里叶变换 
    109.     for (int k = 0; k < 9; k++)
    110.     {
    111.         fft(x1[k]);
    112.     }
    114.     //h的DFT分别与x1,x2,x3的DFT相乘
    115.     for (int k = 0; k < 9; k++)
    116.     {
    117.         for (i = 0; i < size_x; i++)
    118.         {
    119.             y2[k][i].real = x1[k][i].real * h[i].real - x1[k][i].img * h[i].img;
    120.             y2[k][i].img = x1[k][i].img * h[i].real + x1[k][i].real * h[i].img;
    121.         }
    122.     }
    123.     //对乘积求IDFT, 得到子序列和单位脉冲响应的循环卷积
    124.     //调用快速傅里叶变换
    125.     for (int k = 0; k < 9; k++)
    126.     {
    127.         ifft_y1(y2[k]);
    128.     }
    129.     //把每一段循环卷积数据前端的M-1个点去掉之后衔接起来
    130.     for (i = 0; i < size_x; i++)
    131.     {
    132.         for (int k = 0; k < 9; k++)
    133.         {
    134.             y2[k][i].real = y2[k][i].real / size_x;
    135.             y2[k][i].img = y2[k][i].img / size_x;
    136.         }
    137.     }
    139.     for (i = 0; i < 4410; i++)
    140.     {
    141.         for (int k = 0; k < 9; k++)
    142.         {
    143.             if (i < (k + 1) * L && i > k * L)
    144.             {
    145.                 y[i].real = y2[k][i + M - L * k - 1].real;
    146.                 y[i].img = y2[k][i + M - L * k - 1].img;
    147.             }
    148.         }
    150.     }
    151.     output();//结果输出 
    152.     system("pause");
    153.     return 0;
    154. }
    157. //h的fft
    158. void fft(complex h[512])
    159. {
    160.     int i = 0, j = 0, k = 0, l = 0;
    161.     complex up, down, product;
    162.     change1(h);  //调用变址函数  
    163.     for (i = 0; i < log(size_x) / log(2); i++)        /*一级蝶形运算 stage */
    164.     {
    165.         l = 1 << i;
    166.         for (j = 0; j < size_x; j += 2 * l)     /*一组蝶形运算 group,每个group的蝶形因子乘数不同*/
    167.         {
    168.             for (k = 0; k < l; k++)        /*一个蝶形运算 每个group内的蝶形运算*/
    169.             {
    170.                 mul(h[j + k + l], W1[size_x * k / 2 / l], &product);
    171.                 add(h[j + k], product, &up);
    172.                 sub(h[j + k], product, &down);
    173.                 h[j + k] = up;
    174.                 h[j + k + l] = down;
    175.             }
    176.         }
    177.     }
    178. }
    180. //y1的ifft
    181. void ifft_y1(complex y_1[512])
    182. {
    184.     int i = 0, j = 0, k = 0, l = 0;
    185.     complex up, down, product;
    186.     change5(y_1);  //调用变址函数  
    187.     for (i = 0; i < log(size_x) / log(2); i++)        /*一级蝶形运算 stage */
    188.     {
    189.         l = 1 << i;
    190.         for (j = 0; j < size_x; j += 2 * l)     /*一组蝶形运算 group,每个group的蝶形因子乘数不同*/
    191.         {
    192.             for (k = 0; k < l; k++)        /*一个蝶形运算 每个group内的蝶形运算*/
    193.             {
    194.                 mul(y_1[j + k + l], W2[size_x * k / 2 / l], &product);
    195.                 add(y_1[j + k], product, &up);
    196.                 sub(y_1[j + k], product, &down);
    197.                 y_1[j + k] = up;
    198.                 y_1[j + k + l] = down;
    199.             }
    200.         }
    201.     }
    202. }
    204. /*初始化变换核，定义一个变换核，相当于旋转因子WAP*/
    206. void initW1()
    207. {
    208.     int i;
    209.     W1 = (complex*)malloc(sizeof(complex) * size_x);  //生成变换核  
    210.     for (i = 0; i < size_x; i++)
    211.     {
    212.         W1[i].real = cos(2 * PI / size_x * i);   //用欧拉公式计算旋转因子  
    213.         W1[i].img = -1 * sin(2 * PI / size_x * i);
    214.     }
    215. }
    217. void initW2()
    218. {
    219.     int i;
    220.     W2 = (complex*)malloc(sizeof(complex) * size_x);  //生成变换核  
    221.     for (i = 0; i < size_x; i++)
    222.     {
    223.         W2[i].real = cos(2 * PI / size_x * i);   //用欧拉公式计算旋转因子  
    224.         W2[i].img = 1 * sin(2 * PI / size_x * i);
    225.     }
    226. }
    231. /*变址计算，将输入序列的码位倒置*/
    232. void change1(complex H[512])
    233. {
    234.     complex temp;
    235.     unsigned short i = 0, j = 0, k = 0;
    236.     double t;
    237.     for (i = 0; i < size_x; i++)
    238.     {
    239.         k = i; j = 0;
    240.         t = (log(size_x) / log(2));
    241.         while ((t--) > 0)    //利用按位与以及循环实现码位颠倒  
    242.         {
    243.             j = j << 1;
    244.             j |= (k & 1);
    245.             k = k >> 1;
    246.         }
    247.         if (j > i)    //将x(n)的码位互换  
    248.         {
    249.             temp = H[i];
    250.             H[i] = H[j];
    251.             H[j] = temp;
    252.         }
    253.     }
    254. }
    256. void change5(complex y_1[512])
    257. {
    258.     complex temp;
    259.     unsigned short i = 0, j = 0, k = 0;
    260.     double t;
    261.     for (i = 0; i < size_x; i++)
    262.     {
    263.         k = i; j = 0;
    264.         t = (log(size_x) / log(2));
    265.         while ((t--) > 0)    //利用按位与以及循环实现码位颠倒  
    266.         {
    267.             j = j << 1;
    268.             j |= (k & 1);
    269.             k = k >> 1;
    270.         }
    271.         if (j > i)    //将x(n)的码位互换  
    272.         {
    273.             temp = y_1[i];
    274.             y_1[i] = y_1[j];
    275.             y_1[j] = temp;
    276.         }
    277.     }
    279. }
    282. /*输出结果*/
    283. void output()
    284. {
    285.     int i;
    286.     FILE* file;
    287.     if ((file = fopen ("D:\\q.txt","w+" )) == NULL)//打开操作不成功
    288.     {
    289.             printf("The file can not be opened.\n");
    290.             exit(1);//结束程序的执行
    291.     }
    292.     printf("The result are as follows：\n");
    293.     for (i = 0; i < 4410; i++)
    294.     {
    295.         printf("%lf", y[i].real);
    296.         //if (y[i].img >= 0.0001)printf("+%lfj\n", y[i].img);
    297.         //else if (fabs(y[i].img) < 0.0001)printf("\n");
    298.         //else printf("%lfj\n", y[i].img);
    299.         fprintf(file, "%f,", y[i].real);
    300.     }
    301.     fclose(file);
    302. }
    305. void add(complex a, complex b, complex* c)  //复数加法的定义  
    306. {
    307.     c->real = a.real + b.real;
    308.     c->img = a.img + b.img;
    309. }
    312. void mul(complex a, complex b, complex* c)  //复数乘法的定义  
    313. {
    314.     c->real = a.real * b.real - a.img * b.img;
    315.     c->img = a.real * b.img + a.img * b.real;
    316. }
    319. void sub(complex a, complex b, complex* c)  //复数减法的定义  
    320. {
    321.     c->real = a.real - b.real;
    322.     c->img = a.img - b.img;
    323. }