diff --git a/README.md b/README.md new file mode 100644 index 0000000..ce165b9 --- /dev/null +++ b/README.md @@ -0,0 +1,3 @@ +## Fast Fourier Tranform +[FFT Code] (https://github.com/rorysroes/Fast_Algorithm_FFT/blob/Fast_Fourier_Transform/修正版mid.cpp): +This project I implemented Cooley–Tukey FFT algorithm with radix_2, radix_3 and radix_5 diff --git a/fft ver2.cpp b/fft ver2.cpp new file mode 100644 index 0000000..75543da --- /dev/null +++ b/fft ver2.cpp @@ -0,0 +1,431 @@ +#include +#include +#include +#include +int SFT(double *x_r, double *x_i, double *y_r, double *y_i, int N); +int FFT(double *x_r, double *x_i, double *y_r, double *y_i, int N); +int Generate_N(int p, int q, int r); +int Initial(double *x, double *y, int N); +int Print_Complex_Vector(double *x, double *y, int N); +int Bit_Increase(int *D,int b, int N); +int Bit_Reserve(int *D, int b, int N); +int Bit_Reserve_Integer(int N); +int FFTver3(double *x_r, double *x_i, double *y_r, double *y_i, int N); +int Group(int N); + +int main() +{ + int D[4]; + //Bit_Reserve(D, 5, 4); + //Bit_Reserve_Integer(125); + //Group(8); + /* + Bit_Increase(D, 2, 4); + Bit_Reserve(D, 2, 4); // 2 ¶i¦ì + */ + // y_k = sum(x_n * w^{-kn}, n=0..N-1) + // w = cos(2*pi/N)+isin(2*pi/N) + int k ,n ,N, p, q, r; + double *y_r, *y_i, *x_r, *x_i, w_r, w_i; + clock_t t1,t2; + //input N + printf("Please input p q r="); + scanf("%d %d %d", &p, &q, &r); + N = Generate_N(p, q, r); + printf("N=2^%d 3^%d 5^%d = %d\n", p, q, r, N); + //N = 1 << p; //²¾¦ìºâ¤l 2^p ,¥ª²¾ p ­Ó¦ì¸m(ex:p=4 => 1= 0001 => 1000 = 16) + //printf("N=%d\n",N); + //create memory for x and y + x_r = (double * ) malloc(N*sizeof(double)); + x_i = (double * ) malloc(N*sizeof(double)); + y_r = (double * ) malloc(N*sizeof(double)); + y_i = (double * ) malloc(N*sizeof(double)); + /* + Initial(x_r, x_i, N); + t1 = clock(); + SFT(x_r, x_i, y_r, y_i, N); + t2 = clock(); + printf("%f secs\n",1.0*(t2 -t1)/CLOCKS_PER_SEC); + //Print_Complex_Vector(y_r,y_i, N); + */ + Initial(x_r, x_i, N); + t1 = clock(); + FFT(x_r, x_i, y_r, y_i, N); + t2 = clock(); + printf("%f secs\n",1.0*(t2 -t1)/CLOCKS_PER_SEC); + //Print_Complex_Vector(y_r,y_i, N); + + Initial(x_r, x_i, N); + t1 = clock(); + FFTver3(x_r, x_i, y_r, y_i, N); + t2 = clock(); + printf("%f secs\n",1.0*(t2 -t1)/CLOCKS_PER_SEC); +} +int SFT(double *x_r, double *x_i, double *y_r, double *y_i, int N) +{ + int k, n; + double w_r, w_i; + for(k=0;k=0) printf("%d : %f +%f i\n", n, x[n], y[n]); + else printf("%d : %f %f i\n", n, x[n], y[n]); + } +} +int Bit_Increase(int *D, int b, int N) +{ + int i; + // D[3] D[2] D[1] D[0] + //print 0000 0001 0010 0011 0100 0101 ... + for(i=0;i=0;i--) printf("%d", D[i]); + printf("\n"); + D[0] = D[0] + 1; + //check every bit , if D[i] = b => D[i] = 0 , D[i+1]+1 0002 => 0010 + i = 0; + while(D[i]==b & i < N-1) + { + D[i] = 0; + i = i + 1; + D[i] = D[i] + 1; + } + system("pause"); + } + return 0; +} + +int Bit_Reserve(int *D, int b, int N) +{ + int i; + // D[3] D[2] D[1] D[0] + //print 0000 1000 0100 1100 0010 1010 ... + for(i=0;i=0;i--) printf("%d", D[i]); + printf("\n"); + + //check every bit , if D[i] = b => D[i] = 0 , D[i+1]+1 2000 => 0100 + i = N-1; + while(D[i]==b-1 & i > 0) + { + D[i] = 0; + i = i - 1; + //D[i] = D[i] + 1; + } + D[i] = D[i] + 1; + system("pause"); + } + return 0; +} +int Bit_Reserve_Integer(int N) +{ + + // N = 8 + // D[2] D[1] D[0] + // x print 000 100 010 110 001 101 ... + // print 0 4 2 6 1 5 + // N = 27 + // 0 9 18 3 12 ..... + // N = 125 + // 0 25 .... + int i=0, j=0, M ; + + while(i < N) + { + printf("%d <-> %d", i, j); + + M = N/5 ; //check 25 + + while( j >= 4 * M & M > 0) + { + j = j - 4 * M ; + M = M / 5 ; + } + j = j + M ; + i = i + 1 ; + + system("pause"); + } + + /*3ªºª©¥» + while(i < N) + { + printf("%d <-> %d", i, j); + + M = N/3 ; //check 9 + + while( j >= 2 * M & M > 0) + { + j = j - 2*M ; + M = M / 3 ; + } + j = j + M ; + i = i + 1 ; + + system("pause"); + } + */ + /*2ªºª©¥» + while(i< N) + { + printf("%d <-> %d", i, j); + + //check every bit , if D[i] = b => D[i] = 0 , D[i+1]+1 2000 => 0100 + //i = N-1; + M = N/2 ; + // while(D[i]==b-1 & i>0) + while( j >= M & M > 0) + { + //D[i] = 0; + j = j - M ; + //i = i - 1; + M = M / 2 ; + } + //D[i] = D[i] + 1; + j = j + M ; + i = i + 1 ; + + system("pause"); + } + */ + return 0; +} + +int FFTver3(double *x_r, double *x_i, double *y_r, double *y_i, int N) +{ + + // input : x = x_r + i * x_i + // output : y = y_r + i * y_i + int n,k; + for(n=0;n= M & M > 0) + { + j = j - M ; + M = M / 2 ; + } + j = j + M ; + i = i + 1 ; + } + n=1; + double u_r,u_i,w_r,w_i; + while(n < N) + { + u_r = cos(-2.0*M_PI/(2*n)); + u_i = sin(-2.0*M_PI/(2*n)); + w_r = 1; + w_i = 0; + for(i=0;i +#include +#include +#include +int SFT(double *x_r, double *x_i, double *y_r, double *y_i, int N); +int FFT_radix_2(double *x_r, double *x_i, double *y_r, double *y_i, int N); +int FFT_radix_3(double *x_r, double *x_i, double *y_r, double *y_i, int N); +int FFT_radix_5(double *x_r, double *x_i, double *y_r, double *y_i, int N); +int Initial(double *x, double *y, int N); +int Print_Complex_Vector(double *x, double *y, int N); + +int main() +{ + // y_k = sum(x_n * w^{-kn}, n=0..N-1) + // w = cos(2*pi/N)+isin(2*pi/N) + int k ,n ,N; + double *y_r, *y_i, *x_r, *x_i; + clock_t t1,t2; + printf("Please input N="); + scanf("%d", &N); + printf("N = %d\n", N); + + x_r = (double * ) malloc(N*sizeof(double)); + x_i = (double * ) malloc(N*sizeof(double)); + y_r = (double * ) malloc(N*sizeof(double)); + y_i = (double * ) malloc(N*sizeof(double)); + /* + Initial(x_r, x_i, N); + t1 = clock(); + SFT(x_r, x_i, y_r, y_i, N); + t2 = clock(); + printf("%f secs\n",1.0*(t2 -t1)/CLOCKS_PER_SEC); + //Print_Complex_Vector(y_r,y_i, N); + */ + Initial(x_r, x_i, N); + if((N%2) == 0) + { + t1 = clock(); + FFT_radix_2(x_r, x_i, y_r, y_i, N); + t2 = clock(); + printf("%f secs\n", 1.0*(t2-t1)/CLOCKS_PER_SEC); + system("pause"); + Print_Complex_Vector(y_r, y_i, N); + } + else if((N%3) == 0){ + t1 = clock(); + FFT_radix_3(x_r, x_i, y_r, y_i, N); + t2 = clock(); + printf("%f secs\n", 1.0*(t2-t1)/CLOCKS_PER_SEC); + system("pause"); + Print_Complex_Vector(y_r, y_i, N); + } + else if((N%5) == 0){ + t1 = clock(); + FFT_radix_5(x_r, x_i, y_r, y_i, N); + t2 = clock(); + printf("%f secs\n", 1.0*(t2-t1)/CLOCKS_PER_SEC); + system("pause"); + Print_Complex_Vector(y_r, y_i, N); + } + else { + printf("Error"); + return 0; + } + return 0; +} +int SFT(double *x_r, double *x_i, double *y_r, double *y_i, int N) +{ + int k, n; + double w_r, w_i; + for(k=0;k=0) printf("%d : %f +%f i\n", n, x[n], y[n]); + else printf("%d : %f %f i\n", n, x[n], y[n]); + } +} + + + + + + + + + + + + + + + + + + + + + + + + diff --git "a/\344\277\256\346\255\243\347\211\210mid.cpp" "b/\344\277\256\346\255\243\347\211\210mid.cpp" new file mode 100644 index 0000000..3c0257e --- /dev/null +++ "b/\344\277\256\346\255\243\347\211\210mid.cpp" @@ -0,0 +1,477 @@ +#include +#include +#include +#include +int SFT(double *x_r, double *x_i, double *y_r, double *y_i, int N); +int FFT_radix_2(double *x_r, double *x_i, double *y_r, double *y_i, int N); +int FFT_radix_3(double *x_r, double *x_i, double *y_r, double *y_i, int N); +int FFT_radix_5(double *x_r, double *x_i, double *y_r, double *y_i, int N); +int Initial(double *x, double *y, int N); +int Print_Complex_Vector(double *x, double *y, int N); + +int main() +{ + // y_k = sum(x_n * w^{-kn}, n=0..N-1) + // w = cos(2*pi/N)+isin(2*pi/N) + int k ,n ,N=1,q,p,r; + double *y_r, *y_i, *x_r, *x_i; + clock_t t1,t2; + printf("Please input p q r="); + scanf("%d %d %d", &p,&q,&r); + for(int i = 0;i

=0) printf("%d : %f +%f i\n", n, x[n], y[n]); + else printf("%d : %f %f i\n", n, x[n], y[n]); + } +} +