1D Fermat Number Transform deconvolution

Project.toml

@@ -5,6 +5,7 @@ version = "0.3.0"
 
 [deps]
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+NumberTheoreticTransforms = "8497c1d1-af72-4391-8d22-bdd566511a1c"
 
 [compat]
 FFTW = "0.0, 0.1, 0.2, 0.3, 1"
@@ -13,6 +14,7 @@ julia = "1"
 [extras]
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+SpecialMatrices = "928aab9d-ef52-54ac-8ca1-acd7ca42c160"
 
 [targets]
-test = ["Test", "Random"]
+test = ["Test", "Random", "SpecialMatrices"]

src/Deconvolution.jl

@@ -17,5 +17,6 @@ using FFTW
 
 include("wiener.jl")
 include("lucy.jl")
+include("fermat.jl")
 
 end # module

src/fermat.jl

@@ -0,0 +1,39 @@
+
+using NumberTheoreticTransforms
+
+export fermat
+
+"""
+    fermat(convolved, h, g, q)
+
+Calculates deconvolution with Number Theoretic Transform.
+"""
+function fermat(convolved::AbstractArray{T, 1}, h::AbstractArray{T, 1}, g::T, q::T) where {T<:Integer}
+    N = length(convolved)
+    bias = div(q-1, 2) + 1 #negative numbers are represeted by the upper half of [0, q) range
+    H = fnt(h, g, q)
+    H_inv = invmod.(H, q)
+    D = fft(h) |> prod |> real |> T #TODO: avoid floats
+    Dm = mod(D, q)
+    xm = ifnt(H_inv .* fnt(convolved, g, q), g, q)
+    x0 = mod.(Dm * xm, q)
+    x0[x0 .>= bias] .-= q
+    x = x0
+    y_prim = fft(h) .* fft(x0) |> ifft .|> real .|> round.|> T #TODO: avoid floats
+    @assert all(iszero, rem.(convolved * D - y_prim, q))
+    y = div.(convolved * D - y_prim, q)
+    ym = mod.(y, q)
+    j = 1
+    while !all(iszero, y)
+        xj = ifnt(H_inv .* fnt(ym, g, q), g, q)
+        xj[xj .>= bias] .-= q
+        x = x .+ (xj * q^j)
+        y_prim = fft(h) .* fft(xj) |> ifft .|> real .|> round .|> T #TODO: avoid floats
+        @assert all(iszero, rem.(y - y_prim, q))
+        y = div.(y - y_prim, q)
+        ym = mod.(y, q)
+        j = j + 1
+    end
+
+    return x // D
+end

test/runtests.jl

@@ -50,3 +50,25 @@ end
 
     @test sum(abs.(s .- estimated)) < sum(abs.(s .- blurred_s))
 end
+
+##### Fermat Number Transform deconvolution
+
+@testset "Fermat Number Transform deconvolution tests" begin
+    using SpecialMatrices
+
+    h = [3, 2, 0, 0]
+    y = [3, 5, 3, 0]
+    x = fermat(y, h, 4, 17)
+    @test Circulant(h) * x == y
+
+    h = [4, 3, 1, 1, 1, 0, 0, 0]
+    y = [4, 2, 5, 1, 4, 2, 1, 3]
+    x = fermat(y, h, 2, 17)
+    @test Circulant(h) * x == y
+
+    h = [3, 2, 0, 0]
+    y = [11, 8, 13, 18]
+    x = fermat(y, h, 4, 17)
+    @test Circulant(h) * x == y
+
+end