summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--NaiveFFT/fft.c271
-rw-r--r--NaiveFFT/fft.h20
-rw-r--r--NaiveFFT/main.c276
3 files changed, 294 insertions, 273 deletions
diff --git a/NaiveFFT/fft.c b/NaiveFFT/fft.c
new file mode 100644
index 0000000..cfcf3fa
--- /dev/null
+++ b/NaiveFFT/fft.c
@@ -0,0 +1,271 @@
+//
+// fft.c
+// NaiveFFT
+//
+// Created by Jacky Jack on 26/08/2021.
+//
+
+#include "fft.h"
+
+
+
+
+//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;
+
+ k=n/2;
+ for (j=1,i=0;i<n,j<n; i++) {
+ if (i<j) {
+
+ //swap values
+ ti = x_i[j-1];
+ tq = x_q[j-1];
+
+ x_i[j-1] = x_i[i];
+ x_q[j-1] = x_q[i];
+
+ x_i[i] = ti;
+ x_q[i] = tq;
+ printf("i=%d j=%d\n",i,j);
+ //find middle point
+ k = n/2;
+ while (k<j) {
+ j = j-k;
+ k=k/2;
+ }
+ } else {
+ printf("i=%d j=%d\n",i,j);
+ k = n/2;
+ while (k<j) {
+ j = j - k;
+ k = k/2;
+ }
+ }
+ j = j+k;
+ }
+ //calculate number of stages:
+ //m=log2(n)
+ m=0;
+ irem = n;
+ while(irem>1) {
+ irem = irem/2;
+ m = m+1;
+ }
+
+ //FFT or IFFT
+ if (inv==1) {
+ sign = 1;
+ } else {
+ sign = -1;
+ }
+
+ // transform
+ for (i=1; i<m; i++) {
+ int le=pow(2,i);
+ int le1=le/2;
+ double ui=1.0;
+ double uq=0.0;
+ double wi = cos(M_PI/le1);
+ double wq = sign*sin(M_PI/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[k-1] - ti;
+ x_q[ip] = x_q[k-1] - tq;
+
+ x_i[k-1] = x_i[k-1] + ti;
+ x_q[k-1] = x_q[k-1] + tq;
+
+ k = k+le;
+ printf("k=%d\n",k-1);
+ }
+ double temp = ui*wi - uq*wq;
+ uq = uq*wi + ui*wq;
+ ui = temp;
+ }
+ }
+
+ //if inverse
+ if (inv==1) {
+ for (i=0;i<n;i++) {
+ x_i[i] = x_i[i]/n;
+ x_q[i] = x_q[i]/n;
+ }
+ }
+
+}
+
+#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;
+ //static array
+ //double Xi[DATA_SIZE],Xq[DATA_SIZE];
+ //dynamic array
+ double *Xi, *Xq;
+ double c,s;
+ int i,j;
+
+ Xi = malloc(n*sizeof(double));
+ Xq = malloc(n*sizeof(double));
+
+ Wn = 2*M_PI/n;
+
+ if (inv==1) {
+ Wn=-Wn;
+ }
+
+ for (i=0;i<n;i++) {
+ Xi[i] = 0.0f;
+ Xq[i] = 0.0f;
+
+ Wk = i*Wn;
+ for (j=0;j<n;j++) {
+ c = cos(j*Wk);
+ s = sin(j*Wk);
+
+ //i - real, q - imaginary
+ Xi[i] = Xi[i] + x_i[j]*c + x_q[j]*s;
+ Xq[i] = Xq[i] - x_i[j]*s + x_q[j]*c;
+ }
+
+ if (inv==1) {
+ Xi[i] = Xi[i]/n;
+ Xq[i] = Xq[i]/n;
+ }
+ }
+
+ for (i=0;i<n;i++) {
+ x_i[i] = Xi[i];
+ x_q[i] = Xq[i];
+ }
+
+}
diff --git a/NaiveFFT/fft.h b/NaiveFFT/fft.h
new file mode 100644
index 0000000..f1ee783
--- /dev/null
+++ b/NaiveFFT/fft.h
@@ -0,0 +1,20 @@
+//
+// fft.h
+// NaiveFFT
+//
+// Created by Jacky Jack on 26/08/2021.
+//
+
+#ifndef fft_h
+#define fft_h
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+void fft_if(double *x_i, double *x_q, int n, int inv);
+void ffti_shuffle_1(double *x_i, double *x_q, uint64_t n);
+void fft_1(double *x_i, double *x_q, uint64_t n, uint64_t inv);
+void dft(double *x_i, double *x_q, int n, int inv);
+
+#endif /* fft_h */
diff --git a/NaiveFFT/main.c b/NaiveFFT/main.c
index 1d71ac5..718e5d0 100644
--- a/NaiveFFT/main.c
+++ b/NaiveFFT/main.c
@@ -14,286 +14,16 @@
#include <time.h>
#include <stdint.h>
-//data array
+#include "fft.h"
+
+//test data array
#define DATA_SIZE (8)
//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;
-
- k=n/2;
- for (j=1,i=0;i<n,j<n; i++) {
- if (i<j) {
-
- //swap values
- ti = x_i[j-1];
- tq = x_q[j-1];
-
- x_i[j-1] = x_i[i];
- x_q[j-1] = x_q[i];
-
- x_i[i] = ti;
- x_q[i] = tq;
- printf("i=%d j=%d\n",i,j);
- //find middle point
- k = n/2;
- while (k<j) {
- j = j-k;
- k=k/2;
- }
- } else {
- printf("i=%d j=%d\n",i,j);
- k = n/2;
- while (k<j) {
- j = j - k;
- k = k/2;
- }
- }
- j = j+k;
- }
- //print rearranged buffer
- /*
- printf("FFT rearranged\n");
- for (i=0;i<DATA_SIZE;i++) {
- printf("%.2f ",x_i[i]);
- }
- printf("\n");
- for (i=0;i<DATA_SIZE;i++) {
- printf("%.2f ",x_q[i]);
- }
- printf("\n");
- */
-
- //calculate number of stages:
- //m=log2(n)
- m=0;
- irem = n;
- while(irem>1) {
- irem = irem/2;
- m = m+1;
- }
-
-
- //FFT or IFFT
- if (inv==1) {
- sign = 1;
- } else {
- sign = -1;
- }
-
- // transform
- for (i=1; i<m; i++) {
- int le=pow(2,i);
- int le1=le/2;
- double ui=1.0;
- double uq=0.0;
- double wi = cos(M_PI/le1);
- double wq = sign*sin(M_PI/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[k-1] - ti;
- x_q[ip] = x_q[k-1] - tq;
-
- x_i[k-1] = x_i[k-1] + ti;
- x_q[k-1] = x_q[k-1] + tq;
-
- k = k+le;
- printf("k=%d\n",k-1);
- }
- double temp = ui*wi - uq*wq;
- uq = uq*wi + ui*wq;
- ui = temp;
- }
- }
-
- //if inverse
- if (inv==1) {
- for (i=0;i<n;i++) {
- x_i[i] = x_i[i]/n;
- x_q[i] = x_q[i]/n;
- }
- }
-
-}
-
-#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;
- //static array
- //double Xi[DATA_SIZE],Xq[DATA_SIZE];
- //dynamic array
- double *Xi, *Xq;
- double c,s;
- int i,j;
-
- Xi = malloc(n*sizeof(double));
- Xq = malloc(n*sizeof(double));
-
- Wn = 2*M_PI/n;
-
- if (inv==1) {
- Wn=-Wn;
- }
-
- for (i=0;i<n;i++) {
- Xi[i] = 0.0f;
- Xq[i] = 0.0f;
-
- Wk = i*Wn;
- for (j=0;j<n;j++) {
- c = cos(j*Wk);
- s = sin(j*Wk);
-
- //i - real, q - imaginary
- Xi[i] = Xi[i] + x_i[j]*c + x_q[j]*s;
- Xq[i] = Xq[i] - x_i[j]*s + x_q[j]*c;
- }
-
- if (inv==1) {
- Xi[i] = Xi[i]/n;
- Xq[i] = Xq[i]/n;
- }
- }
-
- for (i=0;i<n;i++) {
- x_i[i] = Xi[i];
- x_q[i] = Xq[i];
- }
-
-}
int main(int argc, const char * argv[]) {
int i;
generated by cgit v1.2.3 (git 2.39.1) at 2024-12-25 17:38:53 +0000