From 944d2c5b83f4214d386f68ae3b693191c5630c9f Mon Sep 17 00:00:00 2001 From: Benjamin von Ardenne Date: Wed, 5 Feb 2014 11:40:23 +0100 Subject: [PATCH 1/2] 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 2/2] 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}