离散傅里叶变换C++代码

admin

https://blog.csdn.net/jfu22/article/details/50580777

1.
/*

2. * myfft.h

3. */

4.  

5. #ifndef __MYFFT_H__

6. #define __MYFFT_H__

7.  

8. #include <windows.h>

9.  

10.

11.typedef struct _my_complex

12.{

13.        double r; //复数实部

14.        double i; //复数虚部

15.

16.        _my_complex(){}

17.        _my_complex(double _r, double _i){r = _r; i = _i;}

18.        ~_my_complex(){}

19.

20.        _my_complex operator + (const _my_complex & c)//重载加号运算符

21.        {

22.                return _my_complex(r + c.r,i + c.i);

23.        }

24.        

25.        _my_complex operator - (const _my_complex & c)//重载减号运算符

26.        {

27.                return _my_complex(r - c.r,i - c.i);

28.        }

29.

30.        _my_complex operator * (const _my_complex & c)//重载乘号运算符

31.        {

32.                return _my_complex(r * c.r- i * c.i, r * c.i + i * c.r);

33.        }

34.}my_complex;

35.

36.

37.int myfft_by_define_c2c_1d(my_complex * in, int n1, my_complex *out, int n2, int line_width = 1, int is_inv = 0); //按照定义进行一维复数离散傅里叶变换

38.int myfft_by_define_c2c_2d(my_complex * in, int m, int n, int is_inv = 0); //按照定义进行二维复数离散傅里叶变换

39.

40.int myfft_c2c_1d_base2(my_complex * in, int n, int line_width = 1, int is_inv = 0); //进行一维复数离散快速傅里叶变换(基2)

41.

42.int print_complex_array_1d(my_complex * in, int n);

43.int print_complex_array_2d(my_complex * in, int m, int n);

44.

45.

46.#endif //__MYFFT_H__

1. /*

2. * myfft.cpp

3. */

4.  

5. #include <stdio.h>

6. #include <stdlib.h>

7. #include <math.h>

8. #include "myfft.h"

9.  

10.#define PI 3.1415926

11.

12.

13./*------------按照定义进行一维复数离散傅里叶变换----------------------------------------

14. * 原理：

15. *     一维离散傅里叶变换公式为： X[k] = ∑x[n]*e^(-2πi *nk/N)，其中n=0,1,2,...,N-1

16. *     一维离散傅里叶逆变换公式为： x[n] = 1/N * ∑X[k]*e^(2πi* nk/N)，其中k=0,1,2,...,N-1

17. *     欧拉公式：e^ix = cos(x) + isin(x)

18. *

19. * 参数1：in            输入的一维复数数组

20. * 参数2：n1            输入的一维复数数组的长度

21. * 参数3：out           输出的一维复数数组

22. * 参数4：n2            输出的一维复数数组的长度(必须和n1相等)

23. * 参数5：line_width    数组下标的间隔(默认为1)，对于行变换等于1，对于列变换等于一行的宽度

24. * 参数6：is_inv        是否是逆变换(默认不是)

25. *

26. * 过客& 386520874@qq.com & 2016.01.25

27. */

28.int myfft_by_define_c2c_1d(my_complex * in, int n1, my_complex *out, int n2, int line_width/*=1*/, int is_inv/*=0*/)

29.{

30.        if(n1 != n2){return 0;}

31.

32.        double _2PI_N = 2.0 * PI / n1;

33.

34.        if(is_inv == 0) //正变换

35.        {

36.                for(int m = 0; m < n1; m++)

37.                {

38.                         int M = m * line_width;

39.                         for(int k = 0; k < n1; k++ )

40.                         {

41.                                  int K = k * line_width;

42.                                  out[M].r+= in[K].r * cos(_2PI_N * m * k) + in[K].i * sin(_2PI_N * m * k); //计算实部

43.                                  out[M].i+= -in[K].r * sin(_2PI_N * m * k) + in[K].i * cos(_2PI_N * m * k); //计算虚部

44.                         }

45.                }

46.        }else if(is_inv == 1) //逆变换

47.        {

48.                for(int m = 0; m < n1; m++)

49.                {

50.                         int M = m * line_width;

51.                         for(int k = 0; k < n1; k++ )

52.                         {

53.                                  int K = k * line_width;

54.                                  out[M].r+= in[K].r * cos(_2PI_N * m * k) - in[K].i * sin(_2PI_N * m * k); //计算实部

55.                                  out[M].i+= in[K].r * sin(_2PI_N * m * k) + in[K].i * cos(_2PI_N * m * k); //计算虚部

56.                         }

57.                         out[M].r/= n1; //归一化

58.                         out[M].i/= n1; //归一化

59.                }

60.        }

61.

62.        return 1;

63.}

64.

65.

66.//按照定义进行二维复数离散傅里叶变换

67.int myfft_by_define_c2c_2d(my_complex * in, int m, int n, int is_inv/*=0*/)

68.{

69.        my_complex *temp = new my_complex[m * n];

70.        memset(temp, 0, sizeof(my_complex) * m *n);

71.

72.        for(int y = 0; y < n; y++) //先行变换

73.        {

74.                myfft_by_define_c2c_1d(in+ y * n, n, temp + y * n, n, 1, is_inv);

75.        }

76.

77.        memset(in, 0, sizeof(my_complex) * m *n);

78.        for(int x = 0; x < m; x++) //再列变换

79.        {

80.                myfft_by_define_c2c_1d(temp+ x, m, in + x, m, n, is_inv);

81.        }

82.

83.        delete [] temp; temp = NULL;

84.

85.        return 1;

86.}

87.

88.

89./*------------进行一维复数离散快速傅里叶变换(基2)----------------------------------------

90. * 原理：

91. *     对于FFT算法，一般是基2型的，即要求被变换的数组长度必须是2的m幂次方，

92. *     这种要求过于苛刻，当然相比于按照定义进行计算来说，速度是很快的。此外，

93. *     还有基4的算法，以及混合基算法等，但都对数组长度N有严格的限制。

94. *     所以想自己开发一个适合任意N的FFT算法，还是比较难的，目前还是调用FFTW库暂时凑合着用吧。

95. *

96. * 参数1：in            输入的一维复数数组

97. * 参数2：n             输入/出的一维复数数组的长度

98. * 参数3：line_width    数组下标的间隔(默认为1)，对于行变换等于1，对于列变换等于一行的宽度

99. * 参数4：is_inv        是否是逆变换(默认不是)

100.          *

101.          * 过客& 386520874@qq.com & 2016.01.25

102.          */

103.         int myfft_c2c_1d_base2(my_complex * in, int n, int line_width/*=1*/, int is_inv/*=0*/)

104.         {

105.                  int ttt = n;

106.                  int L = 1; //蝶形运算的层数

107.                  while(ttt = ttt >> 1){L++;}

108.                  

109.                  my_complex* X = in;

110.         

111.                  int N = n;

112.                  int len = 1 << L; //将输入的一维数组对齐为 N=2^L

113.                  if(n * 2 == len) //刚好是2的幂

114.                  {

115.                          L-= 1;

116.                  }else if(n < len) //说明n不是2的幂

117.                  {

118.         //               printf("Erorr:n is not 2^m format.\n"); //非基2算法还有待开发

119.         //               return0;

120.         

121.                          N =len;

122.                          X =new my_complex[N]; //扩充数组大小

123.                          int aaa = sizeof(my_complex);

124.                          memset(X, 0, sizeof(my_complex) * N); //全部初始化为0

125.                          memcpy(X, in, sizeof(my_complex) * n); //复制数据

126.                  }

127.         

128.                  //--------数组倒位序 二进制(n0n1n2) => (n2n0n1)-------------------

129.                  //例如： 3的二进制为 011，对应的倒位序为 110，即为6

130.                  my_complex* XX = new my_complex[N];

131.                  for(int i = 0; i < N; i++)

132.                  {

133.                          int p = 0;

134.                          for(int f = 0; f < L; f++)

135.                          {

136.                                   if(i & (1 << f))

137.                                   {

138.                                            p+= 1 << (L - f - 1);

139.                                   }

140.                          }

141.                          XX= X[p];

142.                  }

143.                  memcpy(X, XX, sizeof(my_complex) * N); //复制数据

144.                  delete [] XX; XX = NULL;

145.         

146.                  //---------旋转因子-------------------

147.                  my_complex* W = new my_complex[N / 2];

148.                  memset(W, 0, sizeof(my_complex) * N / 2); //全部初始化为0

149.                  for(int i = 0; i < N / 2; i++)

150.                  {

151.                          W= my_complex(cos(-2.0 * PI * i / N), sin(-2.0 * PI * i / N));

152.                  }

153.         

154.                  //---------FFT算法(基2)----------------

155.                  for(int f = 0; f < L; f++) //蝶形运算层数

156.                  {

157.                          int G = 1 << (L - f);

158.                          int num = 1 << f;

159.         

160.                          for(int u = 0; u < num; u++) //每组的元素个数

161.                          {

162.                                   for(int g = 0; g < G; g += 2) //每层的组数

163.                                   {

164.                                            my_complexX1 = X[g * num + u];

165.                                           my_complex X2 = X[(g + 1) * num + u];

166.         

167.                                            X[g* num + u]       = X1 + W[u * (1 << (L - f - 1))] * X2; //蝶形运算

168.                                            X[(g+ 1) * num + u] = X1 -W[u * (1 << (L - f - 1))] * X2; //蝶形运算

169.                                   }

170.                          }

171.                  }

172.         

173.                  if(X != NULL && X != in) //说明另外申请了内存空间

174.                  {

175.                          memcpy(in, X, sizeof(my_complex) * n);

176.                          if(X){delete [] X; X = NULL;}

177.                  }

178.                  if(W != NULL){delete [] W; W = NULL;}

179.         

180.                  return 1;

181.         }

182.         

183.         

184.         int print_complex_array_1d(my_complex * in, int n)

185.         {

186.                  for(int i = 0; i < n; i++)

187.                  {

188.                          printf("%d:%f+i%f", i, in.r, in.i);

189.                  }

190.                  printf("\n");

191.         

192.                  return 1;

193.         }

194.         

195.         

196.         int print_complex_array_2d(my_complex * in, int m, int n)

197.         {

198.                  for(int y = 0; y < n; y++)

199.                  {

200.                          for(int x = 0; x < m; x++)

201.                          {

202.                                  printf("%d_%d:%f+i%f", y, x, (*(in + y * m + x)).r, (*(in + y * m + x)).i);

203.                          }

204.                          printf("\n");

205.                  }

206.                  printf("\n");

207.         

208.                  return 1;

209.         }

1.       /*

2.       * main.cpp

3.       */

4.        

5.       #include <stdio.h>

6.       #include <stdlib.h>

7.       #include "myfft.h"

8.        

9.        

10.      int main(int argc, char * argv[]) //DFT测试 一维

11.      {

12.              int n = 8;

13.              int n1 = n;

14.              int n2 = n;

15.              my_complex * in = new my_complex[n1];

16.              my_complex * out = new my_complex[n2];

17.      

18.              double data[] = {4, 3, 2, 6, 7, 8, 9, 0};

19.      

20.              for(int i = 0; i < n; i++)

21.              {

22.                     in.r = data;in.i = 0;

23.                     out.r = 0; out.i = 0;

24.              }

25.      

26.              print_complex_array_1d(in,n1);

27.      

28.              //----------正变换----------------------

29.              myfft_by_define_c2c_1d(in,n1, out, n2, 1, 0); //按照定义进行一维复数离散傅里叶变换

30.      

31.              print_complex_array_1d(out,n2);

32.              

33.              //----------逆变换----------------------

34.              memset(in, 0, sizeof(my_complex) * n);

35.      

36.              myfft_by_define_c2c_1d(out,n2, in, n1, 1, 1); //按照定义进行一维复数离散傅里叶变换

37.      

38.              print_complex_array_1d(in,n1);

39.      

40.              if(in){delete [] in; in = NULL;}

41.              if(out){delete [] out; out = NULL;}

42.      

43.              system("pause");

44.      

45.              return 0;

46.      }

47.      /*

48.      0:4.000000+i0.000000 1:3.000000+i0.000000 2:2.000000+i0.0000003:6.000000+i0.000000 4:7.000000+i0.000000 5:8.000000+i0.0000006:9.000000+i0.000000 7:0.000000+i0.000000

49.      0:39.000000+i0.000000 1:-10.778175+i6.292892 2:0.000001+i-5.0000013:4.778176+i-7.707106 4:5.000000+i0.000001 5:4.778172+i7.7071086:-0.000002+i4.999998 7:-10.778169+i-6.292899

50.      0:4.000000+i-0.000001 1:3.000000+i-0.000001 2:2.000000+i-0.0000003:6.000000+i-0.000000 4:7.000000+i0.000000 5:8.000000+i0.0000006:9.000000+i0.000001 7:0.000000+i0.000001

51.      请按任意键继续. . .

52.      */

53.      

54.      

55.      int main2(int argc, char * argv[]) //DFT测试 二维

56.      {

57.              int m = 3;

58.              int n = 3;

59.              my_complex * in = new my_complex[m * n];

60.      

61.              in[0] = my_complex(0, 0); in[1] = my_complex(2, 0); in[2] = my_complex(4, 0);

62.              in[3] = my_complex(6, 0); in[4] = my_complex(1, 0); in[5] = my_complex(3, 0);

63.              in[6] = my_complex(5, 0); in[7] = my_complex(7, 0); in[8] = my_complex(4, 0);

64.      

65.              print_complex_array_2d(in, m,n);

66.      

67.              //----------正变换----------------------

68.              myfft_by_define_c2c_2d(in, m,n, 0); //按照定义进行二维复数离散傅里叶变换

69.              

70.              print_complex_array_2d(in, m,n);

71.              

72.              //----------逆变换----------------------

73.              myfft_by_define_c2c_2d(in, m,n, 1); //按照定义进行二维复数离散傅里叶变换

74.              

75.              print_complex_array_2d(in, m,n);

76.      

77.              if(in){delete [] in; in = NULL;}

78.      

79.              system("pause");

80.      

81.              return 0;

82.      }

83.      /*

84.      0_0:0.000000+i0.000000 0_1:2.000000+i0.000000 0_2:4.000000+i0.000000

85.      1_0:6.000000+i0.000000 1_1:1.000000+i0.000000 1_2:3.000000+i0.000000

86.      2_0:5.000000+i0.000000 2_1:7.000000+i0.000000 2_2:4.000000+i0.000000

87.        

88.      0_0:32.000000+i0.000000 0_1:0.500000+i0.866025 0_2:0.500001+i-0.866027

89.      1_0:-7.000001+i5.196152 1_1:-1.000000+i-1.732051 1_2:-8.499999+i-6.062178

90.      2_0:-6.999999+i-5.196154 2_1:-8.500001+i6.062177 2_2:-1.000000+i1.732051

91.        

92.      0_0:0.000000+i-0.000000 0_1:2.000000+i-0.000000 0_2:4.000000+i-0.000000

93.      1_0:6.000000+i-0.000000 1_1:1.000000+i-0.000000 1_2:3.000000+i0.000000

94.      2_0:5.000000+i-0.000000 2_1:7.000000+i0.000000 2_2:4.000000+i0.000001

95.        

96.      请按任意键继续. . .

97.      */

98.      

99.      

100.     int main3(int argc, char * argv[]) //FFT测试 一维(基2)

101.     {

102.             int n = 8;

103.             my_complex * in = new my_complex[n];

104.             double data[] = {4, 3, 2, 6, 7, 8, 9, 0};

105.     

106.             for(int i = 0; i < n; i++)

107.             {

108.                    in.r = data;

109.                    in.i = 0;

110.             }

111.     

112.             print_complex_array_1d(in,n);

113.     

114.             //----------正变换----------------------

115.             myfft_c2c_1d_base2(in, n, 1, 0); //进行一维复数离散快速傅里叶变换(基2)

116.     

117.             print_complex_array_1d(in,n);

118.     

119.             if(in){delete [] in; in = NULL;}

120.     

121.             system("pause");

122.     

123.             return 0;

124.     }

125.     /*

126.     0:4.000000+i0.000000 1:3.000000+i0.000000 2:2.000000+i0.0000003:6.000000+i0.000000 4:7.000000+i0.000000 5:8.000000+i0.0000006:9.000000+i0.000000 7:0.000000+i0.000000

127.     0:39.000000+i0.000000 1:-10.778175+i6.292893 2:0.000000+i-5.0000003:4.778175+i-7.707106 4:5.000000+i0.000000 5:4.778174+i7.7071076:-0.000000+i5.000000 7:-10.778175+i-6.292894

128.     请按任意键继续. . .