From 6f5d2facf25fe75396f31b46bbf7e759845bacd2 Mon Sep 17 00:00:00 2001 From: FreeArtMan Date: Wed, 25 Aug 2021 09:30:47 +0100 Subject: Added FFT implementation post --- md/writeup.md | 1 + md/writeup/naive_fft_implementation_in_c.md | 241 ++++++++++++++++++++++++++++ 2 files changed, 242 insertions(+) create mode 100644 md/writeup/naive_fft_implementation_in_c.md diff --git a/md/writeup.md b/md/writeup.md index 9987142..846bfbf 100644 --- a/md/writeup.md +++ b/md/writeup.md @@ -4,6 +4,7 @@ title: Writeup page ## Writeup's +[Naive FFT implementation in C](writeup/naive_fft_implementation_in_c.md) [Web assembly audio with fir filter](writeup/web_assembly_audio_with_fir_filter.md) [Calculate fir coefficients with C](writeup/calculate_fir_coefficients_with_c.md) [ARM64 assembly crc32](writeup/arm64_assembly_crc32.md) diff --git a/md/writeup/naive_fft_implementation_in_c.md b/md/writeup/naive_fft_implementation_in_c.md new file mode 100644 index 0000000..a6ca9dc --- /dev/null +++ b/md/writeup/naive_fft_implementation_in_c.md @@ -0,0 +1,241 @@ +title:Naive FFT implementation in C +keywords:c,linux,fir,dsp,octave,fft + +# Naive FFT implementation in C + +## Intro + + + + + +Time to implement DFT/FFT by myself. Just followed Ifeachor book on theory part and +IOWA FFT source on implementation side. + + +## Implementation + +### DFT + + +$$X(k) = F_{D}[x(nT)] = $$ + +$$ \sum_{n=0}^{N-1} x(nT) e^{-jk\omega nT} , k=0,1,...,N-1 $$ + + +```c +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>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; + } +} +``` + + +#### Implementation + +```c +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