diff options
Diffstat (limited to 'NaiveFFT/main.c')
| -rw-r--r-- | NaiveFFT/main.c | 209 |
1 files changed, 190 insertions, 19 deletions
diff --git a/NaiveFFT/main.c b/NaiveFFT/main.c index fec893b..f8b56cd 100644 --- a/NaiveFFT/main.c +++ b/NaiveFFT/main.c @@ -11,28 +11,31 @@ #include <string.h> #include <math.h> #include <complex.h> +#include <time.h> //data array #define DATA_SIZE (8) -double data_i[DATA_SIZE] = {1.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f}; +//double data_i[DATA_SIZE] = {0.0f,1.0f,2.0f,3.0f,4.0f,5.0f,6.0f,7.0f}; +//double data_q[DATA_SIZE] = {0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f}; +double data_i[DATA_SIZE] = {1.0f,1.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f}; double data_q[DATA_SIZE] = {0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f,0.0f}; + //calc the fft void fft_if(double *x_i, double *x_q, int n, int inv) { int i=0,j=0,k=0,m=0, irem=0, sign; double tq,ti; - printf("FFT\n"); k=n/2; - for (j=1,i=1;i<n; i++) { + for (j=1,i=0;i<n,j<n; i++) { if (i<j) { //swap values - ti = x_i[j]; - tq = x_q[j]; + ti = x_i[j-1]; + tq = x_q[j-1]; - x_i[j] = x_i[i]; - x_q[j] = x_q[i]; + x_i[j-1] = x_i[i]; + x_q[j-1] = x_q[i]; x_i[i] = ti; x_q[i] = tq; @@ -53,7 +56,9 @@ void fft_if(double *x_i, double *x_q, int n, int inv) { } j = j+k; } - + //print rearranged buffer + /* + printf("FFT rearranged\n"); for (i=0;i<DATA_SIZE;i++) { printf("%.2f ",x_i[i]); } @@ -62,6 +67,7 @@ void fft_if(double *x_i, double *x_q, int n, int inv) { printf("%.2f ",x_q[i]); } printf("\n"); + */ //calculate number of stages: //m=log2(n) @@ -81,7 +87,7 @@ void fft_if(double *x_i, double *x_q, int n, int inv) { } // transform - for (i=0; i<m; i++) { + for (i=1; i<m; i++) { int le=pow(2,i); int le1=le/2; double ui=1.0; @@ -89,21 +95,24 @@ void fft_if(double *x_i, double *x_q, int n, int inv) { double wi = cos(M_PI/le1); double wq = sign*sin(M_PI/le1); - for (j=1; j<le1;j++) { - i=j; - while (i<n) { - int ip = i+le1; + printf("le1=%d\n",le1); + for (j=0; j<le1;j++) { + + k=j; + while (k<n) { + int ip = k+le1; ti = x_i[ip]*ui - x_q[ip]*uq; tq = x_q[ip]*ui + x_i[ip]*uq; - x_i[ip] = x_i[i] - ti; - x_q[ip] = x_q[i] - tq; + x_i[ip] = x_i[k-1] - ti; + x_q[ip] = x_q[k-1] - tq; - x_i[i] = x_i[i] + ti; - x_q[i] = x_q[i] + tq; + x_i[k-1] = x_i[k-1] + ti; + x_q[k-1] = x_q[k-1] + tq; - i = i+le; + k = k+le; + printf("k=%d\n",k-1); } double temp = ui*wi - uq*wq; uq = uq*wi + ui*wq; @@ -121,6 +130,124 @@ void fft_if(double *x_i, double *x_q, int n, int inv) { } +#define complex_mul_re(a_re, a_im, b_re, b_im) (a_re * b_re - a_im * b_im) +#define complex_mul_im(a_re, a_im, b_re, b_im) (a_re * b_im + a_im * b_re) +//https://github.com/rshuston/FFT-C/blob/master/libfft/fft.c +void ffti_shuffle_1(double *x_i, double *x_q, uint64_t n) { + int Nd2 = n>>1; + int Nm1 = n-1; + int i,j; + + + for (i = 0, j = 0; i < n; i++) { + if (j > i) { + double tmp_r = x_i[i]; + double tmp_i = x_q[i]; + //data[i] = data[j]; + x_i[i] = x_q[j]; + x_q[i] = x_q[j]; + //data[j] = tmp; + x_i[j] = tmp_r; + x_q[j] = tmp_i; + } + + /* + * Find least significant zero bit + */ + + unsigned lszb = ~i & (i + 1); + + /* + * Use division to bit-reverse the single bit so that we now have + * the most significant zero bit + * + * N = 2^r = 2^(m+1) + * Nd2 = N/2 = 2^m + * if lszb = 2^k, where k is within the range of 0...m, then + * mszb = Nd2 / lszb + * = 2^m / 2^k + * = 2^(m-k) + * = bit-reversed value of lszb + */ + + unsigned mszb = Nd2 / lszb; + + /* + * Toggle bits with bit-reverse mask + */ + + unsigned bits = Nm1 & ~(mszb - 1); + j ^= bits; + } +} +void fft_1(double *x_i, double *x_q, uint64_t n, uint64_t inv) { + uint64_t n_log2; + uint64_t r; + uint64_t m, md2; + uint64_t i,j,k; + uint64_t i_e, i_o; + double theta_2pi; + double theta; + + double Wm_r, Wm_i, Wmk_r, Wmk_i; + double u_r, u_i, t_r, t_i; + + //find log of n + i=n; + n_log2 = 0; + while (i>1) { + i=i/2; + n_log2+=1; + } + + if (inv==1) { + theta_2pi = -2*M_PI; + } else { + theta_2pi = 2*M_PI; + } + + + for (i=1; i<= n_log2; i++) { + m = 1 << i; + md2 = m >> 1; + theta = theta_2pi / m; + Wm_r = cos(theta); + Wm_i = sin(theta); + + + for (j=0; j<n; j+=m) { + Wmk_r = 1.0f; + Wmk_i = 0.0f; + + for (k=0; k<md2; k++) { + i_e = j+k; + i_o = i_e + md2; + + u_r = x_i[i_e]; + u_i = x_q[i_e]; + + //t_r = Wmk_r * x_i[i_o] - Wmk_i * x_q[i_o]; + //t_i = Wmk_r * x_q[i_o] - Wmk_i * x_i[i_o]; + t_r = complex_mul_re(Wmk_r, Wmk_i, x_i[i_o], x_q[i_o]); + t_i = complex_mul_im(Wmk_r, Wmk_i, x_i[i_o], x_q[i_o]); + + + x_i[i_e] = u_r + t_r; + x_q[i_e] = u_i + t_i; + + x_i[i_o] = u_r - t_r; + x_q[i_o] = u_i - t_i; + + t_r = complex_mul_re(Wmk_r, Wmk_i, Wm_r, Wm_i); + t_i = complex_mul_im(Wmk_r, Wmk_i, Wm_r, Wm_i); + + Wmk_r = t_r; + Wmk_i = t_i; + } + } + } +} + //dft works fine void dft(double *x_i, double *x_q, int n, int inv) { double Wn,Wk; @@ -181,10 +308,13 @@ int main(int argc, const char * argv[]) { printf("%.2f ",data_q[i]); } printf("\n"); - fft_if(data_i, data_q, DATA_SIZE, 0); + ffti_shuffle_1(data_i, data_q, DATA_SIZE); + fft_1(data_i, data_q, DATA_SIZE, 0); //dft(data_i, data_q, DATA_SIZE,0); + printf("after FFT\n"); + for (i=0;i<DATA_SIZE;i++) { printf("%.2f ",data_i[i]); } @@ -194,5 +324,46 @@ int main(int argc, const char * argv[]) { } printf("\n"); + //big buffer test + { + time_t timer; + char buffer[32]; + struct tm *tm_info; + + timer = time(NULL); + tm_info = localtime(&timer); + + strftime(buffer, 32, "%H:%M:%S", tm_info); + puts(buffer); + + } + + //test big big fft + { + double *d_i=malloc(sizeof(double)*4*8192+16); + double *d_q=malloc(sizeof(double)*4*8192+16); + for (i=0;i<4*8192;i++) { + d_i[i] = sin(i*1.0/128); + d_q[i] = 0.0f; + } + + dft(d_i, d_q, 4*8192,0); + free(d_i); d_i = NULL; + free(d_q); d_q = NULL; + } + + + { + time_t timer; + char buffer[32]; + struct tm *tm_info; + + timer = time(NULL); + tm_info = localtime(&timer); + + strftime(buffer, 32, "%H:%M:%S", tm_info); + puts(buffer); + + } return 0; } |