From 944d2c5b83f4214d386f68ae3b693191c5630c9f Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Wed, 5 Feb 2014 11:40:23 +0100 Subject: [PATCH 01/20] Added 2D Fourier transformation as num.fft2 and num.fft2_inv. --- doc/user-manual/fft.rst | 15 ++++ fft-init.lua | 161 +++++++++++++++++++++++++++++++++++++++- tests/fft-tests.lua | 24 +++++- 3 files changed, 198 insertions(+), 2 deletions(-) diff --git a/doc/user-manual/fft.rst b/doc/user-manual/fft.rst index 5bdce72ff..0b4703667 100644 --- a/doc/user-manual/fft.rst +++ b/doc/user-manual/fft.rst @@ -130,6 +130,21 @@ As shown in the example above, you can use the Lua operator '#' to obtain the si vt = num.fftinv(ft) -- we perform the inverse Fourier transform -- now vt is a vector of the same size of v +.. function:: fft2(mat) + + Calculates the 2D Fourier transformation of the real or complex matrix ``mat``. The result is returned as a complex matrix. + Typical usage is:: + + --creating a real or complex matrix + mat = matrix.new(5,5, |i,j| i+j) + + --calling the fft2 function + num.fft2(mat) + +.. function:: fft2_inv(mat) + + Calculates the inverse 2D Fourier transformation of the real or complex matrix ``mat``. The result is returned as a complex matrix. See :func:`fft2` for usage. + FFT example ----------- diff --git a/fft-init.lua b/fft-init.lua index d7bad58b3..bf3df372b 100644 --- a/fft-init.lua +++ b/fft-init.lua @@ -82,6 +82,72 @@ ffi.cdef [[ const gsl_fft_halfcomplex_wavetable * wavetable, gsl_fft_real_workspace * work); +typedef double * gsl_complex_packed_array ; +typedef float * gsl_complex_packed_array_float ; +typedef long double * gsl_complex_packed_array_long_double ; + + int gsl_fft_complex_radix2_forward (gsl_complex_packed_array data, + const size_t stride, + const size_t n); + + int gsl_fft_complex_radix2_backward (gsl_complex_packed_array data, + const size_t stride, + const size_t n); + + int gsl_fft_complex_radix2_inverse (gsl_complex_packed_array data, + const size_t stride, + const size_t n); + +typedef struct + { + size_t n; + size_t nf; + size_t factor[64]; + gsl_complex *twiddle[64]; + gsl_complex *trig; + } +gsl_fft_complex_wavetable; + +typedef struct +{ + size_t n; + double *scratch; +} +gsl_fft_complex_workspace; + + +gsl_fft_complex_wavetable *gsl_fft_complex_wavetable_alloc (size_t n); + +void gsl_fft_complex_wavetable_free (gsl_fft_complex_wavetable * wavetable); + +gsl_fft_complex_workspace *gsl_fft_complex_workspace_alloc (size_t n); + +void gsl_fft_complex_workspace_free (gsl_fft_complex_workspace * workspace); + +int gsl_fft_complex_memcpy (gsl_fft_complex_wavetable * dest, + gsl_fft_complex_wavetable * src); + + +int gsl_fft_complex_forward (gsl_complex_packed_array data, + const size_t stride, + const size_t n, + const gsl_fft_complex_wavetable * wavetable, + gsl_fft_complex_workspace * work); + +int gsl_fft_complex_backward (gsl_complex_packed_array data, + const size_t stride, + const size_t n, + const gsl_fft_complex_wavetable * wavetable, + gsl_fft_complex_workspace * work); + +int gsl_fft_complex_inverse (gsl_complex_packed_array data, + const size_t stride, + const size_t n, + const gsl_fft_complex_wavetable * wavetable, + gsl_fft_complex_workspace * work); + + + ]] local fft_hc = ffi.typeof('fft_hc') @@ -110,7 +176,9 @@ end local cache_allocator = { real_wavetable = res_allocator('real_wavetable'), halfcomplex_wavetable = res_allocator('halfcomplex_wavetable'), - real_workspace = res_allocator('real_workspace') + real_workspace = res_allocator('real_workspace'), + complex_wavetable = res_allocator('complex_wavetable'), + complex_workspace = res_allocator('complex_workspace') } local function get_resource(name, n) @@ -180,6 +248,97 @@ function num.fftinv(ft, ip) return gsl_matrix(n, 1, stride, data, b, 1) end +function num.fft2(mat) + local n1 = tonumber(mat.size1) + local n2 = tonumber(mat.size2) + + if ffi.istype(gsl_matrix, mat) then + return num.fft2(matrix.cnew(n1, n2, |i,j| mat:get(i,j))) + else + mat = mat:copy() + end + + local newm = {} + local finm = {} + local wt,ws + + for i = 1,n2 do + local vec = mat:col(i) + local b,data,stride = vec.block, vec.data, vec.tda + if is_two_power(n1) then + gsl_check(gsl.gsl_fft_complex_radix2_forward(data, stride, n1 )) + else + wt = get_resource('complex_wavetable', n1) + ws = get_resource('complex_workspace', n1) + gsl_check(gsl.gsl_fft_complex_forward(data, stride, n1, wt, ws)) + end + table.insert(newm, vec) + end + newm = matrix.cdef(newm) + + for i = 1,n1 do + local vec = newm:col(i) + local b, data, stride = vec.block, vec.data, vec.tda + if is_two_power(n2) then + gsl_check(gsl.gsl_fft_complex_radix2_forward(data, stride, n2 )) + else + if n1 ~= n2 then + wt = get_resource('complex_wavetable', n2) + ws = get_resource('complex_workspace', n2) + end + gsl_check(gsl.gsl_fft_complex_forward(data, stride, n2, wt, ws)) + end + table.insert(finm, vec) + end + return matrix.cdef(finm) +end + +function num.fft2_inv(mat) + local n1 = tonumber(mat.size1) + local n2 = tonumber(mat.size2) + + if ffi.istype(gsl_matrix, mat) then + return num.fft2_inv(matrix.cnew(n1, n2, |i,j| mat:get(i,j))) + else + mat = mat:copy() + end + + local newm = {} + local finm = {} + local wt,ws + + for i = 1,n2 do + local vec = mat:col(i) + local b,data,stride = vec.block, vec.data, vec.tda + if is_two_power(n1) then + gsl_check(gsl.gsl_fft_complex_radix2_inverse(data, stride, n1 )) + else + wt = get_resource('complex_wavetable', n1) + ws = get_resource('complex_workspace', n1) + gsl_check(gsl.gsl_fft_complex_inverse(data, stride, n1, wt, ws)) + end + table.insert(newm, vec) + end + newm = matrix.cdef(newm) + + for i = 1,n1 do + local vec = newm:col(i) + local b, data, stride = vec.block, vec.data, vec.tda + if is_two_power(n2) then + gsl_check(gsl.gsl_fft_complex_radix2_inverse(data, stride, n2 )) + else + if n1 ~= n2 then + wt = get_resource('complex_wavetable', n2) + ws = get_resource('complex_workspace', n2) + end + gsl_check(gsl.gsl_fft_complex_inverse(data, stride, n2, wt, ws)) + end + table.insert(finm, vec) + end + return matrix.cdef(finm) +end + + local function halfcomplex_radix2_index(n, stride, k) if k < 0 or k >= n then error('invalid halfcomplex index', 2) end local half_n = n/2 diff --git a/tests/fft-tests.lua b/tests/fft-tests.lua index b777797e3..ff5fea709 100644 --- a/tests/fft-tests.lua +++ b/tests/fft-tests.lua @@ -164,9 +164,31 @@ local function test1_stride() end end +local function test_fft2() + local realmat_radix = matrix.new(4,4, |i,j|i) + local realmat_notradix = matrix.new(5,5, |i,j|i) + + local compmat_radix = matrix.cnew(4,4, |i,j|i+j*1i) + local compmat_notradix = matrix.cnew(5,5, |i,j|i+j*1i) + + print(num.fft2(realmat_radix)) + print(num.fft2(realmat_notradix)) + + print(num.fft2(compmat_radix)) + print(num.fft2(compmat_notradix)) + + print(num.fft2_inv(realmat_radix)) + print(num.fft2_inv(realmat_notradix)) + + print(num.fft2_inv(compmat_radix)) + print(num.fft2_inv(compmat_notradix)) + +end + return {test1= test1, test2= test1_radix2, test3= test1_ip_radix2, test4= test1_ip, test5= test2, - test6= test1_stride} + test6= test1_stride, + test7 = test_fft2} From 84d781fb0f3cf77bb26823f2e24e0f90b97ec05c Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Thu, 6 Feb 2014 12:55:26 +0100 Subject: [PATCH 02/20] Revert "Added 2D Fourier transformation as num.fft2 and num.fft2_inv." This reverts commit 944d2c5b83f4214d386f68ae3b693191c5630c9f. --- doc/user-manual/fft.rst | 15 ---- fft-init.lua | 161 +--------------------------------------- tests/fft-tests.lua | 24 +----- 3 files changed, 2 insertions(+), 198 deletions(-) diff --git a/doc/user-manual/fft.rst b/doc/user-manual/fft.rst index 0b4703667..5bdce72ff 100644 --- a/doc/user-manual/fft.rst +++ b/doc/user-manual/fft.rst @@ -130,21 +130,6 @@ As shown in the example above, you can use the Lua operator '#' to obtain the si vt = num.fftinv(ft) -- we perform the inverse Fourier transform -- now vt is a vector of the same size of v -.. function:: fft2(mat) - - Calculates the 2D Fourier transformation of the real or complex matrix ``mat``. The result is returned as a complex matrix. - Typical usage is:: - - --creating a real or complex matrix - mat = matrix.new(5,5, |i,j| i+j) - - --calling the fft2 function - num.fft2(mat) - -.. function:: fft2_inv(mat) - - Calculates the inverse 2D Fourier transformation of the real or complex matrix ``mat``. The result is returned as a complex matrix. See :func:`fft2` for usage. - FFT example ----------- diff --git a/fft-init.lua b/fft-init.lua index bf3df372b..d7bad58b3 100644 --- a/fft-init.lua +++ b/fft-init.lua @@ -82,72 +82,6 @@ ffi.cdef [[ const gsl_fft_halfcomplex_wavetable * wavetable, gsl_fft_real_workspace * work); -typedef double * gsl_complex_packed_array ; -typedef float * gsl_complex_packed_array_float ; -typedef long double * gsl_complex_packed_array_long_double ; - - int gsl_fft_complex_radix2_forward (gsl_complex_packed_array data, - const size_t stride, - const size_t n); - - int gsl_fft_complex_radix2_backward (gsl_complex_packed_array data, - const size_t stride, - const size_t n); - - int gsl_fft_complex_radix2_inverse (gsl_complex_packed_array data, - const size_t stride, - const size_t n); - -typedef struct - { - size_t n; - size_t nf; - size_t factor[64]; - gsl_complex *twiddle[64]; - gsl_complex *trig; - } -gsl_fft_complex_wavetable; - -typedef struct -{ - size_t n; - double *scratch; -} -gsl_fft_complex_workspace; - - -gsl_fft_complex_wavetable *gsl_fft_complex_wavetable_alloc (size_t n); - -void gsl_fft_complex_wavetable_free (gsl_fft_complex_wavetable * wavetable); - -gsl_fft_complex_workspace *gsl_fft_complex_workspace_alloc (size_t n); - -void gsl_fft_complex_workspace_free (gsl_fft_complex_workspace * workspace); - -int gsl_fft_complex_memcpy (gsl_fft_complex_wavetable * dest, - gsl_fft_complex_wavetable * src); - - -int gsl_fft_complex_forward (gsl_complex_packed_array data, - const size_t stride, - const size_t n, - const gsl_fft_complex_wavetable * wavetable, - gsl_fft_complex_workspace * work); - -int gsl_fft_complex_backward (gsl_complex_packed_array data, - const size_t stride, - const size_t n, - const gsl_fft_complex_wavetable * wavetable, - gsl_fft_complex_workspace * work); - -int gsl_fft_complex_inverse (gsl_complex_packed_array data, - const size_t stride, - const size_t n, - const gsl_fft_complex_wavetable * wavetable, - gsl_fft_complex_workspace * work); - - - ]] local fft_hc = ffi.typeof('fft_hc') @@ -176,9 +110,7 @@ end local cache_allocator = { real_wavetable = res_allocator('real_wavetable'), halfcomplex_wavetable = res_allocator('halfcomplex_wavetable'), - real_workspace = res_allocator('real_workspace'), - complex_wavetable = res_allocator('complex_wavetable'), - complex_workspace = res_allocator('complex_workspace') + real_workspace = res_allocator('real_workspace') } local function get_resource(name, n) @@ -248,97 +180,6 @@ function num.fftinv(ft, ip) return gsl_matrix(n, 1, stride, data, b, 1) end -function num.fft2(mat) - local n1 = tonumber(mat.size1) - local n2 = tonumber(mat.size2) - - if ffi.istype(gsl_matrix, mat) then - return num.fft2(matrix.cnew(n1, n2, |i,j| mat:get(i,j))) - else - mat = mat:copy() - end - - local newm = {} - local finm = {} - local wt,ws - - for i = 1,n2 do - local vec = mat:col(i) - local b,data,stride = vec.block, vec.data, vec.tda - if is_two_power(n1) then - gsl_check(gsl.gsl_fft_complex_radix2_forward(data, stride, n1 )) - else - wt = get_resource('complex_wavetable', n1) - ws = get_resource('complex_workspace', n1) - gsl_check(gsl.gsl_fft_complex_forward(data, stride, n1, wt, ws)) - end - table.insert(newm, vec) - end - newm = matrix.cdef(newm) - - for i = 1,n1 do - local vec = newm:col(i) - local b, data, stride = vec.block, vec.data, vec.tda - if is_two_power(n2) then - gsl_check(gsl.gsl_fft_complex_radix2_forward(data, stride, n2 )) - else - if n1 ~= n2 then - wt = get_resource('complex_wavetable', n2) - ws = get_resource('complex_workspace', n2) - end - gsl_check(gsl.gsl_fft_complex_forward(data, stride, n2, wt, ws)) - end - table.insert(finm, vec) - end - return matrix.cdef(finm) -end - -function num.fft2_inv(mat) - local n1 = tonumber(mat.size1) - local n2 = tonumber(mat.size2) - - if ffi.istype(gsl_matrix, mat) then - return num.fft2_inv(matrix.cnew(n1, n2, |i,j| mat:get(i,j))) - else - mat = mat:copy() - end - - local newm = {} - local finm = {} - local wt,ws - - for i = 1,n2 do - local vec = mat:col(i) - local b,data,stride = vec.block, vec.data, vec.tda - if is_two_power(n1) then - gsl_check(gsl.gsl_fft_complex_radix2_inverse(data, stride, n1 )) - else - wt = get_resource('complex_wavetable', n1) - ws = get_resource('complex_workspace', n1) - gsl_check(gsl.gsl_fft_complex_inverse(data, stride, n1, wt, ws)) - end - table.insert(newm, vec) - end - newm = matrix.cdef(newm) - - for i = 1,n1 do - local vec = newm:col(i) - local b, data, stride = vec.block, vec.data, vec.tda - if is_two_power(n2) then - gsl_check(gsl.gsl_fft_complex_radix2_inverse(data, stride, n2 )) - else - if n1 ~= n2 then - wt = get_resource('complex_wavetable', n2) - ws = get_resource('complex_workspace', n2) - end - gsl_check(gsl.gsl_fft_complex_inverse(data, stride, n2, wt, ws)) - end - table.insert(finm, vec) - end - return matrix.cdef(finm) -end - - local function halfcomplex_radix2_index(n, stride, k) if k < 0 or k >= n then error('invalid halfcomplex index', 2) end local half_n = n/2 diff --git a/tests/fft-tests.lua b/tests/fft-tests.lua index ff5fea709..b777797e3 100644 --- a/tests/fft-tests.lua +++ b/tests/fft-tests.lua @@ -164,31 +164,9 @@ local function test1_stride() end end -local function test_fft2() - local realmat_radix = matrix.new(4,4, |i,j|i) - local realmat_notradix = matrix.new(5,5, |i,j|i) - - local compmat_radix = matrix.cnew(4,4, |i,j|i+j*1i) - local compmat_notradix = matrix.cnew(5,5, |i,j|i+j*1i) - - print(num.fft2(realmat_radix)) - print(num.fft2(realmat_notradix)) - - print(num.fft2(compmat_radix)) - print(num.fft2(compmat_notradix)) - - print(num.fft2_inv(realmat_radix)) - print(num.fft2_inv(realmat_notradix)) - - print(num.fft2_inv(compmat_radix)) - print(num.fft2_inv(compmat_notradix)) - -end - return {test1= test1, test2= test1_radix2, test3= test1_ip_radix2, test4= test1_ip, test5= test2, - test6= test1_stride, - test7 = test_fft2} + test6= test1_stride} From 6684c13189330931a7d1dcc1b77a47f70c3f2969 Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Thu, 6 Feb 2014 18:12:40 +0100 Subject: [PATCH 03/20] Implemented fft, fft2, fft3 and fftn as well as their inverse and real counterparts based on the FFTW3 library via FFI. --- fft.lua | 252 +++++++++++++++++++++++++++++++++ fftw3/cdefs.lua | 346 ++++++++++++++++++++++++++++++++++++++++++++++ fftw3/defines.lua | 36 +++++ fftw3/init.lua | 249 +++++++++++++++++++++++++++++++++ 4 files changed, 883 insertions(+) create mode 100644 fft.lua create mode 100644 fftw3/cdefs.lua create mode 100644 fftw3/defines.lua create mode 100644 fftw3/init.lua diff --git a/fft.lua b/fft.lua new file mode 100644 index 000000000..2b9ecc6fa --- /dev/null +++ b/fft.lua @@ -0,0 +1,252 @@ +fftw = require'fftw3.init' +local gsl_matrix = ffi.typeof('gsl_matrix') +local gsl_matrix_complex = ffi.typeof('gsl_matrix_complex') + +function num.fft(vec) + if ffi.istype(gsl_matrix, vec) then + return num.rfft(vec) + else + local n = tonumber(vec.size1) + local outvec = matrix.cnew(n, 1) + local output = ffi.cast("fftw_complex*",outvec.data) + local input = ffi.cast("fftw_complex*",vec.data) + local plan = fftw.plan_dft_1d(n, input, output, fftw.FORWARD, fftw.MEASURE) + fftw.execute(plan) + fftw.destroy_plan(plan) + return outvec + end +end + +function num.fft_inv(vec) + if ffi.istype(gsl_matrix, vec) then + return num.rfft_inv(vec) + else + local n = tonumber(vec.size1) + local outvec = matrix.cnew(n, 1) + local output = ffi.cast("fftw_complex*",outvec.data) + local input = ffi.cast("fftw_complex*",vec.data) + local plan = fftw.plan_dft_1d(n, input, output, fftw.BACKWARD, fftw.MEASURE) + fftw.execute(plan) + fftw.destroy_plan(plan) + return outvec/n + end +end + +function num.rfft(vec) + if ffi.istype(gsl_matrix, vec) then + local n = tonumber(vec.size1) + local outvec = matrix.cnew(math.floor(n/2)+1, 1) + local output = ffi.cast("fftw_complex*",outvec.data) + local plan = fftw.plan_dft_r2c_1d(n, vec.data, output, fftw.MEASURE ) + fftw.execute(plan) + fftw.destroy_plan(plan) + return outvec + else + error("Cannot handle nonreal input data.") + end +end + +function num.rfft_inv(vec) + if ffi.istype(gsl_matrix_complex, vec) then + local n = tonumber(vec.size1) + local nnew = (n-1)*2 + local outvec = matrix.new(nnew, 1) + local input = ffi.cast("fftw_complex*",vec.data) + local plan = fftw.plan_dft_c2r_1d(nnew, input, outvec.data, fftw.MEASURE ) + fftw.execute(plan) + fftw.destroy_plan(plan) + return outvec/nnew + else + error("Cannot handle noncomplex input data.") + end +end + +---------------------------------- + +function num.fft2(mat) + if ffi.istype(gsl_matrix, mat) then + return num.rfft2(mat) + else + local n1 = tonumber(mat.size1) + local n2 = tonumber(mat.size2) + local outmat = matrix.cnew(n1, n2) + local output = ffi.cast("fftw_complex*",outmat.data) + local input = ffi.cast("fftw_complex*", mat.data) + local plan = fftw.plan_dft_2d(n1,n2, input, output, fftw.FORWARD, fftw.MEASURE) + fftw.execute(plan) + fftw.destroy_plan(plan) + return outmat + end +end + +function num.fft2_inv(mat) + if ffi.istype(gsl_matrix, mat) then + return num.rfft2_inv(mat) + else + local n1 = tonumber(mat.size1) + local n2 = tonumber(mat.size2) + local outmat = matrix.cnew(n1, n2) + local output = ffi.cast("fftw_complex*",outmat.data) + local input = ffi.cast("fftw_complex*", mat.data) + local plan = fftw.plan_dft_2d(n1,n2, input, output, fftw.BACKWARD, fftw.MEASURE) + fftw.execute(plan) + fftw.destroy_plan(plan) + return outmat + end +end + +function num.rfft2(mat) + if ffi.istype(gsl_matrix, mat) then + local n1 = tonumber(mat.size1) + local n2 = tonumber(mat.size2) + local outmat = matrix.cnew(n1, math.floor(n2/2)+1) + local output = ffi.cast("fftw_complex*",outmat.data) + local plan = fftw.plan_dft_r2c_2d(n1,n2, mat.data, output,fftw.MEASURE) + fftw.execute(plan) + fftw.destroy_plan(plan) + return outmat + else + error("Cannot handle nonreal input matrix.") + end +end + +function num.rfft2_inv(mat) + if ffi.istype(gsl_matrix_complex, mat) then + local n1 = tonumber(mat.size1) + local n2 = tonumber(mat.size2) + local n2new = (n2-1)*2 + local outmat = matrix.new(n1, n2new) + local input = ffi.cast("fftw_complex*",mat.data) + local plan = fftw.plan_dft_c2r_2d(n1,n2new, input, outmat.data,fftw.MEASURE) + fftw.execute(plan) + fftw.destroy_plan(plan) + return outmat/(n1*n2new) + else + error("Cannot handle noncomplex input matrix.") + end +end + +-------------------------------------- + +function num.fft3(datavec, n1, n2, n3) + if ffi.istype(gsl_matrix, datavec) then + return num.rfft3(datavec, n1, n2, n3) + else + local outvec = matrix.cnew(n1*n2*n3, 1) + local output = ffi.cast("fftw_complex*",outvec.data) + local input = ffi.cast("fftw_complex*", datavec.data) + local plan = fftw.plan_dft_3d(n1,n2,n3, input, output, fftw.FORWARD, fftw.MEASURE) + fftw.execute(plan) + fftw.destroy_plan(plan) + return outvec + end +end + +function num.fft3_inv(datavec, n1, n2, n3) + if ffi.istype(gsl_matrix, datavec) then + return num.rfft3_inv(datavec, n1, n2, n3) + else + local outvec = matrix.cnew(n1*n2*n3, 1) + local output = ffi.cast("fftw_complex*",outvec.data) + local input = ffi.cast("fftw_complex*", datavec.data) + local plan = fftw.plan_dft_3d(n1,n2,n3, input, output, fftw.BACKWARD, fftw.MEASURE) + fftw.execute(plan) + fftw.destroy_plan(plan) + return outvec/(n1*n2*n3) + end +end + +function num.rfft3(datavec, n1, n2, n3) + if ffi.istype(gsl_matrix, datavec) then + local outvec = matrix.cnew(n1*n2*(math.floor(n3/2)+1), 1) + local output = ffi.cast("fftw_complex*",outvec.data) + local plan = fftw.plan_dft_r2c_3d(n1,n2,n3, datavec.data, output, fftw.MEASURE) + fftw.execute(plan) + fftw.destroy_plan(plan) + return outvec + else + error("Cannot handle nonreal input data.") + end +end + +function num.rfft3_inv(datavec, n1, n2, n3) + if ffi.istype(gsl_matrix, datavec) then + local outvec = matrix.new(n1*n2*n3, 1) + local input = ffi.cast("fftw_complex*",datavec.data) + local plan = fftw.plan_dft_r2c_3d(n1,n2,n3, input, outvec.data, fftw.MEASURE) + fftw.execute(plan) + fftw.destroy_plan(plan) + return outvec/(n1*n2*n3) + else + error("Cannot handle nonreal input data.") + end +end + +-------------------------------------- + +function num.fftn(datavec, dimensionvec) + if ffi.istype(gsl_matrix, datavec) then + return num.rfftn(datavec, dimensionvec) + else + local outvec = matrix.cnew(tonumber(datavec.size1), 1) + local output = ffi.cast("fftw_complex*",outvec.data) + local input = ffi.cast("fftw_complex*", datavec.data) + local plan = fftw.plan_dft(#dimensionvec, dimensionvec.data, input, output, fftw.FORWARD, fftw.MEASURE) + fftw.execute(plan) + fftw.destroy_plan(plan) + return outvec + end +end + +function num.fftn_inv(datavec, dimensionvec) + if ffi.istype(gsl_matrix, datavec) then + return num.rfftn_inv(datavec, dimensionvec) + else + local outvec = matrix.cnew(tonumber(datavec.size1), 1) + local output = ffi.cast("fftw_complex*",outvec.data) + local input = ffi.cast("fftw_complex*", datavec.data) + local plan = fftw.plan_dft(#dimensionvec, dimensionvec.data, input, output, fftw.BACKWARD, fftw.MEASURE) + fftw.execute(plan) + fftw.destroy_plan(plan) + return outvec + end +end + +function num.rfftn(datavec, dimensionvec) + if ffi.istype(gsl_matrix, datavec) then + local d = #dimensionvec + local newsize = 1 + for i=1,d-1 do + newsize = newsize * dimensionvec[i] + end + newsize = newsize * math.floor(dimensionvec[d]/2)+1 + + local outvec = matrix.cnew(newsize, 1) + local output = ffi.cast("fftw_complex*",outvec.data) + local plan = fftw.plan_dft_r2c(d, dimensionvec.data, datavec.data, output, fftw.MEASURE) + fftw.execute(plan) + fftw.destroy_plan(plan) + return outvec + else + error("Cannot handle nonreal input data.") + end +end + +function num.rfftn_inv(datavec, dimensionvec) + if ffi.istype(gsl_matrix, datavec) then + local d = #dimensionvec + local newsize = 1 + for i=1,d do + newsize = newsize * dimensionvec[i] + end + + local outvec = matrix.new(newsize, 1) + local output = ffi.cast("fftw_complex*",outvec.data) + local plan = fftw.plan_dft_r2c(d, dimensionvec.data, datavec.data, output, fftw.MEASURE) + fftw.execute(plan) + fftw.destroy_plan(plan) + return outvec + else + error("Cannot handle nonreal input data.") + end +end diff --git a/fftw3/cdefs.lua b/fftw3/cdefs.lua new file mode 100644 index 000000000..e1512f151 --- /dev/null +++ b/fftw3/cdefs.lua @@ -0,0 +1,346 @@ +-- Cut and paste from the C preprocessor output +-- Removed inline/defined functions which are not supported by luajit +-- Instead, those are defined into defines.lua +-- Note there are some tests here and there to stay cross-platform + +local ffi = require 'ffi' + +ffi.cdef[[ +typedef struct _FILE FILE; +typedef long double __float128; +]] + + +ffi.cdef[[ + +enum fftw_r2r_kind_do_not_use_me { + FFTW_R2HC=0, FFTW_HC2R=1, FFTW_DHT=2, + FFTW_REDFT00=3, FFTW_REDFT01=4, FFTW_REDFT10=5, FFTW_REDFT11=6, + FFTW_RODFT00=7, FFTW_RODFT01=8, FFTW_RODFT10=9, FFTW_RODFT11=10 +}; + +struct fftw_iodim_do_not_use_me { + int n; + int is; + int os; +}; + +struct fftw_iodim64_do_not_use_me { + ptrdiff_t n; + ptrdiff_t is; + ptrdiff_t os; +}; + +typedef void (*fftw_write_char_func_do_not_use_me)(char c, void *); +typedef int (*fftw_read_char_func_do_not_use_me)(void *); + +typedef double fftw_complex[2]; +typedef struct fftw_plan_s *fftw_plan; +typedef struct fftw_iodim_do_not_use_me fftw_iodim; +typedef struct fftw_iodim64_do_not_use_me fftw_iodim64; +typedef enum fftw_r2r_kind_do_not_use_me fftw_r2r_kind; +typedef fftw_write_char_func_do_not_use_me fftw_write_char_func; +typedef fftw_read_char_func_do_not_use_me fftw_read_char_func; +extern void fftw_execute(const fftw_plan p); +extern fftw_plan fftw_plan_dft(int rank, const int *n, fftw_complex *in, fftw_complex *out, int sign, unsigned flags); +extern fftw_plan fftw_plan_dft_1d(int n, fftw_complex *in, fftw_complex *out, int sign, unsigned flags); +extern fftw_plan fftw_plan_dft_2d(int n0, int n1, fftw_complex *in, fftw_complex *out, int sign, unsigned flags); +extern fftw_plan fftw_plan_dft_3d(int n0, int n1, int n2, fftw_complex *in, fftw_complex *out, int sign, unsigned flags); +extern fftw_plan fftw_plan_many_dft(int rank, const int *n, int howmany, fftw_complex *in, const int *inembed, int istride, int idist, fftw_complex *out, const int *onembed, int ostride, int odist, int sign, unsigned flags); +extern fftw_plan fftw_plan_guru_dft(int rank, const fftw_iodim *dims, int howmany_rank, const fftw_iodim *howmany_dims, fftw_complex *in, fftw_complex *out, int sign, unsigned flags); +extern fftw_plan fftw_plan_guru_split_dft(int rank, const fftw_iodim *dims, int howmany_rank, const fftw_iodim *howmany_dims, double *ri, double *ii, double *ro, double *io, unsigned flags); +extern fftw_plan fftw_plan_guru64_dft(int rank, const fftw_iodim64 *dims, int howmany_rank, const fftw_iodim64 *howmany_dims, fftw_complex *in, fftw_complex *out, int sign, unsigned flags); +extern fftw_plan fftw_plan_guru64_split_dft(int rank, const fftw_iodim64 *dims, int howmany_rank, const fftw_iodim64 *howmany_dims, double *ri, double *ii, double *ro, double *io, unsigned flags); +extern void fftw_execute_dft(const fftw_plan p, fftw_complex *in, fftw_complex *out); +extern void fftw_execute_split_dft(const fftw_plan p, double *ri, double *ii, double *ro, double *io); +extern fftw_plan fftw_plan_many_dft_r2c(int rank, const int *n, int howmany, double *in, const int *inembed, int istride, int idist, fftw_complex *out, const int *onembed, int ostride, int odist, unsigned flags); +extern fftw_plan fftw_plan_dft_r2c(int rank, const int *n, double *in, fftw_complex *out, unsigned flags); +extern fftw_plan fftw_plan_dft_r2c_1d(int n,double *in,fftw_complex *out,unsigned flags); +extern fftw_plan fftw_plan_dft_r2c_2d(int n0, int n1, double *in, fftw_complex *out, unsigned flags); +extern fftw_plan fftw_plan_dft_r2c_3d(int n0, int n1, int n2, double *in, fftw_complex *out, unsigned flags); +extern fftw_plan fftw_plan_many_dft_c2r(int rank, const int *n, int howmany, fftw_complex *in, const int *inembed, int istride, int idist, double *out, const int *onembed, int ostride, int odist, unsigned flags); +extern fftw_plan fftw_plan_dft_c2r(int rank, const int *n, fftw_complex *in, double *out, unsigned flags); +extern fftw_plan fftw_plan_dft_c2r_1d(int n,fftw_complex *in,double *out,unsigned flags); +extern fftw_plan fftw_plan_dft_c2r_2d(int n0, int n1, fftw_complex *in, double *out, unsigned flags); +extern fftw_plan fftw_plan_dft_c2r_3d(int n0, int n1, int n2, fftw_complex *in, double *out, unsigned flags); +extern fftw_plan fftw_plan_guru_dft_r2c(int rank, const fftw_iodim *dims, int howmany_rank, const fftw_iodim *howmany_dims, double *in, fftw_complex *out, unsigned flags); +extern fftw_plan fftw_plan_guru_dft_c2r(int rank, const fftw_iodim *dims, int howmany_rank, const fftw_iodim *howmany_dims, fftw_complex *in, double *out, unsigned flags); +extern fftw_plan fftw_plan_guru_split_dft_r2c( int rank, const fftw_iodim *dims, int howmany_rank, const fftw_iodim *howmany_dims, double *in, double *ro, double *io, unsigned flags); +extern fftw_plan fftw_plan_guru_split_dft_c2r( int rank, const fftw_iodim *dims, int howmany_rank, const fftw_iodim *howmany_dims, double *ri, double *ii, double *out, unsigned flags); +extern fftw_plan fftw_plan_guru64_dft_r2c(int rank, const fftw_iodim64 *dims, int howmany_rank, const fftw_iodim64 *howmany_dims, double *in, fftw_complex *out, unsigned flags); +extern fftw_plan fftw_plan_guru64_dft_c2r(int rank, const fftw_iodim64 *dims, int howmany_rank, const fftw_iodim64 *howmany_dims, fftw_complex *in, double *out, unsigned flags); +extern fftw_plan fftw_plan_guru64_split_dft_r2c( int rank, const fftw_iodim64 *dims, int howmany_rank, const fftw_iodim64 *howmany_dims, double *in, double *ro, double *io, unsigned flags); +extern fftw_plan fftw_plan_guru64_split_dft_c2r( int rank, const fftw_iodim64 *dims, int howmany_rank, const fftw_iodim64 *howmany_dims, double *ri, double *ii, double *out, unsigned flags); +extern void fftw_execute_dft_r2c(const fftw_plan p, double *in, fftw_complex *out); +extern void fftw_execute_dft_c2r(const fftw_plan p, fftw_complex *in, double *out); +extern void fftw_execute_split_dft_r2c(const fftw_plan p, double *in, double *ro, double *io); +extern void fftw_execute_split_dft_c2r(const fftw_plan p, double *ri, double *ii, double *out); +extern fftw_plan fftw_plan_many_r2r(int rank, const int *n, int howmany, double *in, const int *inembed, int istride, int idist, double *out, const int *onembed, int ostride, int odist, const fftw_r2r_kind *kind, unsigned flags); +extern fftw_plan fftw_plan_r2r(int rank, const int *n, double *in, double *out, const fftw_r2r_kind *kind, unsigned flags); +extern fftw_plan fftw_plan_r2r_1d(int n, double *in, double *out, fftw_r2r_kind kind, unsigned flags); +extern fftw_plan fftw_plan_r2r_2d(int n0, int n1, double *in, double *out, fftw_r2r_kind kind0, fftw_r2r_kind kind1, unsigned flags); +extern fftw_plan fftw_plan_r2r_3d(int n0, int n1, int n2, double *in, double *out, fftw_r2r_kind kind0, fftw_r2r_kind kind1, fftw_r2r_kind kind2, unsigned flags); +extern fftw_plan fftw_plan_guru_r2r(int rank, const fftw_iodim *dims, int howmany_rank, const fftw_iodim *howmany_dims, double *in, double *out, const fftw_r2r_kind *kind, unsigned flags); +extern fftw_plan fftw_plan_guru64_r2r(int rank, const fftw_iodim64 *dims, int howmany_rank, const fftw_iodim64 *howmany_dims, double *in, double *out, const fftw_r2r_kind *kind, unsigned flags); +extern void fftw_execute_r2r(const fftw_plan p, double *in, double *out); +extern void fftw_destroy_plan(fftw_plan p); +extern void fftw_forget_wisdom(void); +extern void fftw_cleanup(void); +extern void fftw_set_timelimit(double t); +extern void fftw_plan_with_nthreads(int nthreads); +extern int fftw_init_threads(void); +extern void fftw_cleanup_threads(void); +extern int fftw_export_wisdom_to_filename(const char *filename); +extern void fftw_export_wisdom_to_file(FILE *output_file); +extern char *fftw_export_wisdom_to_string(void); +extern void fftw_export_wisdom(fftw_write_char_func write_char, void *data); +extern int fftw_import_system_wisdom(void); +extern int fftw_import_wisdom_from_filename(const char *filename); +extern int fftw_import_wisdom_from_file(FILE *input_file); +extern int fftw_import_wisdom_from_string(const char *input_string); +extern int fftw_import_wisdom(fftw_read_char_func read_char, void *data); +extern void fftw_fprint_plan(const fftw_plan p, FILE *output_file); +extern void fftw_print_plan(const fftw_plan p); +extern void *fftw_malloc(size_t n); +extern double *fftw_alloc_real(size_t n); +extern fftw_complex *fftw_alloc_complex(size_t n); +extern void fftw_free(void *p); +extern void fftw_flops(const fftw_plan p, double *add, double *mul, double *fmas); +extern double fftw_estimate_cost(const fftw_plan p); +extern double fftw_cost(const fftw_plan p); +extern const char fftw_version[]; +extern const char fftw_cc[]; +extern const char fftw_codelet_optim[]; +typedef float fftwf_complex[2]; +typedef struct fftwf_plan_s *fftwf_plan; +typedef struct fftw_iodim_do_not_use_me fftwf_iodim; +typedef struct fftw_iodim64_do_not_use_me fftwf_iodim64; +typedef enum fftw_r2r_kind_do_not_use_me fftwf_r2r_kind; +typedef fftw_write_char_func_do_not_use_me fftwf_write_char_func; +typedef fftw_read_char_func_do_not_use_me fftwf_read_char_func; +extern void fftwf_execute(const fftwf_plan p); +extern fftwf_plan fftwf_plan_dft(int rank, const int *n, fftwf_complex *in, fftwf_complex *out, int sign, unsigned flags); +extern fftwf_plan fftwf_plan_dft_1d(int n, fftwf_complex *in, fftwf_complex *out, int sign, unsigned flags); +extern fftwf_plan fftwf_plan_dft_2d(int n0, int n1, fftwf_complex *in, fftwf_complex *out, int sign, unsigned flags); +extern fftwf_plan fftwf_plan_dft_3d(int n0, int n1, int n2, fftwf_complex *in, fftwf_complex *out, int sign, unsigned flags); +extern fftwf_plan fftwf_plan_many_dft(int rank, const int *n, int howmany, fftwf_complex *in, const int *inembed, int istride, int idist, fftwf_complex *out, const int *onembed, int ostride, int odist, int sign, unsigned flags); +extern fftwf_plan fftwf_plan_guru_dft(int rank, const fftwf_iodim *dims, int howmany_rank, const fftwf_iodim *howmany_dims, fftwf_complex *in, fftwf_complex *out, int sign, unsigned flags); +extern fftwf_plan fftwf_plan_guru_split_dft(int rank, const fftwf_iodim *dims, int howmany_rank, const fftwf_iodim *howmany_dims, float *ri, float *ii, float *ro, float *io, unsigned flags); +extern fftwf_plan fftwf_plan_guru64_dft(int rank, const fftwf_iodim64 *dims, int howmany_rank, const fftwf_iodim64 *howmany_dims, fftwf_complex *in, fftwf_complex *out, int sign, unsigned flags); +extern fftwf_plan fftwf_plan_guru64_split_dft(int rank, const fftwf_iodim64 *dims, int howmany_rank, const fftwf_iodim64 *howmany_dims, float *ri, float *ii, float *ro, float *io, unsigned flags); +extern void fftwf_execute_dft(const fftwf_plan p, fftwf_complex *in, fftwf_complex *out); +extern void fftwf_execute_split_dft(const fftwf_plan p, float *ri, float *ii, float *ro, float *io); +extern fftwf_plan fftwf_plan_many_dft_r2c(int rank, const int *n, int howmany, float *in, const int *inembed, int istride, int idist, fftwf_complex *out, const int *onembed, int ostride, int odist, unsigned flags); +extern fftwf_plan fftwf_plan_dft_r2c(int rank, const int *n, float *in, fftwf_complex *out, unsigned flags); +extern fftwf_plan fftwf_plan_dft_r2c_1d(int n,float *in,fftwf_complex *out,unsigned flags); +extern fftwf_plan fftwf_plan_dft_r2c_2d(int n0, int n1, float *in, fftwf_complex *out, unsigned flags); +extern fftwf_plan fftwf_plan_dft_r2c_3d(int n0, int n1, int n2, float *in, fftwf_complex *out, unsigned flags); +extern fftwf_plan fftwf_plan_many_dft_c2r(int rank, const int *n, int howmany, fftwf_complex *in, const int *inembed, int istride, int idist, float *out, const int *onembed, int ostride, int odist, unsigned flags); +extern fftwf_plan fftwf_plan_dft_c2r(int rank, const int *n, fftwf_complex *in, float *out, unsigned flags); +extern fftwf_plan fftwf_plan_dft_c2r_1d(int n,fftwf_complex *in,float *out,unsigned flags); +extern fftwf_plan fftwf_plan_dft_c2r_2d(int n0, int n1, fftwf_complex *in, float *out, unsigned flags); +extern fftwf_plan fftwf_plan_dft_c2r_3d(int n0, int n1, int n2, fftwf_complex *in, float *out, unsigned flags); +extern fftwf_plan fftwf_plan_guru_dft_r2c(int rank, const fftwf_iodim *dims, int howmany_rank, const fftwf_iodim *howmany_dims, float *in, fftwf_complex *out, unsigned flags); +extern fftwf_plan fftwf_plan_guru_dft_c2r(int rank, const fftwf_iodim *dims, int howmany_rank, const fftwf_iodim *howmany_dims, fftwf_complex *in, float *out, unsigned flags); +extern fftwf_plan fftwf_plan_guru_split_dft_r2c( int rank, const fftwf_iodim *dims, int howmany_rank, const fftwf_iodim *howmany_dims, float *in, float *ro, float *io, unsigned flags); +extern fftwf_plan fftwf_plan_guru_split_dft_c2r( int rank, const fftwf_iodim *dims, int howmany_rank, const fftwf_iodim *howmany_dims, float *ri, float *ii, float *out, unsigned flags); +extern fftwf_plan fftwf_plan_guru64_dft_r2c(int rank, const fftwf_iodim64 *dims, int howmany_rank, const fftwf_iodim64 *howmany_dims, float *in, fftwf_complex *out, unsigned flags); +extern fftwf_plan fftwf_plan_guru64_dft_c2r(int rank, const fftwf_iodim64 *dims, int howmany_rank, const fftwf_iodim64 *howmany_dims, fftwf_complex *in, float *out, unsigned flags); +extern fftwf_plan fftwf_plan_guru64_split_dft_r2c( int rank, const fftwf_iodim64 *dims, int howmany_rank, const fftwf_iodim64 *howmany_dims, float *in, float *ro, float *io, unsigned flags); +extern fftwf_plan fftwf_plan_guru64_split_dft_c2r( int rank, const fftwf_iodim64 *dims, int howmany_rank, const fftwf_iodim64 *howmany_dims, float *ri, float *ii, float *out, unsigned flags); +extern void fftwf_execute_dft_r2c(const fftwf_plan p, float *in, fftwf_complex *out); +extern void fftwf_execute_dft_c2r(const fftwf_plan p, fftwf_complex *in, float *out); +extern void fftwf_execute_split_dft_r2c(const fftwf_plan p, float *in, float *ro, float *io); +extern void fftwf_execute_split_dft_c2r(const fftwf_plan p, float *ri, float *ii, float *out); +extern fftwf_plan fftwf_plan_many_r2r(int rank, const int *n, int howmany, float *in, const int *inembed, int istride, int idist, float *out, const int *onembed, int ostride, int odist, const fftwf_r2r_kind *kind, unsigned flags); +extern fftwf_plan fftwf_plan_r2r(int rank, const int *n, float *in, float *out, const fftwf_r2r_kind *kind, unsigned flags); +extern fftwf_plan fftwf_plan_r2r_1d(int n, float *in, float *out, fftwf_r2r_kind kind, unsigned flags); +extern fftwf_plan fftwf_plan_r2r_2d(int n0, int n1, float *in, float *out, fftwf_r2r_kind kind0, fftwf_r2r_kind kind1, unsigned flags); +extern fftwf_plan fftwf_plan_r2r_3d(int n0, int n1, int n2, float *in, float *out, fftwf_r2r_kind kind0, fftwf_r2r_kind kind1, fftwf_r2r_kind kind2, unsigned flags); +extern fftwf_plan fftwf_plan_guru_r2r(int rank, const fftwf_iodim *dims, int howmany_rank, const fftwf_iodim *howmany_dims, float *in, float *out, const fftwf_r2r_kind *kind, unsigned flags); +extern fftwf_plan fftwf_plan_guru64_r2r(int rank, const fftwf_iodim64 *dims, int howmany_rank, const fftwf_iodim64 *howmany_dims, float *in, float *out, const fftwf_r2r_kind *kind, unsigned flags); +extern void fftwf_execute_r2r(const fftwf_plan p, float *in, float *out); +extern void fftwf_destroy_plan(fftwf_plan p); +extern void fftwf_forget_wisdom(void); +extern void fftwf_cleanup(void); +extern void fftwf_set_timelimit(double t); +extern void fftwf_plan_with_nthreads(int nthreads); +extern int fftwf_init_threads(void); +extern void fftwf_cleanup_threads(void); +extern int fftwf_export_wisdom_to_filename(const char *filename); +extern void fftwf_export_wisdom_to_file(FILE *output_file); +extern char *fftwf_export_wisdom_to_string(void); +extern void fftwf_export_wisdom(fftwf_write_char_func write_char, void *data); +extern int fftwf_import_system_wisdom(void); +extern int fftwf_import_wisdom_from_filename(const char *filename); +extern int fftwf_import_wisdom_from_file(FILE *input_file); +extern int fftwf_import_wisdom_from_string(const char *input_string); +extern int fftwf_import_wisdom(fftwf_read_char_func read_char, void *data); +extern void fftwf_fprint_plan(const fftwf_plan p, FILE *output_file); +extern void fftwf_print_plan(const fftwf_plan p); +extern void *fftwf_malloc(size_t n); +extern float *fftwf_alloc_real(size_t n); +extern fftwf_complex *fftwf_alloc_complex(size_t n); +extern void fftwf_free(void *p); +extern void fftwf_flops(const fftwf_plan p, double *add, double *mul, double *fmas); +extern double fftwf_estimate_cost(const fftwf_plan p); +extern double fftwf_cost(const fftwf_plan p); +extern const char fftwf_version[]; +extern const char fftwf_cc[]; +extern const char fftwf_codelet_optim[]; +typedef long double fftwl_complex[2]; +typedef struct fftwl_plan_s *fftwl_plan; +typedef struct fftw_iodim_do_not_use_me fftwl_iodim; +typedef struct fftw_iodim64_do_not_use_me fftwl_iodim64; +typedef enum fftw_r2r_kind_do_not_use_me fftwl_r2r_kind; +typedef fftw_write_char_func_do_not_use_me fftwl_write_char_func; +typedef fftw_read_char_func_do_not_use_me fftwl_read_char_func; +extern void fftwl_execute(const fftwl_plan p); +extern fftwl_plan fftwl_plan_dft(int rank, const int *n, fftwl_complex *in, fftwl_complex *out, int sign, unsigned flags); +extern fftwl_plan fftwl_plan_dft_1d(int n, fftwl_complex *in, fftwl_complex *out, int sign, unsigned flags); +extern fftwl_plan fftwl_plan_dft_2d(int n0, int n1, fftwl_complex *in, fftwl_complex *out, int sign, unsigned flags); +extern fftwl_plan fftwl_plan_dft_3d(int n0, int n1, int n2, fftwl_complex *in, fftwl_complex *out, int sign, unsigned flags); +extern fftwl_plan fftwl_plan_many_dft(int rank, const int *n, int howmany, fftwl_complex *in, const int *inembed, int istride, int idist, fftwl_complex *out, const int *onembed, int ostride, int odist, int sign, unsigned flags); +extern fftwl_plan fftwl_plan_guru_dft(int rank, const fftwl_iodim *dims, int howmany_rank, const fftwl_iodim *howmany_dims, fftwl_complex *in, fftwl_complex *out, int sign, unsigned flags); +extern fftwl_plan fftwl_plan_guru_split_dft(int rank, const fftwl_iodim *dims, int howmany_rank, const fftwl_iodim *howmany_dims, long double *ri, long double *ii, long double *ro, long double *io, unsigned flags); +extern fftwl_plan fftwl_plan_guru64_dft(int rank, const fftwl_iodim64 *dims, int howmany_rank, const fftwl_iodim64 *howmany_dims, fftwl_complex *in, fftwl_complex *out, int sign, unsigned flags); +extern fftwl_plan fftwl_plan_guru64_split_dft(int rank, const fftwl_iodim64 *dims, int howmany_rank, const fftwl_iodim64 *howmany_dims, long double *ri, long double *ii, long double *ro, long double *io, unsigned flags); +extern void fftwl_execute_dft(const fftwl_plan p, fftwl_complex *in, fftwl_complex *out); +extern void fftwl_execute_split_dft(const fftwl_plan p, long double *ri, long double *ii, long double *ro, long double *io); +extern fftwl_plan fftwl_plan_many_dft_r2c(int rank, const int *n, int howmany, long double *in, const int *inembed, int istride, int idist, fftwl_complex *out, const int *onembed, int ostride, int odist, unsigned flags); +extern fftwl_plan fftwl_plan_dft_r2c(int rank, const int *n, long double *in, fftwl_complex *out, unsigned flags); +extern fftwl_plan fftwl_plan_dft_r2c_1d(int n,long double *in,fftwl_complex *out,unsigned flags); +extern fftwl_plan fftwl_plan_dft_r2c_2d(int n0, int n1, long double *in, fftwl_complex *out, unsigned flags); +extern fftwl_plan fftwl_plan_dft_r2c_3d(int n0, int n1, int n2, long double *in, fftwl_complex *out, unsigned flags); +extern fftwl_plan fftwl_plan_many_dft_c2r(int rank, const int *n, int howmany, fftwl_complex *in, const int *inembed, int istride, int idist, long double *out, const int *onembed, int ostride, int odist, unsigned flags); +extern fftwl_plan fftwl_plan_dft_c2r(int rank, const int *n, fftwl_complex *in, long double *out, unsigned flags); +extern fftwl_plan fftwl_plan_dft_c2r_1d(int n,fftwl_complex *in,long double *out,unsigned flags); +extern fftwl_plan fftwl_plan_dft_c2r_2d(int n0, int n1, fftwl_complex *in, long double *out, unsigned flags); +extern fftwl_plan fftwl_plan_dft_c2r_3d(int n0, int n1, int n2, fftwl_complex *in, long double *out, unsigned flags); +extern fftwl_plan fftwl_plan_guru_dft_r2c(int rank, const fftwl_iodim *dims, int howmany_rank, const fftwl_iodim *howmany_dims, long double *in, fftwl_complex *out, unsigned flags); +extern fftwl_plan fftwl_plan_guru_dft_c2r(int rank, const fftwl_iodim *dims, int howmany_rank, const fftwl_iodim *howmany_dims, fftwl_complex *in, long double *out, unsigned flags); +extern fftwl_plan fftwl_plan_guru_split_dft_r2c( int rank, const fftwl_iodim *dims, int howmany_rank, const fftwl_iodim *howmany_dims, long double *in, long double *ro, long double *io, unsigned flags); +extern fftwl_plan fftwl_plan_guru_split_dft_c2r( int rank, const fftwl_iodim *dims, int howmany_rank, const fftwl_iodim *howmany_dims, long double *ri, long double *ii, long double *out, unsigned flags); +extern fftwl_plan fftwl_plan_guru64_dft_r2c(int rank, const fftwl_iodim64 *dims, int howmany_rank, const fftwl_iodim64 *howmany_dims, long double *in, fftwl_complex *out, unsigned flags); +extern fftwl_plan fftwl_plan_guru64_dft_c2r(int rank, const fftwl_iodim64 *dims, int howmany_rank, const fftwl_iodim64 *howmany_dims, fftwl_complex *in, long double *out, unsigned flags); +extern fftwl_plan fftwl_plan_guru64_split_dft_r2c( int rank, const fftwl_iodim64 *dims, int howmany_rank, const fftwl_iodim64 *howmany_dims, long double *in, long double *ro, long double *io, unsigned flags); +extern fftwl_plan fftwl_plan_guru64_split_dft_c2r( int rank, const fftwl_iodim64 *dims, int howmany_rank, const fftwl_iodim64 *howmany_dims, long double *ri, long double *ii, long double *out, unsigned flags); +extern void fftwl_execute_dft_r2c(const fftwl_plan p, long double *in, fftwl_complex *out); +extern void fftwl_execute_dft_c2r(const fftwl_plan p, fftwl_complex *in, long double *out); +extern void fftwl_execute_split_dft_r2c(const fftwl_plan p, long double *in, long double *ro, long double *io); +extern void fftwl_execute_split_dft_c2r(const fftwl_plan p, long double *ri, long double *ii, long double *out); +extern fftwl_plan fftwl_plan_many_r2r(int rank, const int *n, int howmany, long double *in, const int *inembed, int istride, int idist, long double *out, const int *onembed, int ostride, int odist, const fftwl_r2r_kind *kind, unsigned flags); +extern fftwl_plan fftwl_plan_r2r(int rank, const int *n, long double *in, long double *out, const fftwl_r2r_kind *kind, unsigned flags); +extern fftwl_plan fftwl_plan_r2r_1d(int n, long double *in, long double *out, fftwl_r2r_kind kind, unsigned flags); +extern fftwl_plan fftwl_plan_r2r_2d(int n0, int n1, long double *in, long double *out, fftwl_r2r_kind kind0, fftwl_r2r_kind kind1, unsigned flags); +extern fftwl_plan fftwl_plan_r2r_3d(int n0, int n1, int n2, long double *in, long double *out, fftwl_r2r_kind kind0, fftwl_r2r_kind kind1, fftwl_r2r_kind kind2, unsigned flags); +extern fftwl_plan fftwl_plan_guru_r2r(int rank, const fftwl_iodim *dims, int howmany_rank, const fftwl_iodim *howmany_dims, long double *in, long double *out, const fftwl_r2r_kind *kind, unsigned flags); +extern fftwl_plan fftwl_plan_guru64_r2r(int rank, const fftwl_iodim64 *dims, int howmany_rank, const fftwl_iodim64 *howmany_dims, long double *in, long double *out, const fftwl_r2r_kind *kind, unsigned flags); +extern void fftwl_execute_r2r(const fftwl_plan p, long double *in, long double *out); +extern void fftwl_destroy_plan(fftwl_plan p); +extern void fftwl_forget_wisdom(void); +extern void fftwl_cleanup(void); +extern void fftwl_set_timelimit(double t); +extern void fftwl_plan_with_nthreads(int nthreads); +extern int fftwl_init_threads(void); +extern void fftwl_cleanup_threads(void); +extern int fftwl_export_wisdom_to_filename(const char *filename); +extern void fftwl_export_wisdom_to_file(FILE *output_file); +extern char *fftwl_export_wisdom_to_string(void); +extern void fftwl_export_wisdom(fftwl_write_char_func write_char, void *data); +extern int fftwl_import_system_wisdom(void); +extern int fftwl_import_wisdom_from_filename(const char *filename); +extern int fftwl_import_wisdom_from_file(FILE *input_file); +extern int fftwl_import_wisdom_from_string(const char *input_string); +extern int fftwl_import_wisdom(fftwl_read_char_func read_char, void *data); +extern void fftwl_fprint_plan(const fftwl_plan p, FILE *output_file); +extern void fftwl_print_plan(const fftwl_plan p); +extern void *fftwl_malloc(size_t n); +extern long double *fftwl_alloc_real(size_t n); +extern fftwl_complex *fftwl_alloc_complex(size_t n); +extern void fftwl_free(void *p); +extern void fftwl_flops(const fftwl_plan p, double *add, double *mul, double *fmas); +extern double fftwl_estimate_cost(const fftwl_plan p); +extern double fftwl_cost(const fftwl_plan p); +extern const char fftwl_version[]; +extern const char fftwl_cc[]; +extern const char fftwl_codelet_optim[]; +typedef __float128 fftwq_complex[2]; +typedef struct fftwq_plan_s *fftwq_plan; +typedef struct fftw_iodim_do_not_use_me fftwq_iodim; +typedef struct fftw_iodim64_do_not_use_me fftwq_iodim64; +typedef enum fftw_r2r_kind_do_not_use_me fftwq_r2r_kind; +typedef fftw_write_char_func_do_not_use_me fftwq_write_char_func; +typedef fftw_read_char_func_do_not_use_me fftwq_read_char_func; +extern void fftwq_execute(const fftwq_plan p); +extern fftwq_plan fftwq_plan_dft(int rank, const int *n, fftwq_complex *in, fftwq_complex *out, int sign, unsigned flags); +extern fftwq_plan fftwq_plan_dft_1d(int n, fftwq_complex *in, fftwq_complex *out, int sign, unsigned flags); +extern fftwq_plan fftwq_plan_dft_2d(int n0, int n1, fftwq_complex *in, fftwq_complex *out, int sign, unsigned flags); +extern fftwq_plan fftwq_plan_dft_3d(int n0, int n1, int n2, fftwq_complex *in, fftwq_complex *out, int sign, unsigned flags); +extern fftwq_plan fftwq_plan_many_dft(int rank, const int *n, int howmany, fftwq_complex *in, const int *inembed, int istride, int idist, fftwq_complex *out, const int *onembed, int ostride, int odist, int sign, unsigned flags); +extern fftwq_plan fftwq_plan_guru_dft(int rank, const fftwq_iodim *dims, int howmany_rank, const fftwq_iodim *howmany_dims, fftwq_complex *in, fftwq_complex *out, int sign, unsigned flags); +extern fftwq_plan fftwq_plan_guru_split_dft(int rank, const fftwq_iodim *dims, int howmany_rank, const fftwq_iodim *howmany_dims, __float128 *ri, __float128 *ii, __float128 *ro, __float128 *io, unsigned flags); +extern fftwq_plan fftwq_plan_guru64_dft(int rank, const fftwq_iodim64 *dims, int howmany_rank, const fftwq_iodim64 *howmany_dims, fftwq_complex *in, fftwq_complex *out, int sign, unsigned flags); +extern fftwq_plan fftwq_plan_guru64_split_dft(int rank, const fftwq_iodim64 *dims, int howmany_rank, const fftwq_iodim64 *howmany_dims, __float128 *ri, __float128 *ii, __float128 *ro, __float128 *io, unsigned flags); +extern void fftwq_execute_dft(const fftwq_plan p, fftwq_complex *in, fftwq_complex *out); +extern void fftwq_execute_split_dft(const fftwq_plan p, __float128 *ri, __float128 *ii, __float128 *ro, __float128 *io); +extern fftwq_plan fftwq_plan_many_dft_r2c(int rank, const int *n, int howmany, __float128 *in, const int *inembed, int istride, int idist, fftwq_complex *out, const int *onembed, int ostride, int odist, unsigned flags); +extern fftwq_plan fftwq_plan_dft_r2c(int rank, const int *n, __float128 *in, fftwq_complex *out, unsigned flags); +extern fftwq_plan fftwq_plan_dft_r2c_1d(int n,__float128 *in,fftwq_complex *out,unsigned flags); +extern fftwq_plan fftwq_plan_dft_r2c_2d(int n0, int n1, __float128 *in, fftwq_complex *out, unsigned flags); +extern fftwq_plan fftwq_plan_dft_r2c_3d(int n0, int n1, int n2, __float128 *in, fftwq_complex *out, unsigned flags); +extern fftwq_plan fftwq_plan_many_dft_c2r(int rank, const int *n, int howmany, fftwq_complex *in, const int *inembed, int istride, int idist, __float128 *out, const int *onembed, int ostride, int odist, unsigned flags); +extern fftwq_plan fftwq_plan_dft_c2r(int rank, const int *n, fftwq_complex *in, __float128 *out, unsigned flags); +extern fftwq_plan fftwq_plan_dft_c2r_1d(int n,fftwq_complex *in,__float128 *out,unsigned flags); +extern fftwq_plan fftwq_plan_dft_c2r_2d(int n0, int n1, fftwq_complex *in, __float128 *out, unsigned flags); +extern fftwq_plan fftwq_plan_dft_c2r_3d(int n0, int n1, int n2, fftwq_complex *in, __float128 *out, unsigned flags); +extern fftwq_plan fftwq_plan_guru_dft_r2c(int rank, const fftwq_iodim *dims, int howmany_rank, const fftwq_iodim *howmany_dims, __float128 *in, fftwq_complex *out, unsigned flags); +extern fftwq_plan fftwq_plan_guru_dft_c2r(int rank, const fftwq_iodim *dims, int howmany_rank, const fftwq_iodim *howmany_dims, fftwq_complex *in, __float128 *out, unsigned flags); +extern fftwq_plan fftwq_plan_guru_split_dft_r2c( int rank, const fftwq_iodim *dims, int howmany_rank, const fftwq_iodim *howmany_dims, __float128 *in, __float128 *ro, __float128 *io, unsigned flags); +extern fftwq_plan fftwq_plan_guru_split_dft_c2r( int rank, const fftwq_iodim *dims, int howmany_rank, const fftwq_iodim *howmany_dims, __float128 *ri, __float128 *ii, __float128 *out, unsigned flags); +extern fftwq_plan fftwq_plan_guru64_dft_r2c(int rank, const fftwq_iodim64 *dims, int howmany_rank, const fftwq_iodim64 *howmany_dims, __float128 *in, fftwq_complex *out, unsigned flags); +extern fftwq_plan fftwq_plan_guru64_dft_c2r(int rank, const fftwq_iodim64 *dims, int howmany_rank, const fftwq_iodim64 *howmany_dims, fftwq_complex *in, __float128 *out, unsigned flags); +extern fftwq_plan fftwq_plan_guru64_split_dft_r2c( int rank, const fftwq_iodim64 *dims, int howmany_rank, const fftwq_iodim64 *howmany_dims, __float128 *in, __float128 *ro, __float128 *io, unsigned flags); +extern fftwq_plan fftwq_plan_guru64_split_dft_c2r( int rank, const fftwq_iodim64 *dims, int howmany_rank, const fftwq_iodim64 *howmany_dims, __float128 *ri, __float128 *ii, __float128 *out, unsigned flags); +extern void fftwq_execute_dft_r2c(const fftwq_plan p, __float128 *in, fftwq_complex *out); +extern void fftwq_execute_dft_c2r(const fftwq_plan p, fftwq_complex *in, __float128 *out); +extern void fftwq_execute_split_dft_r2c(const fftwq_plan p, __float128 *in, __float128 *ro, __float128 *io); +extern void fftwq_execute_split_dft_c2r(const fftwq_plan p, __float128 *ri, __float128 *ii, __float128 *out); +extern fftwq_plan fftwq_plan_many_r2r(int rank, const int *n, int howmany, __float128 *in, const int *inembed, int istride, int idist, __float128 *out, const int *onembed, int ostride, int odist, const fftwq_r2r_kind *kind, unsigned flags); +extern fftwq_plan fftwq_plan_r2r(int rank, const int *n, __float128 *in, __float128 *out, const fftwq_r2r_kind *kind, unsigned flags); +extern fftwq_plan fftwq_plan_r2r_1d(int n, __float128 *in, __float128 *out, fftwq_r2r_kind kind, unsigned flags); +extern fftwq_plan fftwq_plan_r2r_2d(int n0, int n1, __float128 *in, __float128 *out, fftwq_r2r_kind kind0, fftwq_r2r_kind kind1, unsigned flags); +extern fftwq_plan fftwq_plan_r2r_3d(int n0, int n1, int n2, __float128 *in, __float128 *out, fftwq_r2r_kind kind0, fftwq_r2r_kind kind1, fftwq_r2r_kind kind2, unsigned flags); +extern fftwq_plan fftwq_plan_guru_r2r(int rank, const fftwq_iodim *dims, int howmany_rank, const fftwq_iodim *howmany_dims, __float128 *in, __float128 *out, const fftwq_r2r_kind *kind, unsigned flags); +extern fftwq_plan fftwq_plan_guru64_r2r(int rank, const fftwq_iodim64 *dims, int howmany_rank, const fftwq_iodim64 *howmany_dims, __float128 *in, __float128 *out, const fftwq_r2r_kind *kind, unsigned flags); +extern void fftwq_execute_r2r(const fftwq_plan p, __float128 *in, __float128 *out); +extern void fftwq_destroy_plan(fftwq_plan p); +extern void fftwq_forget_wisdom(void); +extern void fftwq_cleanup(void); +extern void fftwq_set_timelimit(double t); +extern void fftwq_plan_with_nthreads(int nthreads); +extern int fftwq_init_threads(void); +extern void fftwq_cleanup_threads(void); +extern int fftwq_export_wisdom_to_filename(const char *filename); +extern void fftwq_export_wisdom_to_file(FILE *output_file); +extern char *fftwq_export_wisdom_to_string(void); +extern void fftwq_export_wisdom(fftwq_write_char_func write_char, void *data); +extern int fftwq_import_system_wisdom(void); +extern int fftwq_import_wisdom_from_filename(const char *filename); +extern int fftwq_import_wisdom_from_file(FILE *input_file); +extern int fftwq_import_wisdom_from_string(const char *input_string); +extern int fftwq_import_wisdom(fftwq_read_char_func read_char, void *data); +extern void fftwq_fprint_plan(const fftwq_plan p, FILE *output_file); +extern void fftwq_print_plan(const fftwq_plan p); +extern void *fftwq_malloc(size_t n); +extern __float128 *fftwq_alloc_real(size_t n); +extern fftwq_complex *fftwq_alloc_complex(size_t n); +extern void fftwq_free(void *p); +extern void fftwq_flops(const fftwq_plan p, double *add, double *mul, double *fmas); +extern double fftwq_estimate_cost(const fftwq_plan p); +extern double fftwq_cost(const fftwq_plan p); +extern const char fftwq_version[]; +extern const char fftwq_cc[]; +extern const char fftwq_codelet_optim[]; + +]] diff --git a/fftw3/defines.lua b/fftw3/defines.lua new file mode 100644 index 000000000..a252b2aee --- /dev/null +++ b/fftw3/defines.lua @@ -0,0 +1,36 @@ +local ffi=require 'ffi' +local defines = {} +function defines.register_hashdefs(fftw, C) + fftw.MEASURE = 0 + fftw.DESTROY_INPUT = 1 + fftw.UNALIGNED = 2 + fftw.CONSERVE_MEMORY = 4 + fftw.EXHAUSTIVE = 8 -- /* NO_EXHAUSTIVE is default */ + fftw.PRESERVE_INPUT = 16 -- /* cancels DESTROY_INPUT */ + fftw.PATIENT =32 -- /* IMPATIENT is default */ + fftw.ESTIMATE = 64 + fftw.WISDOM_ONLY = 2097152 -- (1U << 21) + + fftw.FORWARD = -1 + fftw.BACKWARD = 1 + + fftw.NO_TIMELIMIT = (-1.0) + + -- -- typedef + fftw.r2r_kind = ffi.typeof('fftw_r2r_kind') + fftw.R2HC=0 + fftw.HC2R=1 + fftw.DHT=2 + + fftw.REDFT00=3 + fftw.REDFT01=4 + fftw.REDFT10=5 + fftw.REDFT11=6 + + fftw.RODFT00=7 + fftw.RODFT01=8 + fftw.RODFT10=9 + fftw.RODFT11=10 +end + +return defines \ No newline at end of file diff --git a/fftw3/init.lua b/fftw3/init.lua new file mode 100644 index 000000000..5e5e5590d --- /dev/null +++ b/fftw3/init.lua @@ -0,0 +1,249 @@ +-- Do not change this file manually +-- Generated with dev/create-init.lua + +local ffi = require 'ffi' +local C = ffi.load('fftw3') +local fftw = {C=C} + +require 'fftw3.cdefs' + +local defines = require 'fftw3.defines' +defines.register_hashdefs(fftw, C) + +local function register(luafuncname, funcname) + local symexists, msg = pcall(function() + local sym = C[funcname] + end) + if symexists then + fftw[luafuncname] = C[funcname] + end +end + +register('execute', 'fftw_execute') +register('plan_dft', 'fftw_plan_dft') +register('plan_dft_1d', 'fftw_plan_dft_1d') +register('plan_dft_2d', 'fftw_plan_dft_2d') +register('plan_dft_3d', 'fftw_plan_dft_3d') +register('plan_many_dft', 'fftw_plan_many_dft') +register('plan_guru_dft', 'fftw_plan_guru_dft') +register('plan_guru_split_dft', 'fftw_plan_guru_split_dft') +register('plan_guru64_dft', 'fftw_plan_guru64_dft') +register('plan_guru64_split_dft', 'fftw_plan_guru64_split_dft') +register('execute_dft', 'fftw_execute_dft') +register('execute_split_dft', 'fftw_execute_split_dft') +register('plan_many_dft_r2c', 'fftw_plan_many_dft_r2c') +register('plan_dft_r2c', 'fftw_plan_dft_r2c') +register('plan_dft_r2c_1d', 'fftw_plan_dft_r2c_1d') +register('plan_dft_r2c_2d', 'fftw_plan_dft_r2c_2d') +register('plan_dft_r2c_3d', 'fftw_plan_dft_r2c_3d') +register('plan_many_dft_c2r', 'fftw_plan_many_dft_c2r') +register('plan_dft_c2r', 'fftw_plan_dft_c2r') +register('plan_dft_c2r_1d', 'fftw_plan_dft_c2r_1d') +register('plan_dft_c2r_2d', 'fftw_plan_dft_c2r_2d') +register('plan_dft_c2r_3d', 'fftw_plan_dft_c2r_3d') +register('plan_guru_dft_r2c', 'fftw_plan_guru_dft_r2c') +register('plan_guru_dft_c2r', 'fftw_plan_guru_dft_c2r') +register('plan_guru_split_dft_r2c', 'fftw_plan_guru_split_dft_r2c') +register('plan_guru_split_dft_c2r', 'fftw_plan_guru_split_dft_c2r') +register('plan_guru64_dft_r2c', 'fftw_plan_guru64_dft_r2c') +register('plan_guru64_dft_c2r', 'fftw_plan_guru64_dft_c2r') +register('plan_guru64_split_dft_r2c', 'fftw_plan_guru64_split_dft_r2c') +register('plan_guru64_split_dft_c2r', 'fftw_plan_guru64_split_dft_c2r') +register('execute_dft_r2c', 'fftw_execute_dft_r2c') +register('execute_dft_c2r', 'fftw_execute_dft_c2r') +register('execute_split_dft_r2c', 'fftw_execute_split_dft_r2c') +register('execute_split_dft_c2r', 'fftw_execute_split_dft_c2r') +register('plan_many_r2r', 'fftw_plan_many_r2r') +register('plan_r2r', 'fftw_plan_r2r') +register('plan_r2r_1d', 'fftw_plan_r2r_1d') +register('plan_r2r_2d', 'fftw_plan_r2r_2d') +register('plan_r2r_3d', 'fftw_plan_r2r_3d') +register('plan_guru_r2r', 'fftw_plan_guru_r2r') +register('plan_guru64_r2r', 'fftw_plan_guru64_r2r') +register('execute_r2r', 'fftw_execute_r2r') +register('destroy_plan', 'fftw_destroy_plan') +register('forget_wisdom', 'fftw_forget_wisdom') +register('cleanup', 'fftw_cleanup') +register('set_timelimit', 'fftw_set_timelimit') +register('plan_with_nthreads', 'fftw_plan_with_nthreads') +register('init_threads', 'fftw_init_threads') +register('cleanup_threads', 'fftw_cleanup_threads') +register('export_wisdom_to_filename', 'fftw_export_wisdom_to_filename') +register('export_wisdom_to_file', 'fftw_export_wisdom_to_file') +register('export_wisdom_to_string', 'fftw_export_wisdom_to_string') +register('export_wisdom', 'fftw_export_wisdom') +register('import_system_wisdom', 'fftw_import_system_wisdom') +register('import_wisdom_from_filename', 'fftw_import_wisdom_from_filename') +register('import_wisdom_from_file', 'fftw_import_wisdom_from_file') +register('import_wisdom_from_string', 'fftw_import_wisdom_from_string') +register('import_wisdom', 'fftw_import_wisdom') +register('fprint_plan', 'fftw_fprint_plan') +register('print_plan', 'fftw_print_plan') +register('malloc', 'fftw_malloc') +register('alloc_real', 'fftw_alloc_real') +register('alloc_complex', 'fftw_alloc_complex') +register('free', 'fftw_free') +register('flops', 'fftw_flops') +register('estimate_cost', 'fftw_estimate_cost') +register('cost', 'fftw_cost') + +register('r2r_kind_do_not_use_me', 'fftw_r2r_kind_do_not_use_me') +register('iodim_do_not_use_me', 'fftw_iodim_do_not_use_me') +register('iodim64_do_not_use_me', 'fftw_iodim64_do_not_use_me') +register('write_char_func_do_not_use_me', 'fftw_write_char_func_do_not_use_me') +register('read_char_func_do_not_use_me', 'fftw_read_char_func_do_not_use_me') +register('complex', 'fftw_complex') +register('plan_s', 'fftw_plan_s') +register('plan', 'fftw_plan') +register('iodim_do_not_use_me', 'fftw_iodim_do_not_use_me') +register('iodim', 'fftw_iodim') +register('iodim64_do_not_use_me', 'fftw_iodim64_do_not_use_me') +register('iodim64', 'fftw_iodim64') +register('r2r_kind_do_not_use_me', 'fftw_r2r_kind_do_not_use_me') +register('r2r_kind', 'fftw_r2r_kind') +register('write_char_func_do_not_use_me', 'fftw_write_char_func_do_not_use_me') +register('write_char_func', 'fftw_write_char_func') +register('read_char_func_do_not_use_me', 'fftw_read_char_func_do_not_use_me') +register('read_char_func', 'fftw_read_char_func') +register('plan', 'fftw_plan') +register('plan', 'fftw_plan') +register('complex', 'fftw_complex') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('complex', 'fftw_complex') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('complex', 'fftw_complex') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('complex', 'fftw_complex') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('complex', 'fftw_complex') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('iodim', 'fftw_iodim') +register('iodim', 'fftw_iodim') +register('complex', 'fftw_complex') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('iodim', 'fftw_iodim') +register('iodim', 'fftw_iodim') +register('plan', 'fftw_plan') +register('iodim64', 'fftw_iodim64') +register('iodim64', 'fftw_iodim64') +register('complex', 'fftw_complex') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('iodim64', 'fftw_iodim64') +register('iodim64', 'fftw_iodim64') +register('plan', 'fftw_plan') +register('complex', 'fftw_complex') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('plan', 'fftw_plan') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('iodim', 'fftw_iodim') +register('iodim', 'fftw_iodim') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('iodim', 'fftw_iodim') +register('iodim', 'fftw_iodim') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('iodim', 'fftw_iodim') +register('iodim', 'fftw_iodim') +register('plan', 'fftw_plan') +register('iodim', 'fftw_iodim') +register('iodim', 'fftw_iodim') +register('plan', 'fftw_plan') +register('iodim64', 'fftw_iodim64') +register('iodim64', 'fftw_iodim64') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('iodim64', 'fftw_iodim64') +register('iodim64', 'fftw_iodim64') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('iodim64', 'fftw_iodim64') +register('iodim64', 'fftw_iodim64') +register('plan', 'fftw_plan') +register('iodim64', 'fftw_iodim64') +register('iodim64', 'fftw_iodim64') +register('plan', 'fftw_plan') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('plan', 'fftw_plan') +register('plan', 'fftw_plan') +register('r2r_kind', 'fftw_r2r_kind') +register('plan', 'fftw_plan') +register('r2r_kind', 'fftw_r2r_kind') +register('plan', 'fftw_plan') +register('r2r_kind', 'fftw_r2r_kind') +register('plan', 'fftw_plan') +register('r2r_kind', 'fftw_r2r_kind') +register('r2r_kind', 'fftw_r2r_kind') +register('plan', 'fftw_plan') +register('r2r_kind', 'fftw_r2r_kind') +register('r2r_kind', 'fftw_r2r_kind') +register('r2r_kind', 'fftw_r2r_kind') +register('plan', 'fftw_plan') +register('iodim', 'fftw_iodim') +register('iodim', 'fftw_iodim') +register('r2r_kind', 'fftw_r2r_kind') +register('plan', 'fftw_plan') +register('iodim64', 'fftw_iodim64') +register('iodim64', 'fftw_iodim64') +register('r2r_kind', 'fftw_r2r_kind') +register('plan', 'fftw_plan') +register('plan', 'fftw_plan') +register('write_char_func', 'fftw_write_char_func') +register('read_char_func', 'fftw_read_char_func') +register('plan', 'fftw_plan') +register('plan', 'fftw_plan') +register('complex', 'fftw_complex') +register('plan', 'fftw_plan') +register('plan', 'fftw_plan') +register('plan', 'fftw_plan') +register('version', 'fftw_version') +register('cc', 'fftw_cc') +register('codelet_optim', 'fftw_codelet_optim') +register('iodim_do_not_use_me', 'fftw_iodim_do_not_use_me') +register('iodim64_do_not_use_me', 'fftw_iodim64_do_not_use_me') +register('r2r_kind_do_not_use_me', 'fftw_r2r_kind_do_not_use_me') +register('write_char_func_do_not_use_me', 'fftw_write_char_func_do_not_use_me') +register('read_char_func_do_not_use_me', 'fftw_read_char_func_do_not_use_me') +register('iodim_do_not_use_me', 'fftw_iodim_do_not_use_me') +register('iodim64_do_not_use_me', 'fftw_iodim64_do_not_use_me') +register('r2r_kind_do_not_use_me', 'fftw_r2r_kind_do_not_use_me') +register('write_char_func_do_not_use_me', 'fftw_write_char_func_do_not_use_me') +register('read_char_func_do_not_use_me', 'fftw_read_char_func_do_not_use_me') +register('iodim_do_not_use_me', 'fftw_iodim_do_not_use_me') +register('iodim64_do_not_use_me', 'fftw_iodim64_do_not_use_me') +register('r2r_kind_do_not_use_me', 'fftw_r2r_kind_do_not_use_me') +register('write_char_func_do_not_use_me', 'fftw_write_char_func_do_not_use_me') +register('read_char_func_do_not_use_me', 'fftw_read_char_func_do_not_use_me') + +return fftw + From 45b20a32bacf0f2061c0ef6483d077ec35589d69 Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Thu, 6 Feb 2014 18:14:12 +0100 Subject: [PATCH 04/20] Removed old fft implementation. --- fft-init.lua | 306 --------------------------------------------------- gslext.lua | 2 +- 2 files changed, 1 insertion(+), 307 deletions(-) delete mode 100644 fft-init.lua diff --git a/fft-init.lua b/fft-init.lua deleted file mode 100644 index d7bad58b3..000000000 --- a/fft-init.lua +++ /dev/null @@ -1,306 +0,0 @@ - -local ffi = require 'ffi' -local bit = require 'bit' -local gsl = require 'gsl' - -local gsl_check = require 'gsl-check' - -local check = require 'check' -local is_integer = check.is_integer - -local tobit, band, rshift = bit.tobit, bit.band, bit.rshift -local tonumber = tonumber - -ffi.cdef [[ - typedef struct - { - size_t size; - size_t stride; - double * data; - gsl_block * block; - } fft_hc; - - typedef struct - { - size_t size; - size_t stride; - double * data; - gsl_block * block; - } fft_radix2_hc; - - typedef struct - { - size_t n; - size_t nf; - size_t factor[64]; - gsl_complex *twiddle[64]; - gsl_complex *trig; - } gsl_fft_real_wavetable; - - typedef struct - { - size_t n; - double *scratch; - } gsl_fft_real_workspace; - - typedef struct - { - size_t n; - size_t nf; - size_t factor[64]; - gsl_complex *twiddle[64]; - gsl_complex *trig; - } gsl_fft_halfcomplex_wavetable; - - gsl_fft_real_wavetable * gsl_fft_real_wavetable_alloc (size_t n); - - void gsl_fft_real_wavetable_free (gsl_fft_real_wavetable * wavetable); - - gsl_fft_halfcomplex_wavetable * gsl_fft_halfcomplex_wavetable_alloc (size_t n); - - void gsl_fft_halfcomplex_wavetable_free (gsl_fft_halfcomplex_wavetable * wavetable); - - gsl_fft_real_workspace * gsl_fft_real_workspace_alloc (size_t n); - - void gsl_fft_real_workspace_free (gsl_fft_real_workspace * workspace); - - int gsl_fft_real_radix2_transform (double data[], const size_t stride, - const size_t n) ; - - int gsl_fft_halfcomplex_radix2_inverse (double data[], - size_t stride, size_t n); - - int gsl_fft_real_transform (double data[], const size_t stride, const size_t n, - const gsl_fft_real_wavetable * wavetable, - gsl_fft_real_workspace * work); - - int gsl_fft_halfcomplex_inverse (double data[], const size_t stride, const size_t n, - const gsl_fft_halfcomplex_wavetable * wavetable, - gsl_fft_real_workspace * work); - - int gsl_fft_halfcomplex_transform (double data[], const size_t stride, const size_t n, - const gsl_fft_halfcomplex_wavetable * wavetable, - gsl_fft_real_workspace * work); - -]] - -local fft_hc = ffi.typeof('fft_hc') -local fft_radix2_hc = ffi.typeof('fft_radix2_hc') -local gsl_matrix = ffi.typeof('gsl_matrix') - -local function is_two_power(n) - if n > 0 then - local k = tobit(n) - while band(k, 1) == 0 do k = rshift(k, 1) end - return (k == 1) - end -end - -local cache_n = {} -local cache_r = {} - -local function res_allocator(name) - local alloc = gsl['gsl_fft_' .. name .. '_alloc'] - local free = gsl['gsl_fft_' .. name .. '_free'] - return function(n) - return ffi.gc(alloc(n), free) - end -end - -local cache_allocator = { - real_wavetable = res_allocator('real_wavetable'), - halfcomplex_wavetable = res_allocator('halfcomplex_wavetable'), - real_workspace = res_allocator('real_workspace') -} - -local function get_resource(name, n) - local resource - if cache_n[name] ~= n then - resource = cache_allocator[name](n) - cache_n[name] = n - cache_r[name] = resource - else - resource = cache_r[name] - end - return resource -end - -local function get_matrix_block(x, ip) - local n = tonumber(x.size1) - local b, data, stride - if ip then - b, data, stride = x.block, x.data, x.tda - b.ref_count = b.ref_count + 1 - else - b = matrix.block(n) - data, stride = b.data, 1 - for i=0, n-1 do data[i] = x.data[x.tda * i] end - end - return b, data, stride -end - -local function get_hc_block(ft, ip) - local n = tonumber(ft.size) - local b, data, stride - if ip then - b, data, stride = ft.block, ft.data, tonumber(ft.stride) - b.ref_count = b.ref_count + 1 - else - b = matrix.block(n) - data, stride = b.data, 1 - for i=0, n-1 do data[i] = ft.data[tonumber(ft.stride) * i] end - end - return b, data, stride -end - -function num.fft(x, ip) - local n = tonumber(x.size1) - local b, data, stride = get_matrix_block(x, ip) - if is_two_power(n) then - gsl_check(gsl.gsl_fft_real_radix2_transform(data, stride, n)) - return fft_radix2_hc(n, stride, data, b) - else - local wt = get_resource('real_wavetable', n) - local ws = get_resource('real_workspace', n) - gsl_check(gsl.gsl_fft_real_transform(data, stride, n, wt, ws)) - return fft_hc(n, stride, data, b) - end -end - -function num.fftinv(ft, ip) - local n = tonumber(ft.size) - local b, data, stride = get_hc_block(ft, ip) - if is_two_power(n) then - gsl_check(gsl.gsl_fft_halfcomplex_radix2_inverse(data, stride, n)) - else - local wt = get_resource('halfcomplex_wavetable', n) - local ws = get_resource('real_workspace', n) - gsl_check(gsl.gsl_fft_halfcomplex_inverse(data, stride, n, wt, ws)) - end - return gsl_matrix(n, 1, stride, data, b, 1) -end - -local function halfcomplex_radix2_index(n, stride, k) - if k < 0 or k >= n then error('invalid halfcomplex index', 2) end - local half_n = n/2 - if k == 0 then - return 0, 0 - elseif k < half_n then - return 1, k, n-k - elseif k == half_n then - return 0, half_n - elseif k > half_n then - return -1, n-k, k - end -end - -local function halfcomplex_index(n, stride, k) - if k < 0 or k >= n then error('invalid halfcomplex index', 2) end - local half_n = n/2 - if k == 0 then - return 0, 0 - elseif k < half_n then - return 1, 2*k-1, 2*k - elseif k == half_n then - return 0, half_n - elseif k > half_n then - return -1, 2*(n-k)-1, 2*(n-k) - end -end - -local function halfcomplex_get(indexer, data, n, stride, k) - local isign, ridx, iidx = indexer(n, stride, k) - local r = data[stride*ridx] - local i = (isign == 0 and 0 or isign * data[stride*iidx]) - return complex.new(r, i) -end - -local function halfcomplex_set(indexer, data, n, stride, k, z) - local isign, ridx, iidx = indexer(n, stride, k) - local r, i = complex.rect(z) - data[stride*ridx] = r - if isign ~= 0 then - data[stride*iidx] = isign * i - end -end - -local function hc_length(ft) - return tonumber(ft.size) -end - -local function halfcomplex_to_matrix(hc) - return matrix.cnew(tonumber(hc.size), 1, function(i) return hc[i-1] end) -end - -local function hc_tostring(hc) - local m = halfcomplex_to_matrix(hc) - return m:show() -end - -local function hc_radix2_index(ft, k) - if is_integer(k) then - local idx = halfcomplex_radix2_index - local size, stride = tonumber(ft.size), tonumber(ft.stride) - return halfcomplex_get(idx, ft.data, size, stride, k) - elseif k == 'show' then - return hc_tostring - end -end - -local function hc_radix2_newindex(ft, k, z) - if is_integer(k) then - local idx = halfcomplex_radix2_index - local size, stride = tonumber(ft.size), tonumber(ft.stride) - return halfcomplex_set(idx, ft.data, size, stride, k, z) - end -end - -local function hc_index(ft, k) - if is_integer(k) then - local idx = halfcomplex_index - local size, stride = tonumber(ft.size), tonumber(ft.stride) - return halfcomplex_get(idx, ft.data, size, stride, k) - elseif k == 'show' then - return hc_tostring - end -end - -local function hc_newindex(ft, k, z) - if is_integer(k) then - local idx = halfcomplex_index - local size, stride = tonumber(ft.size), tonumber(ft.stride) - return halfcomplex_set(idx, ft.data, size, stride, k, z) - end -end - -local function hc_free(hc) - local b = hc.block - b.ref_count = b.ref_count - 1 - if b.ref_count == 0 then - ffi.C.free(b.data) - ffi.C.free(b) - end -end - -ffi.metatype(fft_hc, { - __gc = hc_free, - __index = hc_index, - __newindex = hc_newindex, - __len = hc_length, --- __tostring = hc_tostring, - } - ) - -ffi.metatype(fft_radix2_hc, { - __gc = hc_free, - __index = hc_radix2_index, - __newindex = hc_radix2_newindex, - __len = hc_length, --- __tostring = hc_tostring, - } - ) - -local register_ffi_type = debug.getregistry().__gsl_reg_ffi_type - -register_ffi_type(fft_radix2_hc, "radix2 half-complex vector") -register_ffi_type(fft_hc, "half-complex vector") diff --git a/gslext.lua b/gslext.lua index 957482dfd..98936fa5b 100644 --- a/gslext.lua +++ b/gslext.lua @@ -7,7 +7,7 @@ require('num') require('rng') require('rnd') require('integ-init') -require('fft-init') +require('fft') require('graph-init') require('randist') require('import') From eb5dda358fa540a556bddbb88ce5d0404cbcaae6 Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Thu, 6 Feb 2014 19:30:37 +0100 Subject: [PATCH 05/20] Changed fft naming. --- fft.lua | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/fft.lua b/fft.lua index 2b9ecc6fa..737e5d704 100644 --- a/fft.lua +++ b/fft.lua @@ -1,4 +1,5 @@ fftw = require'fftw3.init' +ffi = require'ffi' local gsl_matrix = ffi.typeof('gsl_matrix') local gsl_matrix_complex = ffi.typeof('gsl_matrix_complex') @@ -12,12 +13,12 @@ function num.fft(vec) local input = ffi.cast("fftw_complex*",vec.data) local plan = fftw.plan_dft_1d(n, input, output, fftw.FORWARD, fftw.MEASURE) fftw.execute(plan) - fftw.destroy_plan(plan) + --fftw.destroy_plan(plan) return outvec end end -function num.fft_inv(vec) +function num.fftinv(vec) if ffi.istype(gsl_matrix, vec) then return num.rfft_inv(vec) else @@ -39,14 +40,14 @@ function num.rfft(vec) local output = ffi.cast("fftw_complex*",outvec.data) local plan = fftw.plan_dft_r2c_1d(n, vec.data, output, fftw.MEASURE ) fftw.execute(plan) - fftw.destroy_plan(plan) + --fftw.destroy_plan(plan) return outvec else error("Cannot handle nonreal input data.") end end -function num.rfft_inv(vec) +function num.rfftinv(vec) if ffi.istype(gsl_matrix_complex, vec) then local n = tonumber(vec.size1) local nnew = (n-1)*2 @@ -79,7 +80,7 @@ function num.fft2(mat) end end -function num.fft2_inv(mat) +function num.fft2inv(mat) if ffi.istype(gsl_matrix, mat) then return num.rfft2_inv(mat) else @@ -110,7 +111,7 @@ function num.rfft2(mat) end end -function num.rfft2_inv(mat) +function num.rfft2inv(mat) if ffi.istype(gsl_matrix_complex, mat) then local n1 = tonumber(mat.size1) local n2 = tonumber(mat.size2) @@ -142,7 +143,7 @@ function num.fft3(datavec, n1, n2, n3) end end -function num.fft3_inv(datavec, n1, n2, n3) +function num.fft3inv(datavec, n1, n2, n3) if ffi.istype(gsl_matrix, datavec) then return num.rfft3_inv(datavec, n1, n2, n3) else @@ -169,7 +170,7 @@ function num.rfft3(datavec, n1, n2, n3) end end -function num.rfft3_inv(datavec, n1, n2, n3) +function num.rfft3inv(datavec, n1, n2, n3) if ffi.istype(gsl_matrix, datavec) then local outvec = matrix.new(n1*n2*n3, 1) local input = ffi.cast("fftw_complex*",datavec.data) @@ -198,7 +199,7 @@ function num.fftn(datavec, dimensionvec) end end -function num.fftn_inv(datavec, dimensionvec) +function num.fftninv(datavec, dimensionvec) if ffi.istype(gsl_matrix, datavec) then return num.rfftn_inv(datavec, dimensionvec) else @@ -232,7 +233,7 @@ function num.rfftn(datavec, dimensionvec) end end -function num.rfftn_inv(datavec, dimensionvec) +function num.rfftninv(datavec, dimensionvec) if ffi.istype(gsl_matrix, datavec) then local d = #dimensionvec local newsize = 1 @@ -249,4 +250,4 @@ function num.rfftn_inv(datavec, dimensionvec) else error("Cannot handle nonreal input data.") end -end +end \ No newline at end of file From 478e835df7826064944a072730c87d4d1c614e33 Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Thu, 6 Feb 2014 22:41:41 +0100 Subject: [PATCH 06/20] Added plans for optimized FFTs --- fft.lua | 212 ++++++++++++++++++++++---------------------------------- 1 file changed, 82 insertions(+), 130 deletions(-) diff --git a/fft.lua b/fft.lua index 737e5d704..586da958c 100644 --- a/fft.lua +++ b/fft.lua @@ -1,49 +1,48 @@ -fftw = require'fftw3.init' ffi = require'ffi' +fftw = require'fftw3.init' local gsl_matrix = ffi.typeof('gsl_matrix') local gsl_matrix_complex = ffi.typeof('gsl_matrix_complex') function num.fft(vec) - if ffi.istype(gsl_matrix, vec) then - return num.rfft(vec) - else + if ffi.istype(gsl_matrix_complex, vec) then local n = tonumber(vec.size1) - local outvec = matrix.cnew(n, 1) + local outvec = matrix.calloc(n, 1) local output = ffi.cast("fftw_complex*",outvec.data) local input = ffi.cast("fftw_complex*",vec.data) - local plan = fftw.plan_dft_1d(n, input, output, fftw.FORWARD, fftw.MEASURE) + + local plan = ffi.gc(fftw.plan_dft_1d(n, input, output, fftw.FORWARD, fftw.MEASURE), fftw.destroy_plan) fftw.execute(plan) - --fftw.destroy_plan(plan) - return outvec + return outvec, plan + else + error("Input must be complex valued.") end end function num.fftinv(vec) - if ffi.istype(gsl_matrix, vec) then - return num.rfft_inv(vec) - else + if ffi.istype(gsl_matrix_complex, vec) then local n = tonumber(vec.size1) - local outvec = matrix.cnew(n, 1) + local outvec = matrix.calloc(n, 1) local output = ffi.cast("fftw_complex*",outvec.data) local input = ffi.cast("fftw_complex*",vec.data) - local plan = fftw.plan_dft_1d(n, input, output, fftw.BACKWARD, fftw.MEASURE) + + local plan = ffi.gc(fftw.plan_dft_1d(n, input, output, fftw.BACKWARD, fftw.MEASURE), fftw.destroy_plan) fftw.execute(plan) - fftw.destroy_plan(plan) - return outvec/n + return outvec, plan + else + error("Input must be complex valued.") end end function num.rfft(vec) if ffi.istype(gsl_matrix, vec) then local n = tonumber(vec.size1) - local outvec = matrix.cnew(math.floor(n/2)+1, 1) + local outvec = matrix.calloc(math.floor(n/2)+1, 1) local output = ffi.cast("fftw_complex*",outvec.data) local plan = fftw.plan_dft_r2c_1d(n, vec.data, output, fftw.MEASURE ) fftw.execute(plan) - --fftw.destroy_plan(plan) - return outvec + return outvec, plan else - error("Cannot handle nonreal input data.") + error("Input must be real valued.") end end @@ -51,48 +50,45 @@ function num.rfftinv(vec) if ffi.istype(gsl_matrix_complex, vec) then local n = tonumber(vec.size1) local nnew = (n-1)*2 - local outvec = matrix.new(nnew, 1) + local outvec = matrix.alloc(nnew, 1) local input = ffi.cast("fftw_complex*",vec.data) local plan = fftw.plan_dft_c2r_1d(nnew, input, outvec.data, fftw.MEASURE ) fftw.execute(plan) - fftw.destroy_plan(plan) - return outvec/nnew + return outvec,plan else - error("Cannot handle noncomplex input data.") + error("Input must be complex valued.") end end ---------------------------------- function num.fft2(mat) - if ffi.istype(gsl_matrix, mat) then - return num.rfft2(mat) - else + if ffi.istype(gsl_matrix_complex, mat) then local n1 = tonumber(mat.size1) local n2 = tonumber(mat.size2) - local outmat = matrix.cnew(n1, n2) + local outmat = matrix.calloc(n1, n2) local output = ffi.cast("fftw_complex*",outmat.data) local input = ffi.cast("fftw_complex*", mat.data) local plan = fftw.plan_dft_2d(n1,n2, input, output, fftw.FORWARD, fftw.MEASURE) fftw.execute(plan) - fftw.destroy_plan(plan) - return outmat + return outmat, plan + else + error("Input must be complex valued.") end end function num.fft2inv(mat) - if ffi.istype(gsl_matrix, mat) then - return num.rfft2_inv(mat) - else + if ffi.istype(gsl_matrix_complex, mat) then local n1 = tonumber(mat.size1) local n2 = tonumber(mat.size2) - local outmat = matrix.cnew(n1, n2) + local outmat = matrix.calloc(n1, n2) local output = ffi.cast("fftw_complex*",outmat.data) local input = ffi.cast("fftw_complex*", mat.data) local plan = fftw.plan_dft_2d(n1,n2, input, output, fftw.BACKWARD, fftw.MEASURE) fftw.execute(plan) - fftw.destroy_plan(plan) - return outmat + return outmat, plan + else + error("Input must be complex valued.") end end @@ -100,14 +96,13 @@ function num.rfft2(mat) if ffi.istype(gsl_matrix, mat) then local n1 = tonumber(mat.size1) local n2 = tonumber(mat.size2) - local outmat = matrix.cnew(n1, math.floor(n2/2)+1) + local outmat = matrix.calloc(n1, math.floor(n2/2)+1) local output = ffi.cast("fftw_complex*",outmat.data) local plan = fftw.plan_dft_r2c_2d(n1,n2, mat.data, output,fftw.MEASURE) fftw.execute(plan) - fftw.destroy_plan(plan) - return outmat + return outmat, plan else - error("Cannot handle nonreal input matrix.") + error("Input must be real valued.") end end @@ -116,138 +111,95 @@ function num.rfft2inv(mat) local n1 = tonumber(mat.size1) local n2 = tonumber(mat.size2) local n2new = (n2-1)*2 - local outmat = matrix.new(n1, n2new) + local outmat = matrix.alloc(n1, n2new) local input = ffi.cast("fftw_complex*",mat.data) local plan = fftw.plan_dft_c2r_2d(n1,n2new, input, outmat.data,fftw.MEASURE) fftw.execute(plan) - fftw.destroy_plan(plan) - return outmat/(n1*n2new) + return outmat, plan else - error("Cannot handle noncomplex input matrix.") + error("Input must be complex valued.") end end -------------------------------------- -function num.fft3(datavec, n1, n2, n3) - if ffi.istype(gsl_matrix, datavec) then - return num.rfft3(datavec, n1, n2, n3) - else - local outvec = matrix.cnew(n1*n2*n3, 1) +function num.fftn(datavec, dimvec) + if ffi.istype(gsl_matrix_complex, datavec) then + local outvec = matrix.calloc(tonumber(datavec.size1), 1) local output = ffi.cast("fftw_complex*",outvec.data) local input = ffi.cast("fftw_complex*", datavec.data) - local plan = fftw.plan_dft_3d(n1,n2,n3, input, output, fftw.FORWARD, fftw.MEASURE) - fftw.execute(plan) - fftw.destroy_plan(plan) - return outvec - end -end - -function num.fft3inv(datavec, n1, n2, n3) - if ffi.istype(gsl_matrix, datavec) then - return num.rfft3_inv(datavec, n1, n2, n3) - else - local outvec = matrix.cnew(n1*n2*n3, 1) - local output = ffi.cast("fftw_complex*",outvec.data) - local input = ffi.cast("fftw_complex*", datavec.data) - local plan = fftw.plan_dft_3d(n1,n2,n3, input, output, fftw.BACKWARD, fftw.MEASURE) - fftw.execute(plan) - fftw.destroy_plan(plan) - return outvec/(n1*n2*n3) - end -end - -function num.rfft3(datavec, n1, n2, n3) - if ffi.istype(gsl_matrix, datavec) then - local outvec = matrix.cnew(n1*n2*(math.floor(n3/2)+1), 1) - local output = ffi.cast("fftw_complex*",outvec.data) - local plan = fftw.plan_dft_r2c_3d(n1,n2,n3, datavec.data, output, fftw.MEASURE) + local plan = fftw.plan_dft(#dimvec, dimvec.data, input, output, fftw.FORWARD, fftw.MEASURE) fftw.execute(plan) - fftw.destroy_plan(plan) - return outvec + return outvec, plan else - error("Cannot handle nonreal input data.") + error("Input must be complex valued.") end end -function num.rfft3inv(datavec, n1, n2, n3) - if ffi.istype(gsl_matrix, datavec) then - local outvec = matrix.new(n1*n2*n3, 1) - local input = ffi.cast("fftw_complex*",datavec.data) - local plan = fftw.plan_dft_r2c_3d(n1,n2,n3, input, outvec.data, fftw.MEASURE) - fftw.execute(plan) - fftw.destroy_plan(plan) - return outvec/(n1*n2*n3) - else - error("Cannot handle nonreal input data.") - end -end - --------------------------------------- - -function num.fftn(datavec, dimensionvec) - if ffi.istype(gsl_matrix, datavec) then - return num.rfftn(datavec, dimensionvec) - else - local outvec = matrix.cnew(tonumber(datavec.size1), 1) +function num.fftninv(datavec, dimvec) + if ffi.istype(gsl_matrix_complex, datavec) then + local outvec = matrix.calloc(tonumber(datavec.size1), 1) local output = ffi.cast("fftw_complex*",outvec.data) local input = ffi.cast("fftw_complex*", datavec.data) - local plan = fftw.plan_dft(#dimensionvec, dimensionvec.data, input, output, fftw.FORWARD, fftw.MEASURE) + local plan = fftw.plan_dft(#dimvec, dimvec.data, input, output, fftw.BACKWARD, fftw.MEASURE) fftw.execute(plan) - fftw.destroy_plan(plan) - return outvec - end -end - -function num.fftninv(datavec, dimensionvec) - if ffi.istype(gsl_matrix, datavec) then - return num.rfftn_inv(datavec, dimensionvec) + return outvec, plan else - local outvec = matrix.cnew(tonumber(datavec.size1), 1) - local output = ffi.cast("fftw_complex*",outvec.data) - local input = ffi.cast("fftw_complex*", datavec.data) - local plan = fftw.plan_dft(#dimensionvec, dimensionvec.data, input, output, fftw.BACKWARD, fftw.MEASURE) - fftw.execute(plan) - fftw.destroy_plan(plan) - return outvec + error("Input must be complex valued.") end end -function num.rfftn(datavec, dimensionvec) +function num.rfftn(datavec, dimvec) if ffi.istype(gsl_matrix, datavec) then - local d = #dimensionvec + local d = #dimvec local newsize = 1 for i=1,d-1 do - newsize = newsize * dimensionvec[i] + newsize = newsize * dimvec[i] end - newsize = newsize * math.floor(dimensionvec[d]/2)+1 + newsize = newsize * math.floor(dimvec[d]/2)+1 - local outvec = matrix.cnew(newsize, 1) + local outvec = matrix.calloc(newsize, 1) local output = ffi.cast("fftw_complex*",outvec.data) - local plan = fftw.plan_dft_r2c(d, dimensionvec.data, datavec.data, output, fftw.MEASURE) + local plan = fftw.plan_dft_r2c(d, dimvec.data, datavec.data, output, fftw.MEASURE) fftw.execute(plan) - fftw.destroy_plan(plan) - return outvec + return outvec, plan else - error("Cannot handle nonreal input data.") + error("Input must be real valued.") end end -function num.rfftninv(datavec, dimensionvec) - if ffi.istype(gsl_matrix, datavec) then - local d = #dimensionvec +function num.rfftninv(datavec, dimvec) + if ffi.istype(gsl_matrix_complex, datavec) then + local d = #dimvec local newsize = 1 for i=1,d do - newsize = newsize * dimensionvec[i] + newsize = newsize * dimvec[i] end - local outvec = matrix.new(newsize, 1) + local outvec = matrix.alloc(newsize, 1) local output = ffi.cast("fftw_complex*",outvec.data) - local plan = fftw.plan_dft_r2c(d, dimensionvec.data, datavec.data, output, fftw.MEASURE) + local plan = fftw.plan_dft_r2c(d, dimvec.data, datavec.data, output, fftw.MEASURE) fftw.execute(plan) - fftw.destroy_plan(plan) return outvec else - error("Cannot handle nonreal input data.") + error("Input must be complex valued.") end +end + +------------------------------------------ + +function num.fft_plan(plan, input, output) + local fftw_in = ffi.cast("fftw_complex*",input.data) + local fftw_out = ffi.cast("fftw_complex*",output.data) + return fftw.execute_dft(plan, fftw_in, fftw_out) +end + +function num.rfft_plan(plan, input, output) + local fftw_out = ffi.cast("fftw_complex*",output.data) + return fftw.execute_dft_r2c(plan, input.data, fftw_out) +end + +function num.rfftinv_plan(plan, input, output) + local fftw_in = ffi.cast("fftw_complex*",input.data) + return fftw.execute_dft_c2r(plan, fftw_in, output.data) end \ No newline at end of file From 12f28e208ff4e1f22480cdda8cc8450684c2bad9 Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Thu, 6 Feb 2014 22:42:02 +0100 Subject: [PATCH 07/20] Changed demo for usage of new FFT functions --- demos/fft.lua | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/demos/fft.lua b/demos/fft.lua index 3c0f189e0..f4c9ab55b 100644 --- a/demos/fft.lua +++ b/demos/fft.lua @@ -30,13 +30,13 @@ local function demo1() pt:addline(filine(|i| sq[i], n), 'black') - local ft = fft(sq) + local ft = rfft(sq) - local pf = fibars(|k| complex.abs(ft[k]), 0, n/2, 'black') + local pf = fibars(|k| complex.abs(ft[k]), 1, n/2+1, 'black') pf.title = 'FFT Power Spectrum' - for k=ncut, n - ncut do ft[k] = 0 end - sqt = fftinv(ft) + for k=ncut, n/2+1 do ft[k] = 0 end + local sqt = rfftinv(ft)/n pt:addline(filine(|i| sqt[i], n), 'red') pt:show() @@ -54,14 +54,14 @@ local function demo2() pt:addline(filine(|i| sq[i], n), 'black') - ft = fft(sq, true) + ft = rfft(sq) - pf:add(ibars(iter.isample(|k| complex.abs(ft[k]), 0, n/2)), 'black') + pf:add(ibars(iter.isample(|k| complex.abs(ft[k]), 1, #ft)), 'black') - for k=ncut, n - ncut do ft[k] = 0 end - fftinv(ft, true) + for k=ncut, #ft do ft[k] = 0 end + local sqt = rfftinv(ft)/n - pt:addline(filine(|i| sq[i], n), 'red') + pt:addline(filine(|i| sqt[i], n), 'red') pf:show() pt:show() @@ -79,15 +79,15 @@ local function demo3() local p = plot('Original signal / reconstructed') p:addline(filine(|i| bess[i], n), 'black') - local ft = fft(bess) + local ft = rfft(bess) fftplot = plot('FFT power spectrum') - bars = ibars(iter.isample(|k| complex.abs(ft[k]), 0, 60)) + bars = ibars(iter.isample(|k| complex.abs(ft[k]), 1, 60)) fftplot:add(bars, 'black') fftplot:show() - for k=ncut, n/2 do ft[k] = 0 end - local bessr = fftinv(ft) + for k=ncut, #ft do ft[k] = 0 end + local bessr = rfftinv(ft)/n p:addline(filine(|i| bessr[i], n), 'red', {{'dash', 7, 3}}) p:show() From 969f52545fbca3512b91ebba417801c46068911d Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Thu, 6 Feb 2014 23:33:07 +0100 Subject: [PATCH 08/20] Added garbage collection for all FFTW plans --- fft.lua | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/fft.lua b/fft.lua index 586da958c..cc80b3f1d 100644 --- a/fft.lua +++ b/fft.lua @@ -38,7 +38,7 @@ function num.rfft(vec) local n = tonumber(vec.size1) local outvec = matrix.calloc(math.floor(n/2)+1, 1) local output = ffi.cast("fftw_complex*",outvec.data) - local plan = fftw.plan_dft_r2c_1d(n, vec.data, output, fftw.MEASURE ) + local plan = ffi.gc(fftw.plan_dft_r2c_1d(n, vec.data, output, fftw.MEASURE), fftw.destroy_plan) fftw.execute(plan) return outvec, plan else @@ -52,7 +52,7 @@ function num.rfftinv(vec) local nnew = (n-1)*2 local outvec = matrix.alloc(nnew, 1) local input = ffi.cast("fftw_complex*",vec.data) - local plan = fftw.plan_dft_c2r_1d(nnew, input, outvec.data, fftw.MEASURE ) + local plan = ffi.gc(fftw.plan_dft_c2r_1d(nnew, input, outvec.data, fftw.MEASURE), fftw.destroy_plan) fftw.execute(plan) return outvec,plan else @@ -69,7 +69,7 @@ function num.fft2(mat) local outmat = matrix.calloc(n1, n2) local output = ffi.cast("fftw_complex*",outmat.data) local input = ffi.cast("fftw_complex*", mat.data) - local plan = fftw.plan_dft_2d(n1,n2, input, output, fftw.FORWARD, fftw.MEASURE) + local plan = ffi.gc(fftw.plan_dft_2d(n1,n2, input, output, fftw.FORWARD, fftw.MEASURE), fftw.destroy_plan) fftw.execute(plan) return outmat, plan else @@ -84,7 +84,7 @@ function num.fft2inv(mat) local outmat = matrix.calloc(n1, n2) local output = ffi.cast("fftw_complex*",outmat.data) local input = ffi.cast("fftw_complex*", mat.data) - local plan = fftw.plan_dft_2d(n1,n2, input, output, fftw.BACKWARD, fftw.MEASURE) + local plan = ffi.gc(fftw.plan_dft_2d(n1,n2, input, output, fftw.BACKWARD, fftw.MEASURE), fftw.destroy_plan) fftw.execute(plan) return outmat, plan else @@ -98,7 +98,7 @@ function num.rfft2(mat) local n2 = tonumber(mat.size2) local outmat = matrix.calloc(n1, math.floor(n2/2)+1) local output = ffi.cast("fftw_complex*",outmat.data) - local plan = fftw.plan_dft_r2c_2d(n1,n2, mat.data, output,fftw.MEASURE) + local plan = ffi.gc(fftw.plan_dft_r2c_2d(n1,n2, mat.data, output,fftw.MEASURE), fftw.destroy_plan) fftw.execute(plan) return outmat, plan else @@ -113,7 +113,7 @@ function num.rfft2inv(mat) local n2new = (n2-1)*2 local outmat = matrix.alloc(n1, n2new) local input = ffi.cast("fftw_complex*",mat.data) - local plan = fftw.plan_dft_c2r_2d(n1,n2new, input, outmat.data,fftw.MEASURE) + local plan = ffi.gc(fftw.plan_dft_c2r_2d(n1,n2new, input, outmat.data,fftw.MEASURE), fftw.destroy_plan) fftw.execute(plan) return outmat, plan else @@ -128,7 +128,7 @@ function num.fftn(datavec, dimvec) local outvec = matrix.calloc(tonumber(datavec.size1), 1) local output = ffi.cast("fftw_complex*",outvec.data) local input = ffi.cast("fftw_complex*", datavec.data) - local plan = fftw.plan_dft(#dimvec, dimvec.data, input, output, fftw.FORWARD, fftw.MEASURE) + local plan = ffi.gc(fftw.plan_dft(#dimvec, dimvec.data, input, output, fftw.FORWARD, fftw.MEASURE), fftw.destroy_plan) fftw.execute(plan) return outvec, plan else @@ -141,7 +141,7 @@ function num.fftninv(datavec, dimvec) local outvec = matrix.calloc(tonumber(datavec.size1), 1) local output = ffi.cast("fftw_complex*",outvec.data) local input = ffi.cast("fftw_complex*", datavec.data) - local plan = fftw.plan_dft(#dimvec, dimvec.data, input, output, fftw.BACKWARD, fftw.MEASURE) + local plan = ffi.gc(fftw.plan_dft(#dimvec, dimvec.data, input, output, fftw.BACKWARD, fftw.MEASURE), fftw.destroy_plan) fftw.execute(plan) return outvec, plan else @@ -160,7 +160,7 @@ function num.rfftn(datavec, dimvec) local outvec = matrix.calloc(newsize, 1) local output = ffi.cast("fftw_complex*",outvec.data) - local plan = fftw.plan_dft_r2c(d, dimvec.data, datavec.data, output, fftw.MEASURE) + local plan = ffi.gc(fftw.plan_dft_r2c(d, dimvec.data, datavec.data, output, fftw.MEASURE), fftw.destroy_plan) fftw.execute(plan) return outvec, plan else @@ -178,7 +178,7 @@ function num.rfftninv(datavec, dimvec) local outvec = matrix.alloc(newsize, 1) local output = ffi.cast("fftw_complex*",outvec.data) - local plan = fftw.plan_dft_r2c(d, dimvec.data, datavec.data, output, fftw.MEASURE) + local plan = ffi.gc(fftw.plan_dft_r2c(d, dimvec.data, datavec.data, output, fftw.MEASURE), fftw.destroy_plan) fftw.execute(plan) return outvec else From 4f1f3285f8ad81e02b6c5a4e8d00af3d682b8c86 Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Fri, 7 Feb 2014 00:03:27 +0100 Subject: [PATCH 09/20] Fixed dimension list for n-dimensional ffts --- fft.lua | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/fft.lua b/fft.lua index cc80b3f1d..e7a487653 100644 --- a/fft.lua +++ b/fft.lua @@ -123,12 +123,12 @@ end -------------------------------------- -function num.fftn(datavec, dimvec) +function num.fftn(datavec, dimlist) if ffi.istype(gsl_matrix_complex, datavec) then local outvec = matrix.calloc(tonumber(datavec.size1), 1) local output = ffi.cast("fftw_complex*",outvec.data) local input = ffi.cast("fftw_complex*", datavec.data) - local plan = ffi.gc(fftw.plan_dft(#dimvec, dimvec.data, input, output, fftw.FORWARD, fftw.MEASURE), fftw.destroy_plan) + local plan = ffi.gc(fftw.plan_dft(#dimlist, ffi.new("int[?]", #dimlist, dimlist), input, output, fftw.FORWARD, fftw.MEASURE), fftw.destroy_plan) fftw.execute(plan) return outvec, plan else @@ -136,12 +136,13 @@ function num.fftn(datavec, dimvec) end end -function num.fftninv(datavec, dimvec) +function num.fftninv(datavec, dimlist) if ffi.istype(gsl_matrix_complex, datavec) then local outvec = matrix.calloc(tonumber(datavec.size1), 1) local output = ffi.cast("fftw_complex*",outvec.data) local input = ffi.cast("fftw_complex*", datavec.data) - local plan = ffi.gc(fftw.plan_dft(#dimvec, dimvec.data, input, output, fftw.BACKWARD, fftw.MEASURE), fftw.destroy_plan) + + local plan = ffi.gc(fftw.plan_dft(#dimlist, ffi.new("int[?]", #dimlist, dimlist), input, output, fftw.BACKWARD, fftw.MEASURE), fftw.destroy_plan) fftw.execute(plan) return outvec, plan else @@ -149,18 +150,18 @@ function num.fftninv(datavec, dimvec) end end -function num.rfftn(datavec, dimvec) +function num.rfftn(datavec, dimlist) if ffi.istype(gsl_matrix, datavec) then - local d = #dimvec + local d = #dimlist local newsize = 1 for i=1,d-1 do - newsize = newsize * dimvec[i] + newsize = newsize * dimlist[i] end - newsize = newsize * math.floor(dimvec[d]/2)+1 + newsize = newsize * math.floor(dimlist[d]/2)+1 local outvec = matrix.calloc(newsize, 1) local output = ffi.cast("fftw_complex*",outvec.data) - local plan = ffi.gc(fftw.plan_dft_r2c(d, dimvec.data, datavec.data, output, fftw.MEASURE), fftw.destroy_plan) + local plan = ffi.gc(fftw.plan_dft_r2c(d, ffi.new("int[?]", #dimlist, dimlist), datavec.data, output, fftw.MEASURE), fftw.destroy_plan) fftw.execute(plan) return outvec, plan else @@ -168,17 +169,17 @@ function num.rfftn(datavec, dimvec) end end -function num.rfftninv(datavec, dimvec) +function num.rfftninv(datavec, dimlist) if ffi.istype(gsl_matrix_complex, datavec) then - local d = #dimvec + local d = #dimlist local newsize = 1 for i=1,d do - newsize = newsize * dimvec[i] + newsize = newsize * dimlist[i] end local outvec = matrix.alloc(newsize, 1) local output = ffi.cast("fftw_complex*",outvec.data) - local plan = ffi.gc(fftw.plan_dft_r2c(d, dimvec.data, datavec.data, output, fftw.MEASURE), fftw.destroy_plan) + local plan = ffi.gc(fftw.plan_dft_r2c(d, ffi.new("int[?]", #dimlist, dimlist), datavec.data, output, fftw.MEASURE), fftw.destroy_plan) fftw.execute(plan) return outvec else From 18ce5b1a4108e5a3316fb3cd70e9f4a78f6b40b7 Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Thu, 13 Feb 2014 18:35:52 +0100 Subject: [PATCH 10/20] Rewritten Fourier functions. Now 3 fundamental functions exist: dft, rdft, rdftinv which are called by the fft, fft2 and fftn functions. Added test file and removed old fft-tests which was the same as the demo file. --- fft.lua | 339 +++++++++++++++++++++++++------------------- tests/fft-tests.lua | 187 ++++-------------------- 2 files changed, 220 insertions(+), 306 deletions(-) diff --git a/fft.lua b/fft.lua index e7a487653..ca32bcbc8 100644 --- a/fft.lua +++ b/fft.lua @@ -3,204 +3,253 @@ fftw = require'fftw3.init' local gsl_matrix = ffi.typeof('gsl_matrix') local gsl_matrix_complex = ffi.typeof('gsl_matrix_complex') -function num.fft(vec) - if ffi.istype(gsl_matrix_complex, vec) then - local n = tonumber(vec.size1) - local outvec = matrix.calloc(n, 1) - local output = ffi.cast("fftw_complex*",outvec.data) - local input = ffi.cast("fftw_complex*",vec.data) +local FORWARD = 1 +local BACKWARD = 2 - local plan = ffi.gc(fftw.plan_dft_1d(n, input, output, fftw.FORWARD, fftw.MEASURE), fftw.destroy_plan) - fftw.execute(plan) - return outvec, plan - else - error("Input must be complex valued.") + +local function compare_dimensions(dim1, dim2) + for i = 1,#dim1 do + if dim1[i] ~= dim2[i] then + return false + end end + return true end -function num.fftinv(vec) - if ffi.istype(gsl_matrix_complex, vec) then - local n = tonumber(vec.size1) - local outvec = matrix.calloc(n, 1) +local dftplans = {} +local function dft(invec, dimlist, outvec, direction) + + if ffi.istype(gsl_matrix_complex, invec) then + local rank = #dimlist + local size = 1 + for i=1,rank do size = size * dimlist[i] end + local direction = direction or FORWARD + + local outvec = outvec or nil + local output_supplied = true + + --First check if output vector is given by the user, otherwise allocate it + if outvec == nil then + output_supplied = false + outvec = matrix.calloc(size, 1) + end + local output = ffi.cast("fftw_complex*",outvec.data) - local input = ffi.cast("fftw_complex*",vec.data) + local input = ffi.cast("fftw_complex*",invec.data) + + --Check if plan is existing, then use it + if dftplans[rank] ~= nil and dftplans[rank][direction] ~= nil and compare_dimensions(dimlist, dftplans[rank][direction].dimlist) then + local plan = dftplans[rank][direction].plan + fftw.execute_dft(plan, input, output) + if not output_supplied then return outvec end + + --Or create a new plan + else + --allcoate input for plan making but planning can invalidate input, therefore we use extra array + local inputtest = ffi.new("fftw_complex[?]", size) - local plan = ffi.gc(fftw.plan_dft_1d(n, input, output, fftw.BACKWARD, fftw.MEASURE), fftw.destroy_plan) - fftw.execute(plan) - return outvec, plan + --Create plan via MEASURE + local plan = ffi.gc(fftw.plan_dft(rank, ffi.new("int[?]", rank, dimlist), inputtest, output, direction == FORWARD and fftw.FORWARD or fftw.BACKWARD, fftw.MEASURE), fftw.destroy_plan) + + --Save the plan for next time + if dftplans[rank] == nil then dftplans[rank] = {} end + dftplans[rank][direction] = {plan=plan, dimlist=dimlist} + + --Execute the plan with the supplied input + fftw.execute_dft(plan, input, output) + if not output_supplied then return outvec end + end else error("Input must be complex valued.") end end -function num.rfft(vec) - if ffi.istype(gsl_matrix, vec) then - local n = tonumber(vec.size1) - local outvec = matrix.calloc(math.floor(n/2)+1, 1) +local rdftplans = {} +local function rdft(invec, dimlist, outvec) + + if ffi.istype(gsl_matrix, invec) then + local rank = #dimlist + local size = 1 + for i=1,rank do size = size * dimlist[i] end + local outputsize = size/dimlist[rank]*(math.floor(dimlist[rank]/2)+1) + + local direction = FORWARD + + local outvec = outvec or nil + local output_supplied = true + + --First check if output vector is given by the user, otherwise allocate it + if outvec == nil then + output_supplied = false + outvec = matrix.calloc(outputsize, 1) + end + local output = ffi.cast("fftw_complex*",outvec.data) - local plan = ffi.gc(fftw.plan_dft_r2c_1d(n, vec.data, output, fftw.MEASURE), fftw.destroy_plan) - fftw.execute(plan) - return outvec, plan + local input = invec.data + + --Check if plan is existing, then use it + if rdftplans[rank] ~= nil and rdftplans[rank][direction] ~= nil and compare_dimensions(dimlist, rdftplans[rank][direction].dimlist) then + local plan = rdftplans[rank][direction].plan + fftw.execute_dft_r2c(plan, input, output) + if not output_supplied then return outvec end + + --Or create a new plan + else + --allcoate input for plan making but planning can invalidate input, therefore we use extra array + local inputtest = ffi.new("double[?]", size) + + --Create plan via MEASURE + local plan = ffi.gc(fftw.plan_dft_r2c(rank, ffi.new("int[?]", rank, dimlist), inputtest, output, fftw.MEASURE), fftw.destroy_plan) + + --Save the plan for next time + if rdftplans[rank] == nil then rdftplans[rank] = {} end + rdftplans[rank][direction] = {plan=plan, dimlist=dimlist} + + --Execute the plan with the supplied input + fftw.execute_dft_r2c(plan, input, output) + if not output_supplied then return outvec end + end else error("Input must be real valued.") end end -function num.rfftinv(vec) - if ffi.istype(gsl_matrix_complex, vec) then - local n = tonumber(vec.size1) - local nnew = (n-1)*2 - local outvec = matrix.alloc(nnew, 1) - local input = ffi.cast("fftw_complex*",vec.data) - local plan = ffi.gc(fftw.plan_dft_c2r_1d(nnew, input, outvec.data, fftw.MEASURE), fftw.destroy_plan) - fftw.execute(plan) - return outvec,plan +local function rdftinv(invec, dimlist, outvec) + + if ffi.istype(gsl_matrix_complex, invec) then + local rank = #dimlist + local size = 1 + for i=1,rank do size = size * dimlist[i] end + local inputsize = size/dimlist[rank]*(math.floor(dimlist[rank]/2)+1) + + local direction = BACKWARD + + local outvec = outvec or nil + local output_supplied = true + + --First check if output vector is given by the user, otherwise allocate it + if outvec == nil then + output_supplied = false + outvec = matrix.alloc(size, 1) + end + + local output = outvec.data + local input = ffi.cast("fftw_complex*",invec.data) + + --Check if plan is existing, then use it + if rdftplans[rank] ~= nil and rdftplans[rank][direction] ~= nil and compare_dimensions(dimlist, rdftplans[rank][direction].dimlist) then + local plan = rdftplans[rank][direction].plan + fftw.execute_dft_c2r(plan, input, output) + if not output_supplied then return outvec end + + --Or create a new plan + else + --allcoate input for plan making but planning can invalidate input, therefore we use extra array + local inputtest = ffi.new("fftw_complex[?]", inputsize) + + --Create plan via MEASURE + local plan = ffi.gc(fftw.plan_dft_c2r(rank, ffi.new("int[?]", rank, dimlist), inputtest, output, fftw.MEASURE), fftw.destroy_plan) + + --Save the plan for next time + if rdftplans[rank] == nil then rdftplans[rank] = {} end + rdftplans[rank][direction] = {plan=plan, dimlist=dimlist} + + --Execute the plan with the supplied input + fftw.execute_dft_c2r(plan, input, output) + if not output_supplied then return outvec end + end else error("Input must be complex valued.") end end ----------------------------------- +------------------------------------------------ -function num.fft2(mat) - if ffi.istype(gsl_matrix_complex, mat) then - local n1 = tonumber(mat.size1) - local n2 = tonumber(mat.size2) - local outmat = matrix.calloc(n1, n2) - local output = ffi.cast("fftw_complex*",outmat.data) - local input = ffi.cast("fftw_complex*", mat.data) - local plan = ffi.gc(fftw.plan_dft_2d(n1,n2, input, output, fftw.FORWARD, fftw.MEASURE), fftw.destroy_plan) - fftw.execute(plan) - return outmat, plan - else - error("Input must be complex valued.") - end +function num.fft(invec, outvec) + return dft(invec, {#invec}, outvec, FORWARD) end -function num.fft2inv(mat) - if ffi.istype(gsl_matrix_complex, mat) then - local n1 = tonumber(mat.size1) - local n2 = tonumber(mat.size2) - local outmat = matrix.calloc(n1, n2) - local output = ffi.cast("fftw_complex*",outmat.data) - local input = ffi.cast("fftw_complex*", mat.data) - local plan = ffi.gc(fftw.plan_dft_2d(n1,n2, input, output, fftw.BACKWARD, fftw.MEASURE), fftw.destroy_plan) - fftw.execute(plan) - return outmat, plan - else - error("Input must be complex valued.") - end +function num.fftinv(invec, outvec) + return dft(invec, {#invec}, outvec, BACKWARD) end -function num.rfft2(mat) - if ffi.istype(gsl_matrix, mat) then - local n1 = tonumber(mat.size1) - local n2 = tonumber(mat.size2) - local outmat = matrix.calloc(n1, math.floor(n2/2)+1) - local output = ffi.cast("fftw_complex*",outmat.data) - local plan = ffi.gc(fftw.plan_dft_r2c_2d(n1,n2, mat.data, output,fftw.MEASURE), fftw.destroy_plan) - fftw.execute(plan) - return outmat, plan - else - error("Input must be real valued.") - end +function num.rfft(invec, outvec) + return rdft(invec, {#invec}, outvec) end -function num.rfft2inv(mat) - if ffi.istype(gsl_matrix_complex, mat) then - local n1 = tonumber(mat.size1) - local n2 = tonumber(mat.size2) - local n2new = (n2-1)*2 - local outmat = matrix.alloc(n1, n2new) - local input = ffi.cast("fftw_complex*",mat.data) - local plan = ffi.gc(fftw.plan_dft_c2r_2d(n1,n2new, input, outmat.data,fftw.MEASURE), fftw.destroy_plan) - fftw.execute(plan) - return outmat, plan - else - error("Input must be complex valued.") - end +function num.rfftinv(invec, outvec) + return rdftinv(invec, {(#invec-1)*2}, outvec) end --------------------------------------- +function num.fft2(inmat, outmat) + local n1 = tonumber(inmat.size1) + local n2 = tonumber(inmat.size2) -function num.fftn(datavec, dimlist) - if ffi.istype(gsl_matrix_complex, datavec) then - local outvec = matrix.calloc(tonumber(datavec.size1), 1) - local output = ffi.cast("fftw_complex*",outvec.data) - local input = ffi.cast("fftw_complex*", datavec.data) - local plan = ffi.gc(fftw.plan_dft(#dimlist, ffi.new("int[?]", #dimlist, dimlist), input, output, fftw.FORWARD, fftw.MEASURE), fftw.destroy_plan) - fftw.execute(plan) - return outvec, plan + if outmat == nil then + local retvec = dft(inmat, {n1, n2}, outmat, FORWARD) + local b = ffi.cast('gsl_block_complex*', ffi.C.malloc(ffi.sizeof('gsl_block_complex'))) + b.size, b.data, b.ref_count = n1*n2, retvec.data, 2 + return gsl_matrix_complex(n1, n2, n2, b.data, b, 1) else - error("Input must be complex valued.") + dft(inmat, {n1, n2}, outmat, FORWARD) end end -function num.fftninv(datavec, dimlist) - if ffi.istype(gsl_matrix_complex, datavec) then - local outvec = matrix.calloc(tonumber(datavec.size1), 1) - local output = ffi.cast("fftw_complex*",outvec.data) - local input = ffi.cast("fftw_complex*", datavec.data) +function num.fft2inv(inmat, outmat) + local n1 = tonumber(inmat.size1) + local n2 = tonumber(inmat.size2) - local plan = ffi.gc(fftw.plan_dft(#dimlist, ffi.new("int[?]", #dimlist, dimlist), input, output, fftw.BACKWARD, fftw.MEASURE), fftw.destroy_plan) - fftw.execute(plan) - return outvec, plan + if outmat == nil then + local retvec = dft(inmat, {n1, n2}, outmat, BACKWARD) + local b = ffi.cast('gsl_block_complex*', ffi.C.malloc(ffi.sizeof('gsl_block_complex'))) + b.size, b.data, b.ref_count = n1*n2, retvec.data, 2 + return gsl_matrix_complex(n1, n2, n2, b.data, b, 1) else - error("Input must be complex valued.") + dft(inmat, {n1, n2}, outmat, BACKWARD) end end -function num.rfftn(datavec, dimlist) - if ffi.istype(gsl_matrix, datavec) then - local d = #dimlist - local newsize = 1 - for i=1,d-1 do - newsize = newsize * dimlist[i] - end - newsize = newsize * math.floor(dimlist[d]/2)+1 +function num.rfft2(inmat, outmat) + local n1 = tonumber(inmat.size1) + local n2 = tonumber(inmat.size2) - local outvec = matrix.calloc(newsize, 1) - local output = ffi.cast("fftw_complex*",outvec.data) - local plan = ffi.gc(fftw.plan_dft_r2c(d, ffi.new("int[?]", #dimlist, dimlist), datavec.data, output, fftw.MEASURE), fftw.destroy_plan) - fftw.execute(plan) - return outvec, plan + if outmat == nil then + local retvec = rdft(inmat, {n1, n2}, outmat) + local b = ffi.cast('gsl_block_complex*', ffi.C.malloc(ffi.sizeof('gsl_block_complex'))) + local n2new = math.floor(n2/2)+1 + b.size, b.data, b.ref_count = n1*n2new, retvec.data, 2 + return gsl_matrix_complex(n1, n2new, n2new, b.data, b, 1) else - error("Input must be real valued.") + rdft(inmat, {n1, n2}, outmat) end end -function num.rfftninv(datavec, dimlist) - if ffi.istype(gsl_matrix_complex, datavec) then - local d = #dimlist - local newsize = 1 - for i=1,d do - newsize = newsize * dimlist[i] - end +function num.rfft2inv(inmat, outmat) + local n1 = tonumber(inmat.size1) + local n2 = (tonumber(inmat.size2)-1)*2 - local outvec = matrix.alloc(newsize, 1) - local output = ffi.cast("fftw_complex*",outvec.data) - local plan = ffi.gc(fftw.plan_dft_r2c(d, ffi.new("int[?]", #dimlist, dimlist), datavec.data, output, fftw.MEASURE), fftw.destroy_plan) - fftw.execute(plan) - return outvec + if outmat == nil then + local retvec = rdftinv(inmat, {n1, n2}, outmat) + local b = ffi.cast('gsl_block*', ffi.C.malloc(ffi.sizeof('gsl_block'))) + b.size, b.data, b.ref_count = n1*n2, retvec.data, 2 + return gsl_matrix(n1, n2, n2, b.data, b, 1) else - error("Input must be complex valued.") + rdft(inmat, {n1, n2}, outmat) end end ------------------------------------------- +function num.fftn(invec, dimlist, outvec) + return dft(invec, dimlist, outvec, FORWARD) +end -function num.fft_plan(plan, input, output) - local fftw_in = ffi.cast("fftw_complex*",input.data) - local fftw_out = ffi.cast("fftw_complex*",output.data) - return fftw.execute_dft(plan, fftw_in, fftw_out) +function num.fftninv(invec, dimlist, outvec) + return dft(invec, dimlist, outvec, BACKWARD) end -function num.rfft_plan(plan, input, output) - local fftw_out = ffi.cast("fftw_complex*",output.data) - return fftw.execute_dft_r2c(plan, input.data, fftw_out) +function num.rfftn(invec, dimlist, outvec) + return rdft(invec, dimlist, outvec) end -function num.rfftinv_plan(plan, input, output) - local fftw_in = ffi.cast("fftw_complex*",input.data) - return fftw.execute_dft_c2r(plan, fftw_in, output.data) -end \ No newline at end of file +function num.rfftninv(invec, dimlist, outvec) + return rdftinv(invec, dimlist, outvec) +end diff --git a/tests/fft-tests.lua b/tests/fft-tests.lua index b777797e3..161d6aeb7 100644 --- a/tests/fft-tests.lua +++ b/tests/fft-tests.lua @@ -1,172 +1,37 @@ -use 'math' -use 'graph' -use 'gsl' +--This is a test for all the fft functions in the library -local function test1() - local n, ncut = 8*3*5, 16 +local size = 10 +local size2 = 5 - local sq = matrix.new(n, 1, |i| i < n/3 and 0 or (i < 2*n/3 and 1 or 0)) +local r1d = matrix.new(size, 1, |i,j|i) +local c1d = matrix.cnew(size,1, |i,j|i) - local pt = plot('Original signal / reconstructed') - local pf = plot('FFT Power Spectrum') +local r2d = matrix.new(size,size, |i,j|i) +local c2d = matrix.cnew(size,size, |i,j|i) - pt:addline(filine(|i| sq[i], n), 'black') +local r2dasym = matrix.new(size,5, |i,j|i) +local c2dasym = matrix.new(size,5, |i,j|i) - local ft = fft(sq) +local r3d = matrix.new(size2*size2*size2, 1, |i,j|i) +local c3d = matrix.cnew(size2*size2*size2,1, |i,j|i) - pf:add(ibars(isample(|k| complex.abs(ft:get(k)), 0, n/2)), 'black') +-- 1D test - for k=ncut, n - ncut do ft:set(k,0) end - sqt = fftinv(ft) +print(num.fftinv(num.fft(c1d))/size) +print('---') +print(num.rfftinv(num.rfft(r1d))/size) +print('---') - pt:addline(filine(|i| sqt[i], n), 'red') +-- 2D tests - pf:show() - pt:show() +print(num.fft2inv(num.fft2(c2d))/(size*size)) +print('---') +print(num.rfft2inv(num.rfft2(r2d))/(size*size)) +print('---') - return pt, pf -end +-- 3D tests -local function test1_radix2() - local n, ncut = 256, 16 - - local sq = matrix.new(n, 1, |i| i < n/3 and 0 or (i < 2*n/3 and 1 or 0)) - - local pt = plot('Original signal / reconstructed') - local pf = plot('FFT Power Spectrum') - - pt:addline(filine(|i| sq[i], n), 'black') - - local ft = fft(sq) - - pf:add(ibars(isample(|k| complex.abs(ft:get(k)), 0, n/2)), 'black') - - for k=ncut, n - ncut do ft:set(k,0) end - sqt = fftinv(ft) - - pt:addline(filine(|i| sqt[i], n), 'red') - - pf:show() - pt:show() - - return pt, pf -end - -local function test1_ip_radix2() - local hcget, hcset = gsl.halfcomplex_radix2_get, gsl.halfcomplex_radix2_set - - local n, ncut = 256, 16 - - local sq = matrix.new(n, 1, |i| i < n/3 and 0 or (i < 2*n/3 and 1 or 0)) - - local pt = plot('Original signal / reconstructed') - local pf = plot('FFT Power Spectrum') - - pt:addline(filine(|i| sq[i], n), 'black') - - fft_radix2(sq) - - pf:add(ibars(isample(|k| complex.abs(hcget(sq, k)), 0, n/2)), 'black') - - for k=ncut, n - ncut do hcset(sq, k, 0) end - fft_radix2_inverse(sq) - - pt:addline(filine(|i| sq[i], n), 'red') - - pf:show() - pt:show() - - return pt, pf -end - -local function test1_ip() - local hcget, hcset = gsl.halfcomplex_get, gsl.halfcomplex_set - - local n, ncut = 8*3*5, 16 - - local sq = matrix.new(n, 1, |i| i < n/3 and 0 or (i < 2*n/3 and 1 or 0)) - - local pt = plot('Original signal / reconstructed') - local pf = plot('FFT Power Spectrum') - - pt:addline(filine(|i| sq[i], n), 'black') - - fft_real(sq) - - pf:add(ibars(isample(|k| complex.abs(hcget(sq, k)), 0, n/2)), 'black') - - for k=ncut, n - ncut do hcset(sq, k, 0) end - fft_halfcomplex_inverse(sq) - - pt:addline(filine(|i| sq[i], n), 'red') - - pf:show() - pt:show() - - return pt, pf -end - -local function test2() - local n, ncut, order = 512, 11, 8 - local x1 = besselJ_zero(order, 14) - local xsmp = |k| x1*(k-1)/(n-1) - - local bess = matrix.new(n, 1, |i| besselJ(order, xsmp(i))) - - local p = plot('Original signal / reconstructed') - p:addline(filine(|i| bess[i], n), 'black') - - local ft = fft(bess) - - fftplot = plot('FFT power spectrum') - bars = ibars(isample(|k| complex.abs(ft:get(k)), 0, 60)) - fftplot:add(bars, 'darkgreen') - fftplot:addline(bars, 'black') - fftplot:show() - - for k=ncut, n/2 do ft:set(k,0) end - local bessr = fftinv(ft) - - p:addline(filine(|i| bessr[i], n), 'red', {{'dash', 7, 3}}) - p:show() - - return p, fftplot -end - - -local function test1_stride() - local n, ncut, nb = 256, 16, 3 - - local function squaref(i, n1, n2) - return i < n1 and 0 or (i < n2 and 1 or 0) - end - - local sq = matrix.new(n, nb, |i, j| squaref(i, n/(3*j), 2*n/(3*j))) - - local w = window('v' .. string.rep('.', nb)) - local ftt = {} - for j=1, nb do - local ft = fft(sq:col(j)) - local pf = plot() - pf:add(ibars(isample(|k| complex.abs(ft:get(k)), 0, n/2)), 'black') - w:attach(pf, j) - ftt[j] = ft - end - - local w = window('v' .. string.rep('.', nb)) - for j, ft in ipairs(ftt) do - local pt = plot('Original signal / reconstructed') - for k=ncut, n - ncut do ft:set(k,0) end - local sqt = fftinv(ft) - pt:addline(filine(|i| sq:get(i,j), n), 'black') - pt:addline(filine(|i| sqt[i], n)) - w:attach(pt, j) - end -end - -return {test1= test1, - test2= test1_radix2, - test3= test1_ip_radix2, - test4= test1_ip, - test5= test2, - test6= test1_stride} +print('---') +print(num.fftninv(num.fftn(c3d, {size2, size2, size2}), {size2, size2, size2})/(size2*size2*size2)) +print('---') +print(num.rfftninv(num.rfftn(r3d, {size2, size2, size2}), {size2, size2, size2})/(size2*size2*size2)) \ No newline at end of file From cdcc53be5c03d89299e2b3154595dbaec65cd7d8 Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Thu, 13 Feb 2014 22:45:29 +0100 Subject: [PATCH 11/20] Added UNALLIGNED flat to the FFTW calls to prevent segment fault errors. --- fft.lua | 25 ++++++++++--------------- 1 file changed, 10 insertions(+), 15 deletions(-) diff --git a/fft.lua b/fft.lua index ca32bcbc8..348eb9eeb 100644 --- a/fft.lua +++ b/fft.lua @@ -1,5 +1,5 @@ -ffi = require'ffi' -fftw = require'fftw3.init' +local ffi = require'ffi' +local fftw = require'fftw3.init' local gsl_matrix = ffi.typeof('gsl_matrix') local gsl_matrix_complex = ffi.typeof('gsl_matrix_complex') @@ -24,12 +24,10 @@ local function dft(invec, dimlist, outvec, direction) local size = 1 for i=1,rank do size = size * dimlist[i] end local direction = direction or FORWARD - - local outvec = outvec or nil local output_supplied = true --First check if output vector is given by the user, otherwise allocate it - if outvec == nil then + if not outvec then output_supplied = false outvec = matrix.calloc(size, 1) end @@ -49,7 +47,7 @@ local function dft(invec, dimlist, outvec, direction) local inputtest = ffi.new("fftw_complex[?]", size) --Create plan via MEASURE - local plan = ffi.gc(fftw.plan_dft(rank, ffi.new("int[?]", rank, dimlist), inputtest, output, direction == FORWARD and fftw.FORWARD or fftw.BACKWARD, fftw.MEASURE), fftw.destroy_plan) + local plan = ffi.gc(fftw.plan_dft(rank, ffi.new("int[?]", rank, dimlist), inputtest, output, direction == FORWARD and fftw.FORWARD or fftw.BACKWARD, bit.bor(fftw.MEASURE, fftw.UNALIGNED)), fftw.destroy_plan) --Save the plan for next time if dftplans[rank] == nil then dftplans[rank] = {} end @@ -74,12 +72,10 @@ local function rdft(invec, dimlist, outvec) local outputsize = size/dimlist[rank]*(math.floor(dimlist[rank]/2)+1) local direction = FORWARD - - local outvec = outvec or nil local output_supplied = true --First check if output vector is given by the user, otherwise allocate it - if outvec == nil then + if not outvec then output_supplied = false outvec = matrix.calloc(outputsize, 1) end @@ -96,10 +92,11 @@ local function rdft(invec, dimlist, outvec) --Or create a new plan else --allcoate input for plan making but planning can invalidate input, therefore we use extra array - local inputtest = ffi.new("double[?]", size) + local inputtest = ffi.new("double[?]", size) + -- local inputtest = ffi.gc(fftw.malloc(size * ffi.sizeof(double)), fftw.free) --Create plan via MEASURE - local plan = ffi.gc(fftw.plan_dft_r2c(rank, ffi.new("int[?]", rank, dimlist), inputtest, output, fftw.MEASURE), fftw.destroy_plan) + local plan = ffi.gc(fftw.plan_dft_r2c(rank, ffi.new("int[?]", rank, dimlist), inputtest, output, bit.bor(fftw.MEASURE, fftw.UNALIGNED)), fftw.destroy_plan) --Save the plan for next time if rdftplans[rank] == nil then rdftplans[rank] = {} end @@ -123,12 +120,10 @@ local function rdftinv(invec, dimlist, outvec) local inputsize = size/dimlist[rank]*(math.floor(dimlist[rank]/2)+1) local direction = BACKWARD - - local outvec = outvec or nil local output_supplied = true --First check if output vector is given by the user, otherwise allocate it - if outvec == nil then + if not outvec then output_supplied = false outvec = matrix.alloc(size, 1) end @@ -148,7 +143,7 @@ local function rdftinv(invec, dimlist, outvec) local inputtest = ffi.new("fftw_complex[?]", inputsize) --Create plan via MEASURE - local plan = ffi.gc(fftw.plan_dft_c2r(rank, ffi.new("int[?]", rank, dimlist), inputtest, output, fftw.MEASURE), fftw.destroy_plan) + local plan = ffi.gc(fftw.plan_dft_c2r(rank, ffi.new("int[?]", rank, dimlist), inputtest, output, bit.bor(fftw.MEASURE, fftw.UNALIGNED)), fftw.destroy_plan) --Save the plan for next time if rdftplans[rank] == nil then rdftplans[rank] = {} end From e02a6c3e557c7fbf857c98313186919208df9608 Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Fri, 21 Feb 2014 11:07:23 +0100 Subject: [PATCH 12/20] Use advanced FFTW interface to create plans for using stride of input data (tda of gsl blocks)/ --- fft.lua | 64 +++++++++++++++++++++++++++++++++++++++++++++------------ 1 file changed, 51 insertions(+), 13 deletions(-) diff --git a/fft.lua b/fft.lua index 348eb9eeb..93b69d222 100644 --- a/fft.lua +++ b/fft.lua @@ -36,22 +36,34 @@ local function dft(invec, dimlist, outvec, direction) local input = ffi.cast("fftw_complex*",invec.data) --Check if plan is existing, then use it - if dftplans[rank] ~= nil and dftplans[rank][direction] ~= nil and compare_dimensions(dimlist, dftplans[rank][direction].dimlist) then + if dftplans[rank] ~= nil and dftplans[rank][direction] ~= nil and + compare_dimensions(dimlist, dftplans[rank][direction].dimlist) and + invec.tda == dftplans[rank][direction].itda and + outvec.tda == dftplans[rank][direction].otda then + local plan = dftplans[rank][direction].plan fftw.execute_dft(plan, input, output) if not output_supplied then return outvec end - --Or create a new plan else --allcoate input for plan making but planning can invalidate input, therefore we use extra array local inputtest = ffi.new("fftw_complex[?]", size) - --Create plan via MEASURE - local plan = ffi.gc(fftw.plan_dft(rank, ffi.new("int[?]", rank, dimlist), inputtest, output, direction == FORWARD and fftw.FORWARD or fftw.BACKWARD, bit.bor(fftw.MEASURE, fftw.UNALIGNED)), fftw.destroy_plan) + local howmany = 1 + local idist,odist = 0,0 + local istride = invec.tda + local ostride = outvec.tda + local dim = ffi.new("int[?]", rank, dimlist) + local inembed = dim + local onembed = dim + local plan = ffi.gc(fftw.plan_many_dft(rank, dim ,howmany, + inputtest,inembed, istride, idist, + output,onembed,ostride,odist, + direction == FORWARD and fftw.FORWARD or fftw.BACKWARD, bit.bor(fftw.MEASURE, fftw.UNALIGNED)), fftw.destroy_plan) --Save the plan for next time if dftplans[rank] == nil then dftplans[rank] = {} end - dftplans[rank][direction] = {plan=plan, dimlist=dimlist} + dftplans[rank][direction] = {plan=plan, dimlist=dimlist, itda=invec.tda, otda=outvec.tda} --Execute the plan with the supplied input fftw.execute_dft(plan, input, output) @@ -84,7 +96,11 @@ local function rdft(invec, dimlist, outvec) local input = invec.data --Check if plan is existing, then use it - if rdftplans[rank] ~= nil and rdftplans[rank][direction] ~= nil and compare_dimensions(dimlist, rdftplans[rank][direction].dimlist) then + if rdftplans[rank] ~= nil and rdftplans[rank][direction] ~= nil and + compare_dimensions(dimlist, rdftplans[rank][direction].dimlist) and + invec.tda == rdftplans[rank][direction].itda and + outvec.tda == rdftplans[rank][direction].otda then + local plan = rdftplans[rank][direction].plan fftw.execute_dft_r2c(plan, input, output) if not output_supplied then return outvec end @@ -95,12 +111,21 @@ local function rdft(invec, dimlist, outvec) local inputtest = ffi.new("double[?]", size) -- local inputtest = ffi.gc(fftw.malloc(size * ffi.sizeof(double)), fftw.free) - --Create plan via MEASURE - local plan = ffi.gc(fftw.plan_dft_r2c(rank, ffi.new("int[?]", rank, dimlist), inputtest, output, bit.bor(fftw.MEASURE, fftw.UNALIGNED)), fftw.destroy_plan) + local howmany = 1 + local idist,odist = 0,0 + local istride = invec.tda + local ostride = outvec.tda + local dim = ffi.new("int[?]", rank, dimlist) + local inembed = dim + local onembed = dim + local plan = ffi.gc(fftw.plan_many_dft_r2c(rank, dim, howmany, + inputtest,inembed,istride, idist, + output,onembed, ostride, odist, + bit.bor(fftw.MEASURE, fftw.UNALIGNED)), fftw.destroy_plan) --Save the plan for next time if rdftplans[rank] == nil then rdftplans[rank] = {} end - rdftplans[rank][direction] = {plan=plan, dimlist=dimlist} + rdftplans[rank][direction] = {plan=plan, dimlist=dimlist, itda=invec.tda, otda=outvec.tda} --Execute the plan with the supplied input fftw.execute_dft_r2c(plan, input, output) @@ -132,7 +157,11 @@ local function rdftinv(invec, dimlist, outvec) local input = ffi.cast("fftw_complex*",invec.data) --Check if plan is existing, then use it - if rdftplans[rank] ~= nil and rdftplans[rank][direction] ~= nil and compare_dimensions(dimlist, rdftplans[rank][direction].dimlist) then + if rdftplans[rank] ~= nil and rdftplans[rank][direction] ~= nil and + compare_dimensions(dimlist, rdftplans[rank][direction].dimlist) and + invec.tda == rdftplans[rank][direction].itda and + outvec.tda == rdftplans[rank][direction].otda then + local plan = rdftplans[rank][direction].plan fftw.execute_dft_c2r(plan, input, output) if not output_supplied then return outvec end @@ -142,12 +171,21 @@ local function rdftinv(invec, dimlist, outvec) --allcoate input for plan making but planning can invalidate input, therefore we use extra array local inputtest = ffi.new("fftw_complex[?]", inputsize) - --Create plan via MEASURE - local plan = ffi.gc(fftw.plan_dft_c2r(rank, ffi.new("int[?]", rank, dimlist), inputtest, output, bit.bor(fftw.MEASURE, fftw.UNALIGNED)), fftw.destroy_plan) + local howmany = 1 + local idist,odist = 0,0 + local istride = invec.tda + local ostride = outvec.tda + local dim = ffi.new("int[?]", rank, dimlist) + local inembed = dim + local onembed = dim + local plan = ffi.gc(fftw.plan_many_dft_c2r(rank, dim, howmany, + inputtest, inembed, istride, idist, + output, onembed, ostride, odist, + bit.bor(fftw.MEASURE, fftw.UNALIGNED)), fftw.destroy_plan) --Save the plan for next time if rdftplans[rank] == nil then rdftplans[rank] = {} end - rdftplans[rank][direction] = {plan=plan, dimlist=dimlist} + rdftplans[rank][direction] = {plan=plan, dimlist=dimlist, itda=invec.tda, otda=outvec.tda} --Execute the plan with the supplied input fftw.execute_dft_c2r(plan, input, output) From ecf844ffcf950324bb7f6009d9c649a5cac06108 Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Fri, 21 Feb 2014 11:22:07 +0100 Subject: [PATCH 13/20] Removed input/output dimensions because NULL suffices there --- fft.lua | 18 ++++++------------ 1 file changed, 6 insertions(+), 12 deletions(-) diff --git a/fft.lua b/fft.lua index 93b69d222..b78d36cfa 100644 --- a/fft.lua +++ b/fft.lua @@ -54,11 +54,9 @@ local function dft(invec, dimlist, outvec, direction) local istride = invec.tda local ostride = outvec.tda local dim = ffi.new("int[?]", rank, dimlist) - local inembed = dim - local onembed = dim local plan = ffi.gc(fftw.plan_many_dft(rank, dim ,howmany, - inputtest,inembed, istride, idist, - output,onembed,ostride,odist, + inputtest,nil, istride, idist, + output,nil,ostride,odist, direction == FORWARD and fftw.FORWARD or fftw.BACKWARD, bit.bor(fftw.MEASURE, fftw.UNALIGNED)), fftw.destroy_plan) --Save the plan for next time @@ -116,11 +114,9 @@ local function rdft(invec, dimlist, outvec) local istride = invec.tda local ostride = outvec.tda local dim = ffi.new("int[?]", rank, dimlist) - local inembed = dim - local onembed = dim local plan = ffi.gc(fftw.plan_many_dft_r2c(rank, dim, howmany, - inputtest,inembed,istride, idist, - output,onembed, ostride, odist, + inputtest,nil,istride, idist, + output,nil, ostride, odist, bit.bor(fftw.MEASURE, fftw.UNALIGNED)), fftw.destroy_plan) --Save the plan for next time @@ -176,11 +172,9 @@ local function rdftinv(invec, dimlist, outvec) local istride = invec.tda local ostride = outvec.tda local dim = ffi.new("int[?]", rank, dimlist) - local inembed = dim - local onembed = dim local plan = ffi.gc(fftw.plan_many_dft_c2r(rank, dim, howmany, - inputtest, inembed, istride, idist, - output, onembed, ostride, odist, + inputtest, nil, istride, idist, + output, nil, ostride, odist, bit.bor(fftw.MEASURE, fftw.UNALIGNED)), fftw.destroy_plan) --Save the plan for next time From ef519cebe74df3ecb6852b9b2f244bf04edff6ba Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Tue, 25 Feb 2014 11:50:12 +0100 Subject: [PATCH 14/20] Fixed seg fault issue when exiting gsl-shell due to allocation of too few memory for FFTW input plans --- fft.lua | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/fft.lua b/fft.lua index b78d36cfa..b707a4441 100644 --- a/fft.lua +++ b/fft.lua @@ -47,7 +47,7 @@ local function dft(invec, dimlist, outvec, direction) --Or create a new plan else --allcoate input for plan making but planning can invalidate input, therefore we use extra array - local inputtest = ffi.new("fftw_complex[?]", size) + local inputtest = ffi.new("fftw_complex[?]", size * invec.tda) local howmany = 1 local idist,odist = 0,0 @@ -106,8 +106,7 @@ local function rdft(invec, dimlist, outvec) --Or create a new plan else --allcoate input for plan making but planning can invalidate input, therefore we use extra array - local inputtest = ffi.new("double[?]", size) - -- local inputtest = ffi.gc(fftw.malloc(size * ffi.sizeof(double)), fftw.free) + local inputtest = ffi.new("double[?]", size * invec.tda) local howmany = 1 local idist,odist = 0,0 @@ -165,7 +164,7 @@ local function rdftinv(invec, dimlist, outvec) --Or create a new plan else --allcoate input for plan making but planning can invalidate input, therefore we use extra array - local inputtest = ffi.new("fftw_complex[?]", inputsize) + local inputtest = ffi.new("fftw_complex[?]", inputsize * invec.tda) local howmany = 1 local idist,odist = 0,0 From cc10d928b317d7a81b9eaf0cf07af19eb8abeefd Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Mon, 26 May 2014 17:23:25 +0200 Subject: [PATCH 15/20] Removed fft-init from Makefile and added fftw relevant Makefile entries. --- Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Makefile b/Makefile index 3a3c2cb3f..bca4113f4 100644 --- a/Makefile +++ b/Makefile @@ -58,7 +58,7 @@ ifeq ($(strip $(USE_READLINE)),yes) C_SRC_FILES += completion.c endif -LUA_BASE_FILES = bspline.lua fft-init.lua integ-init.lua template.lua check.lua \ +LUA_BASE_FILES = bspline.lua fft.lua fftw3/init.lua fftw3/cdefs.lua fftw3/defines.lua integ-init.lua template.lua check.lua \ graph-init.lua rng.lua rnd.lua randist.lua iter.lua time.lua gsl-check.lua linfit.lua \ roots.lua contour.lua gsl.lua matrix.lua csv.lua gslext.lua num.lua demo-init.lua \ import.lua plot3d.lua sf.lua vegas.lua eigen.lua help.lua cgdt.lua expr-actions.lua \ From ea0a8b777666295c0e16d1a8ebd790a0e08bc4a3 Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Tue, 27 May 2014 11:02:53 +0200 Subject: [PATCH 16/20] Changed documentation to hold information about all available fft functions: fft, fftinv, rfft, rfftinv, fft2, fft2inv, rfft2, rfft2inv, fftn, fftninv, rfftn, rfftninv. --- doc/user-manual/fft.rst | 130 +++++++++++++++++++--------------------- 1 file changed, 60 insertions(+), 70 deletions(-) diff --git a/doc/user-manual/fft.rst b/doc/user-manual/fft.rst index 5bdce72ff..00f9ecc1a 100644 --- a/doc/user-manual/fft.rst +++ b/doc/user-manual/fft.rst @@ -31,106 +31,96 @@ and the definition of the inverse Fourier transform is, The factor of 1/n makes this a true inverse. -In general there are two possible choices for the sign of the exponential in the transform/ inverse-transform pair. -GSL follows the same convention as fftpack, using a negative exponential for the forward transform. +In general there are two possible choices for the sign of the exponential in the transform/ inverse-transform pair. GSL follows the same convention as FFTW, using a negative exponential for the forward transform. The advantage of this convention is that the inverse transform recreates the original function with simple Fourier synthesis. Numerical Recipes uses the opposite convention, a positive exponential in the forward transform. -GSL Shell interface -------------------- -GSL Shell provide a simple interface to perform Fourier transforms of real data with the functions :func:`num.fft` and :func:`num.fftinv`. -The first function performs the Fourier transform of a column matrix and the second is the inverse Fourier transform. - -The function :func:`num.fft` returns a half-complex array. -This latter is similar to a column matrix of complex numbers, but it is actually a different object because the numbers are packed together following some specific rules related to the algorithm. +Methods +------------------------------ -The idea is that you can access the elements of this vector for reading or writing simply by indexing it. -You can also obtain the size of the vector using the operator '#'. -The valid indices for a half-complex object range from 0 to N-1 where N is the size if the vector. -Each element of the vector corresponds to the coefficient :math:`z_k` defined above. +The FFT implementation in gsl-shell is an interface to the `FFTW library `_ which is well tested and very fast for repetitive transformations. -When performing Fourier transforms, it is important to know that the computation speed can be greatly influenced by the size of the vector. If the size is a power of two, a very efficient algorithm can be used and we can talk in this case of a Fast Fourier Transform (FFT). In addition, the algorithm has the advantage that it does not require any additional workspace. When the size of the vector is not a power of two, we can have two different cases: +.. function:: fft(input [, output]) - * the size is a product of small prime numbers - * the size contains a big (> 7) prime number in its factorization + Performs the Fourier transform of the complex-valued column matrix ``input``. + The transformed signal is stored in a complex-valued column matrix ``output`` which is allocated by the function each time. To speed up repetitive Fourier transformations it is recommended to give the output column matrix as a second argument to eliminate the allocation overhead. -This detail is important because if the size is a product of small prime numbers, a fast algorithm is still available but it is still somewhat slower and it does require some additional workspace. -In the worst case when the size cannot be factorized to small prime numbers, the Fourier transform can still be computed but the calculation is slower, especially for large arrays. +.. function:: fftinv(input [, output]) -GSL Shell hides all the details and takes care of choosing the appropriate algorithm based on the size of the vector. -It also transparently provides any additional workspace that may be needed for the algorithm. -In order to avoid repeated allocation of workspace memory, the workspace allocated is kept in memory and reused *if the size of the array does not change*. -This means that the approach of GSL Shell is quite optimal if you perform many Fourier transforms (direct or inverse) of the same size. + Performas the inverse Fourier transformation of the complex-valued column matrix ``input``. The output is stored in a complex-valued column matrix which is allocated by the function each time. To speed up repetitive Fourier transformations it is recommended to give the output column matrix as a second argument to eliminate the allocation overhead. -Even though GSL Shell takes care of the details automatically, you should be aware of these performance notices because it can make a big difference in real applications. -From a practical point of view, it is useful in most cases to always provide samples whose size is a power of two. + This transformation is the inverse of the function :func:`num.fft`, so that if you perform the two transformations consecutively you will obtain a vector identical to the initial one up to a factor ``size=#input``. -Another property of the functions :func:`num.fft` and :func:`num.fftinv` is that they can optionally perform the transformation *in place* by modifying the original data instead of creating a copy. -When a transformation *in place* is requested, the routine still returns a new vector (either a real matrix or a half-complex array) but this latter will point to the same underlying data of the original vector. -The transformation *in place* can be useful in some cases to avoid unnecessary data copying and memory allocation. + A typical usage of :func:`fftinv` is to revert the transformation made with :func:`fft` but by doing some transformations along the way. + So a typical usage path could be:: + -- we assume v is a column matrix with our data + ft = num.fft(v) -- Fourier transform -Fourier Transform of Real Data ------------------------------- + -- here we can manipulate the half-complex array 'ft' + -- using the methods `get' and `set' + {some code here} -For real data, the Fourier coefficients satisfy the relation + vt = num.fftinv(ft)/#ft -- we perform the inverse Fourier transform + -- now vt is a vector of the same size of v -.. math:: - z_k = z_{N-k}^* +.. function:: rfft(input[, output]) -where N is the size of the vector and k is any integer number from 0 to N-1. -Because of this relation, the data is packed in a special type of object called a half-complex array. + Performs the Fourier transformation of real-valued input and complex-valued output. + In the output, the Fourier coefficients satisfy the relation -To access an element in a half-complex array, you can index it with an integer number between 0 and N-1, inclusive. So, for example:: + .. math:: + z_k = z_{N-k}^* - -- get a random number generator - r = rng.new() + where :math:`N` is the size of the vector and :math:`k` is any integer number from 0 to :math:`N-1`. + As a result of this symmetry, half of the output :math:`z` is redundant (being the complex conjugate of the other half), and so the transform only outputs elements :math:`0...\frac{n}{2}` of :math:`z` (:math:`\frac{n}{2}+1` complex numbers), where the division by 2 is rounded down. - -- create a vector with random numbers - x = matrix.new(256, 1, || rnd.gaussian(r, 1)) +.. function:: rfftinv(input[, output]) - -- take the Fourier transform - ft = num.fft(x) + Performs the inverse Fourier transformation with complex-valued input and real-valued output. + For the input of :math:`n`, the output has size :math:`(n-1)\times2`. + This is the direct inversion of :func:`num.rfft` as see in the example:: - -- print all the coefficients of the Fourier transform - for k=0, #ft-1 do print(ft[k]) end + --Forward transformation + ft = num.rfft(matrix.vec{1,2,3,4}) -As shown in the example above, you can use the Lua operator '#' to obtain the size of a half-complex array. + --Backward transformation + orig = num.rfftinv(ft) / #ft -.. function:: fft(v[, in_place]) +.. function:: fft2(input, [output]) +.. function:: fft2inv(input, [output]) - Perform the Fourier transform of the real-valued column matrix ``x``. - If ``in_place`` is ``true`` then the original data is altered and the resulting vector will point to the same underlying data of the original vector. + Performs a 2D forward and backward Fourier transformation of an input matrix ``input``. Giving a preallocated output matrix as a second argument speeds up repetitive transformations. - Please note that the value you obtain is not an ordinary matrix but a half-complex array. - You can access the elements of such an array by indexing the vector. - If you want to have an ordinary matrix you can easily build it with the following instructions:: +.. function:: rfft2(input, [output]) + + Performs a 2D forward Fourier transformation with real-valued input matrix ``input``. Returns the complex-valued output with reduced dimension. Giving a preallocated output matrix as a second argument speeds up repetitive transformations. - -- we suppose that f is an half-complex array - m = matrix.cnew(#f, 1, |i,j| f[i-1]) +.. function:: rfft2inv(input, [output]) -.. function:: fftinv(hc[, in_place]) + Performs the 2D inverse Fourier transformation with complex-valued input matrix ``input``. + Returns the real-valued output with increased dimension. Giving a preallocated output matrix as a second argument speeds up repetitive transformations. - Return a column matrix that contains the inverse Fourier transform of the half-complex vector ``hc``. - If ``in_place`` is ``true`` then the original data is altered and the resulting vector will point to the same underlying data of the original vector. +.. function:: fftn(input, dimlist, [output]) +.. function:: fftninv(input, dimlist, [output]) - This transformation is the inverse of the function :func:`num.fft`, so that if you perform the two transformations consecutively you will obtain a vector identical to the initial one. + Performs the n-dimensional forward and backward Fourier transformation of input data ``input``. + The input is considered to be stored in a `row-major format` meaning that if you had an array with dimensions :math:`n_1 \times n_2 \times n_3` then a point at :math:`(i_1,i_2,i_3)` is stored at index position: :math:`i_3 + n_3\cdot(i_2 + n_2\cdot i_1)`. + The dimensions are given as a table ``dimlist`` as :math:`{n1, n2, n3, ...}`. - A typical usage of :func:`fft_inv` is to revert the transformation made with :func:`fft` but by doing some transformations along the way. - So a typical usage path could be:: + As before, providing the output array as well limits the overhead of allocating the output array each call. - -- we assume v is a column matrix with our data - ft = num.fft(v) -- Fourier transform +.. function:: rfftn(input, dimlist, [output]) +.. function:: rfftninv(input, dimlist, [output]) + + Performs the n-dimensional real forward and backward Fourier transformation. Here the size of the input and output arrays are similar to the 1D and 2D couterparts. - -- here we can manipulate the half-complex array 'ft' - -- using the methods `get' and `set' - some code here + For forward transformations the real-valued input has size :math:`n_1 \times n_2 \times ... \times n_N` and the complex-valued output is of size :math:`n_1 \times n_2 \times ... \times n_N/2+1`. - vt = num.fftinv(ft) -- we perform the inverse Fourier transform - -- now vt is a vector of the same size of v + For the backward transformation the complex-valued input has size :math:`n_1 \times n_2 \times ... \times n_N/2+1` and the real-valued output has size :math:`n_1 \times n_2 \times ... \times n_N` consequently. -FFT example +Examples ----------- In this example we will treat a square pulse in the temporal domain. To illustrate a typical example of FFT usage we perform the Fourier Transform of the signal and we cut the higher order frequencies. Then we perform the inverse transform and we compare the result with the original time signal. @@ -157,13 +147,13 @@ Now we are ready to perform: and plot the results:: - ft = num.fft(y) + ft = num.rfft(y) - pf = graph.fibars(|k| complex.abs(ft[k]), 0, 60) + pf = graph.fibars(|k| complex.abs(ft[k]), 1, 60) pf.title = 'FFT Power Spectrum' - for k=ncut, n/2 do ft[k] = 0 end - ytr = num.fftinv(ft) + for k=ncut, #ft do ft[k] = 0 end + ytr = num.rfftinv(ft)/n pt:addline(graph.filine(|i| ytr[i], n), 'red') From ba2cc35ba22dcb1fb9c8966394fcba86b654d946 Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Tue, 27 May 2014 11:03:54 +0200 Subject: [PATCH 17/20] Added FFT note to the author contributions. --- doc/gsl-shell-index/authors.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/gsl-shell-index/authors.rst b/doc/gsl-shell-index/authors.rst index 04fc660b7..c3e5a0cbf 100644 --- a/doc/gsl-shell-index/authors.rst +++ b/doc/gsl-shell-index/authors.rst @@ -8,6 +8,6 @@ The main author of GSL Shell is **Francesco Abbate**. He was so fearless that he The VEGAS algorithm implementation is from **Lesley De Cruz**. Since he joined the project the author feels less lonely and knows that the user base of GSL Shell is of at least two people (including the authors). -The new implementation of the Special Function and Eigensystem modules are from **Benjamin von Ardenne**. The author suspect it is a pseudonym of Lesley so that the author have the impression of having a huge user base (three users). +The new implementation of the Special Function, Eigensystem and Fast Fourier Transformation modules are contributed by **Benjamin von Ardenne**. The author suspect it is a pseudonym of Lesley so that the author have the impression of having a huge user base (three users). The smart auto-complete function is based on the work of **Jesus Romero Hebrero**. Desperate about the bad English of Francesco he is currently striving to learn Italian and works on an integrated editor in the spare time. From 3d6a5142c17f6b42b4f2cd3e3d20e5b1c3f281b1 Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Tue, 27 May 2014 11:04:51 +0200 Subject: [PATCH 18/20] Minor format changes in fft.lua --- fft.lua | 4 ---- 1 file changed, 4 deletions(-) diff --git a/fft.lua b/fft.lua index b707a4441..cc4460660 100644 --- a/fft.lua +++ b/fft.lua @@ -102,7 +102,6 @@ local function rdft(invec, dimlist, outvec) local plan = rdftplans[rank][direction].plan fftw.execute_dft_r2c(plan, input, output) if not output_supplied then return outvec end - --Or create a new plan else --allcoate input for plan making but planning can invalidate input, therefore we use extra array @@ -121,7 +120,6 @@ local function rdft(invec, dimlist, outvec) --Save the plan for next time if rdftplans[rank] == nil then rdftplans[rank] = {} end rdftplans[rank][direction] = {plan=plan, dimlist=dimlist, itda=invec.tda, otda=outvec.tda} - --Execute the plan with the supplied input fftw.execute_dft_r2c(plan, input, output) if not output_supplied then return outvec end @@ -160,7 +158,6 @@ local function rdftinv(invec, dimlist, outvec) local plan = rdftplans[rank][direction].plan fftw.execute_dft_c2r(plan, input, output) if not output_supplied then return outvec end - --Or create a new plan else --allcoate input for plan making but planning can invalidate input, therefore we use extra array @@ -179,7 +176,6 @@ local function rdftinv(invec, dimlist, outvec) --Save the plan for next time if rdftplans[rank] == nil then rdftplans[rank] = {} end rdftplans[rank][direction] = {plan=plan, dimlist=dimlist, itda=invec.tda, otda=outvec.tda} - --Execute the plan with the supplied input fftw.execute_dft_c2r(plan, input, output) if not output_supplied then return outvec end From b7aa2525b3a583ba56ca399190fec90854198914 Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Tue, 27 May 2014 11:07:26 +0200 Subject: [PATCH 19/20] Added FFTW to the dependency in the INSTALL note. --- INSTALL | 1 + 1 file changed, 1 insertion(+) diff --git a/INSTALL b/INSTALL index a2ed2dcfc..dfde4ba82 100644 --- a/INSTALL +++ b/INSTALL @@ -9,6 +9,7 @@ GSL Shell needs the following external libraries: - GSL Library, version 1.14 or 1.15 - GNU readline library (on Linux only) - Freetype2 library, version 2.4.10 +- FFTW shared library, version 3.3+ To build the Graphical user interface you will need also the FOX 1.6 library. From b4a2858e6a6cb6868b2b9d96497590622e1a2fe0 Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Tue, 27 May 2014 11:51:29 +0200 Subject: [PATCH 20/20] Now copies the input matrix for all real inverse transformations because of strange side effect of all c2r fftw function calls (also documented, see function) --- doc/user-manual/fft.rst | 6 +++--- fft.lua | 18 +++++++++++++++++- 2 files changed, 20 insertions(+), 4 deletions(-) diff --git a/doc/user-manual/fft.rst b/doc/user-manual/fft.rst index 00f9ecc1a..3d5694757 100644 --- a/doc/user-manual/fft.rst +++ b/doc/user-manual/fft.rst @@ -78,7 +78,7 @@ The FFT implementation in gsl-shell is an interface to the `FFTW library