diff options
| -rw-r--r-- | md/writeup.md | 1 | ||||
| -rw-r--r-- | md/writeup/naive_fft_implementation_in_c.md | 241 |
2 files changed, 242 insertions, 0 deletions
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 + +<script> +MathJax = { + tex: { + inlineMath: [['$', '$'], ['\\(', '\\)']] + } +}; +</script> + +<script type="text/javascript" id="MathJax-script" async + src="/js/tex-chtml.js"> +</script> + +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<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]; + } + +} +``` +### FFT +#### Shuffling arrays + +```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; + } +} +``` + + +#### 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<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 = 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; + } + } + } +``` + + + +## Octave verification code + +```matlab +data1 = [1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0] +res1 = fft(data1) + +data2 = [1.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0] +res2 = fft(data2) + +data3 = [1.0 1.0 i i 0.0 0.0 0.0 0.0] +res3 = fft(data3) + +data4 = [1.0+i 1.0+i 1.0+i 1.0+i 0.0 0.0 0.0 0.0] +res4 = fft(data4) + +``` + +### Source + +Source located in main branch of git repo + +``` +git clone http://git.main.lv/cgit.cgi/NaiveFFT.git +git checkout main +``` + +### Build + +#### Linux + +``` +cd NaiveFFT/Build +make +``` + +#### Macos + +Open with Xcode and build + + +## Links + + +[01] https://rosettacode.org/wiki/Fast_Fourier_transform#C +[02] https://www.math.wustl.edu/~victor/mfmm/fourier/fft.c +[03] https://www.strauss-engineering.ch/libdsp.html +[04] http://www.iowahills.com/FFTCode.html +[05] https://github.com/rshuston/FFT-C/blob/master/libfft/fft.c +[06] http://www.guitarscience.net/papers/fftalg.pdf +[07] https://root.cern/doc/master/FFT_8C_source.html +[08] https://lloydrochester.com/post/c/example-fft/ +[09] https://github.com/mborgerding/kissfft/blob/master/kiss_fft.c +[10] https://community.vcvrack.com/t/complete-list-of-c-c-fft-libraries/9153 +[11] https://octave.org/doc/v4.2.1/Signal-Processing.html +[12] https://en.wikipedia.org/wiki/Fast_Fourier_transform +[13] http://www.iowahills.com/FFTCode.html +[14] Digital Signal Processing: A Practical Approach by (Emmanuel Ifeachor, Barrie Jervis) + + |