More progress in FRFFT

Author: roflmaostc

Project.toml

@@ -30,7 +30,8 @@ FractionalTransforms = "e50ca838-b4f0-4a10-ad18-4b920bf1ae5c"
 ImageTransformations = "02fcd773-0e25-5acc-982a-7f6622650795"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Test", "FractionalTransforms", "Random", "ImageTransformations", "Zygote"]
+test = ["Test", "TestImages", "FractionalTransforms", "Random", "ImageTransformations", "Zygote"]

README.md

@@ -28,9 +28,14 @@ The main features are:
 * array/image shifting (including noteworthy subpixel shifts)
 * array/image shearing
 * several tools like `ffts`, `ft`, `fftshift_view` etc. allowing simpler use with Fourier transforms
+* Chirp Z-Transform
+* Fractional Fourier Transform
 
 Have a look in the [examples folder](examples/) for interactive examples. The [documentation](https://bionanoimaging.github.io/FourierTools.jl/dev/) offers a quick overview.
 
+## FFTW Threading
+By default we set 4 Threads. Use `FFTW.set_num_threads(N)` to set `N` threads.
+
 
 ## Cite
 If you use this package in an academic work, please cite us!

examples/fractional_fourier_transform.jl

@@ -70,18 +70,24 @@ img = Float32.(testimage("resolution_test_512"));
 md"
 Fractional order
 
-$(@bind s Slider(-1.5:0.01:2, show_value=true))
+$(@bind s Slider(-3:0.05:3, show_value=true))
 "
 
 # ╔═╡ 7c445baa-d970-4954-a3dc-df828971bfd7
-[gray_show(log1p.(abs.(ft(img)))) gray_show(log1p.(abs.(sqrt(length(img)) .* frfft(img, s))))]
+[gray_show(log1p.(abs.(ft(img)))) gray_show(log1p.(abs.(sqrt(length(img)) .* frfft(img, s, shift=true))))]
 
 # ╔═╡ 1915c023-69cf-4d18-90cb-b47465dbef69
 begin
 	plot(log1p.(abs.(ft(img)[(end+begin)÷2+1,:] ./ sqrt(length(img)))))
 	plot!(log1p.(abs.(frfft(img, s)[(end+begin)÷2+1,:])))
 end
 
+# ╔═╡ 227ae9a3-9387-4ac3-b391-e2a78ce40d49
+begin
+	plot((real.(ft(img)[(end+begin)÷2+1,200:300] ./ sqrt(length(img)))))
+	plot!((real.(frfft(img, s)[(end+begin)÷2+1,200:300])))
+end
+
 # ╔═╡ abff911a-e10d-4311-955a-7afc4e0d344c
 md"## Fractional Fourier Transform on Vector
 Comparison with [FractionalTransforms.jl](https://github.com/SciFracX/FractionalTransforms.jl) roughly matches.
@@ -110,9 +116,10 @@ end
 # ╠═459649ed-ca70-426e-8273-97b146b5bcd5
 # ╟─4371cfbf-a3b3-45dc-847b-019994fbb234
 # ╠═d90b7f67-4166-44fa-aab7-de2c4f38fc00
-# ╟─24901666-4cc4-497f-a6ff-68c3e7ead629
+# ╠═24901666-4cc4-497f-a6ff-68c3e7ead629
 # ╠═7c445baa-d970-4954-a3dc-df828971bfd7
 # ╠═1915c023-69cf-4d18-90cb-b47465dbef69
+# ╠═227ae9a3-9387-4ac3-b391-e2a78ce40d49
 # ╟─abff911a-e10d-4311-955a-7afc4e0d344c
 # ╠═bae3c5b7-8964-493b-9e7b-d343e092219c
 # ╠═07d2b3b6-3584-4c64-9c4a-138beb3d6b88

src/fractional_fourier_transform.jl

@@ -3,16 +3,17 @@ export frfft
 """
     frfft(arr, p; shift=false, method=:garcia)
 
-Calculates the fractional Fourier transform (FRFFT) of the order `p` of `arr`.
+Calculates the fractional Fast Fourier transform (FRFFT) of the order `p` of `arr`.
+No `dims` argument is supported yet.
 
-If `shift=false` the fraction FFT is calculated around the first entry.
-If `shift=true` the FRFT
+If `shift=false` the FRFFT is calculated around the first entry.
+If `shift=true` the FRFFT is calculated aroound the center.
 
 
 ## Methods
 Several implementation exists. The following are implemented: 
 
-* `method=:garcia`: García, J., Mas, D., & Dorsch, R. G. (1996). Fractional-Fourier-transform calculation through the fast-Fourier-transform algorithm. Applied Optics, 35(35), 7013. doi:10.1364/ao.35.007013 
+* `method=:garcia`: A convolutional approach based on 2 FFTs. See García, J., Mas, D., & Dorsch, R. G. (1996). Fractional-Fourier-transform calculation through the fast-Fourier-transform algorithm. Applied Optics, 35(35), 7013. doi:10.1364/ao.35.007013 
 
 """
 function frfft(_arr, p; shift=true, method=:garcia)
@@ -44,7 +45,7 @@ function frfft(_arr, p; shift=true, method=:garcia)
 
     # trivial case
     if iszero(p)
-        return copy(_arr)
+        return complex.(_arr)
     end
 
     # handle shifting
@@ -105,7 +106,36 @@ function _frfft_garcia(arr::AbstractMatrix{T}, _p) where T
     t2 = -1im .* T(π) ./ size(arr) .* sin(p * T(π) / 2)
     φ2 = exp.(t2[1] .* m[1].^2 .+ t2[2] .* m[2]'.^2)
 
-    Mp = exp(-1im * T(π) * sign(sin(p * T(π) / 2)) / 4 + 1im * p * T(π) / 4 + 1im * T(π) /4)
+    Mp = exp(-1im * T(π) * sign(sin(p * T(π) / 2)) / 4 + 1im * p * T(π) / 4 + 1im * T(π) / 4)
+    # can do an in-place fft
+    res = Mp .* φ1 .* ifft!(φ2 .* fft!(arr .* φ1))
+    return res
+end
+
+
+function _frfft_garcia(arr::AbstractArray{T, N}, _p) where {N, T}
+    p = T(_p)
+    
+    
+    m = fftfreq.(size(arr), T.(size(arr)))
+
+    p1 = -1im .* T(π) ./ size(arr) .* tan(p * T(π) / 4)
+    t1 = similar(arr, complex(T))
+    fill!(t1, zero(T))
+
+
+    p2 = -1im .* T(π) ./ size(arr) .* sin(p * T(π) / 2)
+    t2 = similar(arr, complex(T))
+    fill!(t2, zero(T))
+    for d in 1:ndims(arr)
+        ns = ntuple(i -> i == d ? size(arr, i) : 1, Val(ndims(arr)))
+        t1 .+= reshape(p1[d] .* m[d].^2, ns...)
+        t2 .+= reshape(p2[d] .* m[d].^2, ns...)
+    end
+
+    φ1 = exp.(t1)
+    φ2 = exp.(t2)
+    Mp = exp(-1im * T(π) * sign(sin(p * T(π) / 2)) / 4 + 1im * p * T(π) / 4 + 1im * T(π) / 4)
     # can do an in-place fft
     res = Mp .* φ1 .* ifft!(φ2 .* fft!(arr .* φ1))
     return res

test/fractional_fourier_transform.jl

@@ -29,4 +29,17 @@
 
     @test all(.≈(real(frfft(frfft(box1d, 0.5, shift=true), -0.5, shift=true))[30:70] , real(box1d)[30:70], rtol=1e-1))
     @test all(.≈(real(frfft(frfft(box1d_, 0.5, shift=true), -0.5, shift=true))[30:70] , real(box1d_)[30:70], rtol=1e-1))
+
+
+
+    img = Float32.(testimage("resolution_test"))
+
+    @test abs.(ft(img)) ./ sqrt(length(img)) .+ 10 ≈ 10 .+ abs.(frfft(img, 0.999999)) rtol=1e-4
+    @test (real.(ft(img)) ./ sqrt(length(img)))[200:300] ≈ (real.(frfft(img, 0.999999)))[200:300] rtol=0.8
+
+    
+    x = randn((12,))
+    x2 = randn((13,))
+    @test frft(x, 0.5) ≈ frft(reshape(x, 12,1,1,1,1), 0.5)
+    @test frft(x, 0.5) ≈ reshape(frft(reshape(x, 1,12,1,1), 0.5), 12)
 end

test/runtests.jl

@@ -6,6 +6,7 @@ using Zygote
 using NDTools
 using LinearAlgebra # for the assigned nfft function LinearAlgebra.mul!
 using FractionalTransforms
+using TestImages
 
 Random.seed!(42)