--- title: "FFT algorithm and program" subject: "Computational Techniques in Civil Engineering" chapter: 2 key: "semester_8/computational_techniques_in_civil_engineering:fft-algorithm-program" source: https://www.bctengineer.com/q/semester_8/computational_techniques_in_civil_engineering/fft-algorithm-program site: BCT Engineering --- # FFT algorithm and program Computational Techniques in Civil Engineering · Chapter 2: Solutions of linear equations - Asked in: 2070 Magh, 12 marks - Video: [DIT FFT algorithm | Butterfly diagram | Digital signal processing](https://www.youtube.com/watch?v=RioJKiSBlyg) — Smart Engineer - Answer checked by Claude Opus 5.5 (high effort) on 2026-10-14 - Page: https://www.bctengineer.com/q/semester_8/computational_techniques_in_civil_engineering/fft-algorithm-program - Whole chapter: https://www.bctengineer.com/md/semester_8/computational_techniques_in_civil_engineering/chapter-2 ## Question Write an algorithm and a program (C or Fortran or Matlab) for the fast Fourier transform. With a suitable example explain what parameters can be identified with the help of time domain and frequency domain. ## Short answer The radix-2 FFT reorders samples in bit-reversed order and applies butterfly stages, computing the DFT in O(N log N) operations instead of N^2. For x(t)=2 sin(2 pi t)+sin(6 pi t) sampled at 8 Hz, the time domain shows a mixed waveform, while the frequency domain identifies peaks at 1 Hz amplitude 2 and 3 Hz amplitude 1; the Nyquist frequency is 4 Hz. ## Answer #### Fast Fourier transform: idea The DFT of $N$ samples is $X_k=\sum_{n=0}^{N-1}x_n W_N^{kn}$ with $W_N=e^{-i2\pi/N}$, needing $N^2$ operations. The Cooley-Tukey radix-2 FFT (for $N=2^m$) splits the series into even and odd samples: $$ X_k=E_k+W_N^kO_k,\qquad X_{k+N/2}=E_k-W_N^kO_k,\quad k=0,\dots,\tfrac N2-1 $$ and repeats the split $\log_2N$ times. The cost is $\tfrac N2\log_2N$ complex multiplications. #### Algorithm 1. Read $N$ (a power of 2) and the samples $x_n$ at interval $\Delta t$. 2. Reorder the samples in **bit-reversed** order. 3. For stage length $len=2,4,8,\dots,N$: compute $W=e^{-i2\pi/len}$; for each block, for each $k=0\dots len/2-1$ form $u=x_{i+k}$, $v=x_{i+k+len/2}\,W^k$, then set $x_{i+k}=u+v$, $x_{i+k+len/2}=u-v$. 4. After the last stage $x$ holds $X_k$. Amplitude $=2|X_k|/N$, frequency $f_k=k/(N\Delta t)$. ```text x0 --o----o----o-- X0 x4 --o----| | x2 --o----o | (butterflies: x6 --o----| | a+Wb , a-Wb) ... ``` #### Program (C) ```c #include #include #include #define PI 3.14159265358979323846 /* in-place radix-2 FFT, n must be a power of 2 */ void fft(double complex x[], int n) { for (int i = 1, j = 0; i < n; i++) { /* bit reversal */ int bit = n >> 1; for (; j & bit; bit >>= 1) j ^= bit; j ^= bit; if (i < j) { double complex t = x[i]; x[i] = x[j]; x[j] = t; } } for (int len = 2; len <= n; len <<= 1) { /* butterfly stages */ double complex wlen = cexp(-2.0 * I * PI / len); for (int i = 0; i < n; i += len) { double complex w = 1.0; for (int k = 0; k < len / 2; k++) { double complex u = x[i + k]; double complex v = x[i + k + len / 2] * w; x[i + k] = u + v; x[i + k + len / 2] = u - v; w *= wlen; } } } } int main(void) { int n = 8; double fs = 8.0; double complex x[8]; for (int i = 0; i < n; i++) { /* 1 Hz (amp 2) + 3 Hz (amp 1) */ double t = i / fs; x[i] = 2*sin(2*PI*1*t) + 1*sin(2*PI*3*t); } fft(x, n); for (int k = 0; k < n; k++) printf("f=%.1f Hz amp=%.3f\n", k*fs/n, 2*cabs(x[k])/n); return 0; } ``` #### Example: time domain and frequency domain Signal $x(t)=2\sin(2\pi t)+\sin(6\pi t)$ sampled at 8 Hz for 1 s ($N=8$). The program output is: | $f$ (Hz) | 0 | 1 | 2 | 3 | 4 | |---|---|---|---|---|---| | amplitude | 0 | 2.000 | 0 | 1.000 | 0 | (components at $k=5,6,7$ are the mirror image of $k=3,2,1$.) - **Time domain** shows the amplitude of vibration against time, peak value, duration, arrival time and decay of the signal, but the individual components are mixed together. - **Frequency domain** shows the dominant frequencies (1 Hz and 3 Hz), their amplitudes (2 and 1) and phase. For a structure, the peak frequency is the natural frequency, so resonance with earthquake or machine frequency can be checked; noise is separated by filtering. - Frequency resolution is $\Delta f=1/(N\Delta t)=1$ Hz; the highest frequency identified is the Nyquist frequency $f_s/2=4$ Hz. ## Check yourself The answer's key points as multiple-choice questions, with the correct option and why. 1. What is the advantage of the radix-2 FFT over the direct DFT for N samples? A. N^2 operations become N/2 log2 N complex multiplications B. N^2 operations become N^3 C. It needs no complex numbers D. It works for any N, not only powers of 2 Answer: A. The DFT needs N^2 operations; splitting into even and odd parts log2 N times reduces this to (N/2) log2 N multiplications. Radix-2 needs N = 2^m. (Key point: FFT cost) 2. In each FFT butterfly with E_k (even part) and O_k (odd part), what are X_k and X_{k+N/2}? A. X_k = E_k - W^k O_k, X_{k+N/2} = E_k - W^k O_k B. X_k = E_k + O_k, X_{k+N/2} = E_k - O_k C. X_k = E_k + W^k O_k, X_{k+N/2} = E_k - W^k O_k D. X_k = W^k E_k, X_{k+N/2} = O_k Answer: C. The butterfly gives X_k = E_k + W_N^k O_k and X_{k+N/2} = E_k - W_N^k O_k, so one product yields two outputs. (Key point: Butterfly operation) 3. For N = 8 samples at 8 Hz for 1 s, what are the frequency resolution and the highest identifiable frequency? A. 1 Hz and 8 Hz B. 1 Hz and the Nyquist frequency 4 Hz C. 0.5 Hz and 4 Hz D. 8 Hz and 1 Hz Answer: B. Delta f = 1/(N dt) = 1 Hz, and the highest frequency is the Nyquist frequency fs/2 = 4 Hz. Components k = 5, 6, 7 mirror 3, 2, 1. (Key point: Frequency resolution) 4. What does the frequency domain reveal that the time domain does not show clearly? A. The arrival time of the signal B. The dominant frequencies (1 Hz and 3 Hz) and amplitudes (2 and 1), e.g. the natural frequency for resonance checks C. The duration of the signal D. The decay of the signal with time Answer: B. Time domain shows peak, duration, arrival and decay with mixed components; the frequency domain separates components and shows dominant frequencies and amplitudes. (Key point: Time vs frequency domain) ## Related in Computational Techniques in Civil Engineering - [Discrete and Fast Fourier Transform](https://www.bctengineer.com/question.md?key=semester_8%2Fcomputational_techniques_in_civil_engineering%3Adft-fft-brief) — Computational Techniques in Civil Engineering, chapter 2 ## Related in other subjects - [Hilbert transform in frequency domain](https://www.bctengineer.com/question.md?key=semester_6%2Fcommunication_systems%3Ahilbert-define-frequency-domain) — Communication Systems, chapter 2 - [Single tone DSB-FC AM in time and frequency domain](https://www.bctengineer.com/question.md?key=semester_6%2Fcommunication_systems%3Asingle-tone-dsbfc-time-freq) — Communication Systems, chapter 2 - [Angle modulation; single tone FM time and frequency domain](https://www.bctengineer.com/question.md?key=semester_6%2Fcommunication_systems%3Aangle-mod-single-tone-fm-time-freq) — Communication Systems, chapter 3 - [How fast is FFT; X(7) by DIT-FFT](https://www.bctengineer.com/question.md?key=semester_7%2Fdigital_signal_analysis_and_processing%3Afft-speed-x7) — Digital Signal Analysis and Processing, chapter 7, asked 2×