From 505fcaac2ef3e1892dc6badbd0fe4d49e0525cf2 Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Sun, 13 Oct 2019 15:43:49 -0400 Subject: [PATCH 1/5] Start fourier transform implementations in Nim --- contents/cooley_tukey/code/nim/.gitignore | 1 + contents/cooley_tukey/code/nim/fft.nim | 93 +++++++++++++++++++++++ 2 files changed, 94 insertions(+) create mode 100644 contents/cooley_tukey/code/nim/.gitignore create mode 100644 contents/cooley_tukey/code/nim/fft.nim diff --git a/contents/cooley_tukey/code/nim/.gitignore b/contents/cooley_tukey/code/nim/.gitignore new file mode 100644 index 000000000..e91bc4f92 --- /dev/null +++ b/contents/cooley_tukey/code/nim/.gitignore @@ -0,0 +1 @@ +fft diff --git a/contents/cooley_tukey/code/nim/fft.nim b/contents/cooley_tukey/code/nim/fft.nim new file mode 100644 index 000000000..b1fc5805d --- /dev/null +++ b/contents/cooley_tukey/code/nim/fft.nim @@ -0,0 +1,93 @@ +from math import PI +from random import rand, randomize +from sequtils import map, toSeq, zip +from sugar import `=>` +import complex +# $ nimble install fftw3 +import fftw3 + +# For some reason this isn't in the Nim fftw3 bindings. +const FFTW_FORWARD = -1 + +proc allClose[T](reference: openArray[T], calculated: openArray[T], thresh = 1.0e-13): bool = + if reference.len != calculated.len: + return false + result = true + for (r, c) in zip(reference, calculated): + if abs(r - c) > thresh: + result = false + break + +proc skip[T](inp: openArray[T], stride: int): seq[T] = + ## Like `inp[::stride]` in Python. + if stride < 1: + return toSeq(inp) + result = @[] + var i = 0 + while i < inp.len: + result.add(inp[i]) + i += stride + +proc nim_to_fftw(inp: Complex64): fftw_complex = + [inp.re, inp.im] + +proc fftw_to_nim(inp: fftw_complex): Complex64 = + complex64(inp[0], inp[1]) + +proc nim_to_fftw(inp: openArray[Complex64]): seq[fftw_complex] = + inp.map(v => nim_to_fftw(v)) + +proc fftw_to_nim(inp: openArray[fftw_complex]): seq[Complex64] = + inp.map(v => fftw_to_nim(v)) + +proc dft_fftw3(inp: openArray[Complex64]): seq[Complex64] = + var + fftw_inp = nim_to_fftw(inp) + fftw_out = newSeq[fftw_complex](inp.len) + let p = fftw_plan_dft_1d(cint(inp.len), + addr(fftw_inp[low(fftw_inp)]), + addr(fftw_out[low(fftw_out)]), + FFTW_FORWARD, + FFTW_ESTIMATE) + fftw_execute(p) + fftw_destroy_plan(p) + fftw_to_nim(fftw_out) + +proc dft[T: SomeFloat](x: openArray[Complex[T]]): seq[Complex[T]] = + let n = x.len + result = newSeq[Complex[T]](n) + for i in 0..n - 1: + for k in 0..n - 1: + result[i] += x[k] * exp(-complex(0.0, 2.0) * PI * float(i * k / n)) + +proc cooley_tukey[T: SomeFloat](x: openArray[T]): seq[Complex[T]] = + let n = x.len + # if n <= 1: + # return x.map(m => complex(m, m)) + result = newSeq[Complex[T]](n) + let + even = cooley_tukey(skip(x, 2)) + odd = cooley_tukey(skip(x[1..high(x)], 2)) + midpoint = n div 2 + for k in 0..midpoint - 1: + let exp_term = exp(-complex(0.0, 2.0) * PI * float(k / n)) * odd[k] + result[k] = even[k] + exp_term + result[k + midpoint] = even[k] - exp_term + +# proc bit_reverse[T: SomeNumber](x: openArray[T]): seq[T] + +# proc iterative_cooley_tukey[T: SomeFloat](x: openArray[T]): seq[Complex[T]] = +# let n = x.len +# result = bit_reverse(x) + +when isMainModule: + randomize() + let + x = toSeq(1..64).map(i => complex64(rand(1.0), 0.0)) + # y = cooley_tukey(x) + # z = iterative_cooley_tukey(x) + t = dft(x) + reference = dft_fftw3(x) + # assert(allClose(reference, y)) + # assert(allClose(reference, z)) + assert(allClose(reference, t)) From 302653bcc275e4f1630a60fb3f61d0429a59102c Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Sun, 13 Oct 2019 17:52:48 -0400 Subject: [PATCH 2/5] Implement Cooley-Tukey algorithm in Nim --- contents/cooley_tukey/code/nim/fft.nim | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/contents/cooley_tukey/code/nim/fft.nim b/contents/cooley_tukey/code/nim/fft.nim index b1fc5805d..33cd067f4 100644 --- a/contents/cooley_tukey/code/nim/fft.nim +++ b/contents/cooley_tukey/code/nim/fft.nim @@ -9,7 +9,7 @@ import fftw3 # For some reason this isn't in the Nim fftw3 bindings. const FFTW_FORWARD = -1 -proc allClose[T](reference: openArray[T], calculated: openArray[T], thresh = 1.0e-13): bool = +proc allClose[T: SomeNumber | Complex[SomeFloat]](reference: openArray[T], calculated: openArray[T], thresh = 1.0e-13): bool = if reference.len != calculated.len: return false result = true @@ -60,10 +60,10 @@ proc dft[T: SomeFloat](x: openArray[Complex[T]]): seq[Complex[T]] = for k in 0..n - 1: result[i] += x[k] * exp(-complex(0.0, 2.0) * PI * float(i * k / n)) -proc cooley_tukey[T: SomeFloat](x: openArray[T]): seq[Complex[T]] = +proc cooley_tukey[T: SomeFloat](x: openArray[Complex[T]]): seq[Complex[T]] = let n = x.len - # if n <= 1: - # return x.map(m => complex(m, m)) + if n <= 1: + return toSeq(x) result = newSeq[Complex[T]](n) let even = cooley_tukey(skip(x, 2)) @@ -84,10 +84,10 @@ when isMainModule: randomize() let x = toSeq(1..64).map(i => complex64(rand(1.0), 0.0)) - # y = cooley_tukey(x) + y = cooley_tukey(x) # z = iterative_cooley_tukey(x) t = dft(x) reference = dft_fftw3(x) - # assert(allClose(reference, y)) + assert(allClose(reference, y)) # assert(allClose(reference, z)) assert(allClose(reference, t)) From 7907df00d614bc5ce824f2e8d5d69d9db06e1cb9 Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Sun, 13 Oct 2019 19:13:50 -0400 Subject: [PATCH 3/5] Implement iterative Cooley-Tukey in Nim --- contents/cooley_tukey/code/nim/fft.nim | 43 +++++++++++++++++++++----- 1 file changed, 35 insertions(+), 8 deletions(-) diff --git a/contents/cooley_tukey/code/nim/fft.nim b/contents/cooley_tukey/code/nim/fft.nim index 33cd067f4..bd7ad4266 100644 --- a/contents/cooley_tukey/code/nim/fft.nim +++ b/contents/cooley_tukey/code/nim/fft.nim @@ -1,4 +1,5 @@ -from math import PI +from bitops import fastLog2 +from math import PI, `^` from random import rand, randomize from sequtils import map, toSeq, zip from sugar import `=>` @@ -9,7 +10,9 @@ import fftw3 # For some reason this isn't in the Nim fftw3 bindings. const FFTW_FORWARD = -1 -proc allClose[T: SomeNumber | Complex[SomeFloat]](reference: openArray[T], calculated: openArray[T], thresh = 1.0e-13): bool = +type InpNumber = SomeNumber | Complex[SomeFloat] + +proc allClose[T: InpNumber](reference: openArray[T], calculated: openArray[T], thresh = 1.0e-11): bool = if reference.len != calculated.len: return false result = true @@ -74,20 +77,44 @@ proc cooley_tukey[T: SomeFloat](x: openArray[Complex[T]]): seq[Complex[T]] = result[k] = even[k] + exp_term result[k + midpoint] = even[k] - exp_term -# proc bit_reverse[T: SomeNumber](x: openArray[T]): seq[T] +proc bit_reverse[T: InpNumber](x: openArray[T]): seq[T] = + let + n = x.len + l2 = fastLog2(n) + result = newSeq[T](n) + for k in 0..n - 1: + let s = toSeq(0..l2 - 1) + var b = 0 + for i in s: + if (k shr i and 1) == 1: + b += 1 shl (l2 - 1 - i) + result[k] = x[b] + result[b] = x[k] -# proc iterative_cooley_tukey[T: SomeFloat](x: openArray[T]): seq[Complex[T]] = -# let n = x.len -# result = bit_reverse(x) +proc iterative_cooley_tukey[T: SomeFloat](x: openArray[Complex[T]]): seq[Complex[T]] = + let n = x.len + result = bit_reverse(x) + for i in 1..fastLog2(n): + let + stride = 2 ^ i + w = exp(-complex(0.0, 2.0) * PI / float(stride)) + for j in skip(toSeq(0..n - 1), stride): + var v = complex(1.0, 0.0) + let s = stride div 2 + for k in 0..s - 1: + result[k + j + s] = result[k + j] - v * result[k + j + s] + result[k + j] -= result[k + j + s] - result[k + j] + v *= w when isMainModule: randomize() let x = toSeq(1..64).map(i => complex64(rand(1.0), 0.0)) y = cooley_tukey(x) - # z = iterative_cooley_tukey(x) + z = iterative_cooley_tukey(x) t = dft(x) reference = dft_fftw3(x) + assert(bit_reverse(@[1, 2, 3, 4, 7, 11]) == @[7, 3, 11, 4, 1, 3]) assert(allClose(reference, y)) - # assert(allClose(reference, z)) + assert(allClose(reference, z)) assert(allClose(reference, t)) From f6e2d72b9950cecb080e26f6be53ee9779cf1313 Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Sun, 13 Oct 2019 19:28:25 -0400 Subject: [PATCH 4/5] Nim FFT documentation --- contents/cooley_tukey/code/nim/fft.nim | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/contents/cooley_tukey/code/nim/fft.nim b/contents/cooley_tukey/code/nim/fft.nim index bd7ad4266..3f93cb824 100644 --- a/contents/cooley_tukey/code/nim/fft.nim +++ b/contents/cooley_tukey/code/nim/fft.nim @@ -7,12 +7,15 @@ import complex # $ nimble install fftw3 import fftw3 -# For some reason this isn't in the Nim fftw3 bindings. +# For some reason this isn't in the Nim fftw3 bindings, but this is the +# current definition from the C header. const FFTW_FORWARD = -1 type InpNumber = SomeNumber | Complex[SomeFloat] -proc allClose[T: InpNumber](reference: openArray[T], calculated: openArray[T], thresh = 1.0e-11): bool = +proc allClose[T: InpNumber](reference: openArray[T], calculated: openArray[T], + thresh = 1.0e-11): bool = + ## Are all the numbers in the reference equal to those that were calculated within some threshold? if reference.len != calculated.len: return false result = true @@ -43,7 +46,9 @@ proc nim_to_fftw(inp: openArray[Complex64]): seq[fftw_complex] = proc fftw_to_nim(inp: openArray[fftw_complex]): seq[Complex64] = inp.map(v => fftw_to_nim(v)) +# Anything other than 64-bit floats isn't currently supported by the binding. proc dft_fftw3(inp: openArray[Complex64]): seq[Complex64] = + ## Perform a forward Fourier transform using a binding to FFTW3. var fftw_inp = nim_to_fftw(inp) fftw_out = newSeq[fftw_complex](inp.len) @@ -57,6 +62,7 @@ proc dft_fftw3(inp: openArray[Complex64]): seq[Complex64] = fftw_to_nim(fftw_out) proc dft[T: SomeFloat](x: openArray[Complex[T]]): seq[Complex[T]] = + ## Perform a forward Fourier transform using the naive quadratic algorithm. let n = x.len result = newSeq[Complex[T]](n) for i in 0..n - 1: @@ -64,6 +70,7 @@ proc dft[T: SomeFloat](x: openArray[Complex[T]]): seq[Complex[T]] = result[i] += x[k] * exp(-complex(0.0, 2.0) * PI * float(i * k / n)) proc cooley_tukey[T: SomeFloat](x: openArray[Complex[T]]): seq[Complex[T]] = + ## Perform a forward Fourier transform using the recursive Cooley-Tukey algorithm. let n = x.len if n <= 1: return toSeq(x) @@ -92,6 +99,7 @@ proc bit_reverse[T: InpNumber](x: openArray[T]): seq[T] = result[b] = x[k] proc iterative_cooley_tukey[T: SomeFloat](x: openArray[Complex[T]]): seq[Complex[T]] = + ## Perform a forward Fourier transform using the iterative Cooley-Tukey algorithm. let n = x.len result = bit_reverse(x) for i in 1..fastLog2(n): From e45518831b97b207afc84db05a275bd0cec08b69 Mon Sep 17 00:00:00 2001 From: Eric Berquist Date: Sun, 13 Oct 2019 19:34:29 -0400 Subject: [PATCH 5/5] Add Nim FFT to book --- contents/cooley_tukey/cooley_tukey.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/contents/cooley_tukey/cooley_tukey.md b/contents/cooley_tukey/cooley_tukey.md index 1b7066063..4c5bc2af0 100644 --- a/contents/cooley_tukey/cooley_tukey.md +++ b/contents/cooley_tukey/cooley_tukey.md @@ -87,6 +87,8 @@ For some reason, though, putting code to this transformation really helped me fi [import:15-74, lang:"asm-x64"](code/asm-x64/fft.s) {% sample lang="js" %} [import:3-15, lang:"javascript"](code/javascript/fft.js) +{% sample lang="nim" %} +[import:64-70, lang:"nim"](code/nim/fft.nim) {% endmethod %} In this function, we define `n` to be a set of integers from $$0 \rightarrow N-1$$ and arrange them to be a column. @@ -138,6 +140,8 @@ In the end, the code looks like: [import:76-165, lang:"asm-x64"](code/asm-x64/fft.s) {% sample lang="js" %} [import:17-39, lang="javascript"](code/javascript/fft.js) +{% sample lang="nim" %} +[import:72-85, lang:"nim"](code/nim/fft.nim) {% endmethod %} As a side note, we are enforcing that the array must be a power of 2 for the operation to work. @@ -249,6 +253,8 @@ Some rather impressive scratch code was submitted by Jie and can be found here: [import, lang:"asm-x64"](code/asm-x64/fft.s) {% sample lang="js" %} [import, lang:"javascript"](code/javascript/fft.js) +{% sample lang="nim" %} +[import, lang:"nim"](code/nim/fft.nim) {% endmethod %}