Created
May 12, 2022 12:43
-
-
Save amb/4f70bcbea897024452d683b40d18be1f to your computer and use it in GitHub Desktop.
Simple Nim FFT
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| # OTFFT library | |
| # http://wwwa.pikara.ne.jp/okojisan/otfft-en/optimization1.html | |
| # This is +20-50% improvement | |
| const thetaLutSize = 2048 | |
| const thetaLut = static: | |
| var arr: array[thetaLutSize, Complex[float]] | |
| let step = 2.0*PI/float(thetaLutSize) | |
| for k, v in mpairs(arr): | |
| v = complex(cos(step * float(k)), -sin(step * float(k))) | |
| arr | |
| proc fft0(n: int, s: int, eo: bool, x: var seq[Complex[float]], y: var seq[Complex[float]]) = | |
| ## Fast Fourier Transform | |
| ## | |
| ## Inputs: | |
| ## - `n` Sequence length. **Must be power of two** | |
| ## - `s` Stride | |
| ## - `eo` x is output if eo == 0, y is output if eo == 1 | |
| ## - `x` Input sequence(or output sequence if eo == 0) | |
| ## - `y` Work area(or output sequence if eo == 1) | |
| ## | |
| ## Returns: | |
| ## - Output sequence, either `x` or `y` | |
| let m: int = n div 2 | |
| let theta0: float = 2.0*PI/float(n) | |
| # let theta0: float = float(thetaLutSize)/float(n) | |
| if n == 1: | |
| if eo: | |
| for q in 0..<s: | |
| y[q] = x[q] | |
| else: | |
| for p in 0..<m: | |
| let fp = float(p)*theta0 | |
| let wp = complex(cos(fp), -sin(fp)) | |
| # let fp = int(float(p)*theta0) | |
| # let wp = thetaLut[fp] | |
| for q in 0..<s: | |
| let a = x[q + s*(p+0)] | |
| let b = x[q + s*(p+m)] | |
| y[q + s*(2*p+0)] = a + b | |
| y[q + s*(2*p+1)] = (a - b) * wp | |
| fft0(n div 2, 2 * s, not eo, y, x) | |
| proc fft*(x: var seq[Complex[float]]) = | |
| # n : sequence length | |
| # x : input/output sequence | |
| # Input length has to be a power of two | |
| assert x.len > 0 | |
| assert x.len.isPowerOfTwo() | |
| var y: seq[Complex[float]] = newSeq[Complex[float]](x.len) | |
| fft0(x.len, 1, false, x, y) | |
| proc ifft*(x: var seq[Complex[float]]) = | |
| # n : sequence length | |
| # x : input/output sequence | |
| var n: int = x.len | |
| let fn = complex(1.0/float(n)) | |
| for p in 0..<n: | |
| x[p] = (x[p]*fn).conjugate | |
| var y: seq[Complex[float]] = newSeq[Complex[float]](n) | |
| fft0(n, 1, false, x, y) | |
| for p in 0..<n: | |
| x[p] = x[p].conjugate |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment