fft.zip_1D fft_vc fft_visual c

#ifndef FFT_H #define FFT_H #include <complex> // complex<float> #include <stdio.h> #include <iostream> #include <string.h> #include <math.h> #include<time.h> #include<stdlib.h> #include <new> //#include <time> using namespace std; typedef complex<float> Comp; // 复数类型定义 const float PI2 = 2.0f * 3.14159265f; // 常数2PI定义 const int MAX_N = 128; // 最大DFT点数 /*----*----*----*----*----*----*----*----*----*----*----*----* FFT算法模块接口定义 *----*----*----*----*----*----*----*----*----*----*----*----*/ /////////////////////////////////////////// // Function name : BitReverse // Description : 二进制倒序操作 // Return type : int // Argument : int src 待倒读的数 // Argument : int size 二进制位数 const int N = 1024; inline void swap (float &a, float &b) { float t; t = a; a = b; b = t; } void bitrp (float *xreal, float *ximag, int n) { // 位反转置换 Bit-reversal Permutation int i, j, a, b, p; for (i = 1, p = 0; i < n; i *= 2) { p ++; } for (i = 0; i < n; i ++) { a = i; b = 0; for (j = 0; j < p; j ++) { b = (b << 1) + (a & 1); // b = b * 2 + a % 2; a >>= 1; // a = a / 2; } if ( b > i) { swap (xreal [i], xreal [b]); swap (ximag [i], ximag [b]); } } } void FFT(float *xreal, float *ximag, int n) { // 快速傅立叶变换，将复数 x 变换后仍保存在 x 中，xreal, ximag 分别是 x 的实部和虚部 float wreal [N / 2], wimag [N / 2], treal, timag, ureal, uimag, arg; int m, k, j, t, index1, index2; bitrp (xreal, ximag, n); // 计算 1 的前 n / 2 个 n 次方根的共轭复数 W'j = wreal [j] + i * wimag [j] , j = 0, 1, , n / 2 - 1 arg = - 2 * PI / n; treal = float(cos (arg)); timag = float(sin (arg)); wreal [0] = 1.0; wimag [0] = 0.0; for (j = 1; j < n / 2; j ++) { wreal [j] = wreal [j - 1] * treal - wimag [j - 1] * timag; wimag [j] = wreal [j - 1] * timag + wimag [j - 1] * treal; } for (m = 2; m <= n; m *= 2) { for (k = 0; k < n; k += m) { for (j = 0; j < m / 2; j ++) { index1 = k + j; index2 = index1 + m / 2; t = n * j / m; // 旋转因子 w 的实部在 wreal [] 中的下标为 t treal = wreal [t] * xreal [index2] - wimag [t] * ximag [index2]; timag = wreal [t] * ximag [index2] + wimag [t] * xreal [index2]; ureal = xreal [index1]; uimag = ximag [index1]; xreal [index1] = ureal + treal; ximag [index1] = uimag + timag; xreal [index2] = ureal - treal; ximag [index2] = uimag - timag; } } } } void IFFT (float *xreal, float *ximag, int n) { // 快速傅立叶逆变换 float wreal [N / 2], wimag [N / 2], treal, timag, ureal, uimag, arg; int m, k, j, t, index1, index2; bitrp (xreal, ximag, n); // 计算 1 的前 n / 2 个 n 次方根 Wj = wreal [j] + i * wimag [j] , j = 0, 1, , n / 2 - 1 arg = 2 * PI / n; treal = float(cos (arg)); timag = float(sin (arg)); wreal [0] = 1.0f; wimag [0] = 0.0f; for (j = 1; j < n / 2; j ++) { wreal [j] = wreal [j - 1] * treal - wimag [j - 1] * timag; wimag [j] = wreal [j - 1] * timag + wimag [j - 1] * treal; } for (m = 2; m <= n; m *= 2) { for (k = 0; k < n; k += m) { for (j = 0; j < m / 2; j ++) { index1 = k + j; index2 = index1 + m / 2; t = n * j / m; // 旋转因子 w 的实部在 wreal [] 中的下标为 t treal = wreal [t] * xreal [index2] - wimag [t] * ximag [index2]; timag = wreal [t] * ximag [index2] + wimag [t] * xreal [index2]; ureal = xreal [index1]; uimag = ximag [index1]; xreal [index1] = ureal + treal; ximag [index1] = uimag + timag; xreal [index2] = ureal - treal; ximag [index2] = uimag - timag; } } } for (j=0; j < n; j ++) { xreal [j] /= n; ximag [j] /= n; } } ////////////////////////////////////////////////////// // Function name : FFT_2D_Kernel // Description : 2D FFT核心算法 // Return type : void // Argument : Comp x[MAX_N][MAX_N] 二维数据 // Argument : int M 幂次 // Argument : int flag 正逆变换标记 //以下本应由自己完成。 void FFT_2D(complex<float> **x, int n) { float Nx=n,Ny=n; // 先逐行进行 1D-FFT Comp *row; row=(complex<float> *) malloc (sizeof(complex<float>)*Nx); float *row_real,*row_imag; row_real=(float *) malloc (sizeof(float)*Nx); row_imag=(float *) malloc (sizeof(float)*Nx); for (int i=0; i<Ny; i++) // 取得第i行的数据 { row[i]=complex<float>(0.0f,0.0f); row_real[i]=0.0f; row_imag[i]=0.0f; } for (int i=0; i<Ny; i++) // 取得第i行的数据 { for (int j=0; j<Nx; j++) { row[j] = x[i][j]; row_real[j]=real(row[j]); row_imag[j]=imag(row[j]); } FFT(row_real, row_imag, n); // <--- 计算结果再存入矩阵x中 for ( int j=0; j<Nx; j++) { row[j]=complex<float>(row_real[j],row_imag[j]); x[i][j] = row[j]; } } // 再逐列进行 1D-FFT complex<float> *col; col=(complex<float> *) malloc (sizeof(complex<float>)*Ny); float *col_real,*col_imag; col_real=(float *) malloc (sizeof(float)*Nx); col_imag=(float *) malloc (sizeof(float)*Nx); for (int j=0; j<Nx; j++) { col[j]=complex<float>(0.0f,0.0f); col_real[j]=0.0f; col_imag[j]=0.0f; } for (int j=0; j<Nx; j++) { // 取得第j列的数据 for (int i=0; i<Ny; i++) { col[i] = x[i][j]; col_real[i]=real(col[i]); col_imag[i]=imag(col[i]); } // 对第j列数据进行 1D-FFT FFT(col_real, col_imag, n); // <--- 计算结果在数组col中 // 将结果放回矩阵第j列中 for (int i=0; i<Ny; i++) { col[i]=complex<float>(col_real[i],col_imag[i]); x[i][j] = col[i]; } } } void IFFT_2D(complex<float> **x, int n) { float Nx=n,Ny=n; // 先逐行进行 1D-FFT Comp *row; row=(complex<float> *) malloc (sizeof(complex<float>)*Nx); float *row_real,*row_imag; row_real=(float *) malloc (sizeof(float)*Nx); row_imag=(float *) malloc (sizeof(float)*Nx); for (int i=0; i<Ny; i++) // 取得第i行的数据 { for (int j=0; j<Nx; j++) { row[j] = x[i][j]; row_real[j]=real(row[j]); row_imag[j]=imag(row[j]); } IFFT(row_real, row_imag, n); // <--- 计算结果再存入矩阵x中 for ( int j=0; j<Nx; j++) { row[j]=complex<float>(row_real[j],row_imag[j]); x[i][j] = row[j]; } } // 再逐列进行 1D-FFT complex<float> *col; col=(complex<float> *) malloc (sizeof(complex<float>)*Ny); float *col_real,*col_imag; col_real=(float *) malloc (sizeof(float)*Nx); col_imag=(float *) malloc (sizeof(float)*Nx); for (int j=0; j<Nx; j++) { // 取得第j列的数据 for (int i=0; i<Ny; i++) { col[i] = x[i][j]; col_real[i]=real(col[i]); col_imag[i]=imag(col[i]); } // 对第j列数据进行 1D-FFT IFFT(col_real, col_imag, n); // <--- 计算结果在数组col中 // 将结果放回矩阵第j列中 for (int i=0; i<Ny; i++) { col[i]=complex<float>(col_real[i],col_imag[i]); x[i][j] = col[i]; } } } // <--- End of [FFT.H] #endif