Add convolution operator ApproxFun.jl#940 (#127)

* Add Convolution operator

* typo in description

* overload getindex instead of defaultgetindex

* overload getindex instead of ApproxFunBase.getindex

---------

Co-authored-by: Jilhg <[email protected]>

Author: Jilhg

src/ApproxFunFourier.jl

@@ -776,4 +776,6 @@ function show(io::IO, S::SpaceTypes)
     print(io, ")")
 end
 
+include("Convolution.jl")
+
 end #module

src/Convolution.jl

@@ -0,0 +1,66 @@
+export Convolution
+
+"""
+	Convolution(G)
+
+Represents a convolution operator with function `G`. When applied to a function `f` it computes `∫_D G(x-y) f(y) dy for x ∈ D`. For now `D` is a `PeriodigSegment(a,b)` with `a<b∈R` and `space(G)` is either `Laurent(D)` or `Fourier(D)`.
+"""
+
+abstract type Convolution{S,T} <: Operator{T} end
+
+struct ConcreteConvolution{S<:Space,T} <: Convolution{S,T}
+    G::Fun{S,T}
+end
+
+function Convolution(G::Fun{S,T}) where {S,T}
+    @assert isfinite(arclength(domain(G)))
+    ConcreteConvolution(G)
+end
+ApproxFunBase.domain(C::ConcreteConvolution) = ApproxFunBase.domain(C.G)
+ApproxFunBase.domainspace(C::ConcreteConvolution) = ApproxFunBase.space(C.G)
+ApproxFunBase.rangespace(C::ConcreteConvolution) = ApproxFunBase.space(C.G)
+ApproxFunBase.bandwidths(C::ConcreteConvolution) = error("Please implement convolution bandwidths on "*string(space(C.G)))
+getindex(C::ConcreteConvolution,k::Integer,j::Integer) = error("Please implement convolution getindex on "*string(space(C.G)))
+
+ApproxFunBase.bandwidths(C::ConcreteConvolution{Laurent{PeriodicSegment{R},T}}) where {R<:Real,T} = 0,0
+ApproxFunBase.bandwidths(C::ConcreteConvolution{Fourier{PeriodicSegment{R},T}}) where {R<:Real,T} = 1,1
+
+function getindex(C::ConcreteConvolution{Laurent{PeriodicSegment{R},T1},T2},k::Integer,j::Integer) where {R<:Real,T1,T2}
+    fourier_index::Integer = if isodd(k) div(k-1,2) else -div(k,2) end
+    if k == j && k ≤ ncoefficients(C.G)
+        return (exp(-2pi*1im/arclength(domain(C.G))*fourier_index*first(domain(C.G)))*arclength(domain(C.G))*C.G.coefficients[k])::T2
+    else
+        return zero(T2)
+    end
+end
+
+function getindex(C::ConcreteConvolution{Fourier{PeriodicSegment{R},T1},T2},k::Integer,j::Integer) where {R<:Real,T1,T2}
+    fourier_index::Integer = if isodd(k) div(k-1,2) else div(k,2) end
+    if k<1 || j<1 || ncoefficients(C.G)==0
+        return zero(T2)
+    elseif k == 1
+        if j==k
+            return (arclength(domain(C.G))*C.G.coefficients[1])::T2
+        else
+            return zero(T2)
+        end
+    elseif 2*fourier_index ≤ ncoefficients(C.G)
+        Gs = if 2*fourier_index ≤ ncoefficients(C.G) C.G.coefficients[2*fourier_index] else zero(T2) end # sine coefficient
+        Gc = if 2*fourier_index+1 ≤ ncoefficients(C.G) C.G.coefficients[2*fourier_index+1] else zero(T2) end # cosine coefficient
+        phase = 2pi/arclength(domain(C.G))*fourier_index*first(domain(C.G))
+        if iseven(k) && j==k
+            return (arclength(domain(C.G))*(Gc*cos(phase)-Gs*sin(phase))/2)::T2
+        elseif iseven(k) && j==k+1
+            return (arclength(domain(C.G))*(Gc*sin(phase)+Gs*cos(phase))/2)::T2
+        elseif isodd(k) && j==k
+            return (arclength(domain(C.G))*(Gc*cos(phase)-Gs*sin(phase))/2)::T2
+        elseif isodd(k) && j==k-1
+            return (arclength(domain(C.G))*(-Gc*sin(phase)-Gs*cos(phase))/2)::T2
+        else
+            return zero(T2)
+        end
+    else
+        return zero(T2)
+    end
+end
+

test/runtests.jl

@@ -806,3 +806,42 @@ end
         end
     end
 end
+
+@testset "Convolution" begin
+    @test bandwidths(Convolution(ones(Laurent())))==(0,0)
+    @test bandwidths(Convolution(ones(Fourier())))==(1,1)
+
+    atol=1e-6
+    rtol=1e-6
+
+    for S=[Laurent(-15..3),Fourier(-15..3)]
+        f1=ones(S)
+        f2=ones(S)
+        f1f2_theory=18*ones(S)
+        f1f2=Convolution(f1)*f2
+        @test ≈(f1f2,f1f2_theory,atol=atol,rtol=rtol)
+        @test ≈(DefiniteIntegral()*f1f2,(DefiniteIntegral()*f1)*(DefiniteIntegral()*f2),atol=atol,rtol=rtol)
+    end
+
+    for S=[Laurent(-2..3),Fourier(-2..3)]
+        f1=Fun(x->exp(-5*(x-0.5)^2),S)
+        f2=Fun(x->exp(-5*x^2),S)
+        f1f2_theory=Fun(x->exp(-5*(x-0.5)^2/2)*sqrt(2pi)/10*sqrt(5),S)
+        f1f2=Convolution(f1)*f2
+        f2f1=Convolution(f2)*f1
+        @test ≈(f1f2,f1f2_theory,atol=atol,rtol=rtol)
+        @test ≈(f2f1,f1f2_theory,atol=atol,rtol=rtol)
+        @test ≈(DefiniteIntegral()*f1f2,(DefiniteIntegral()*f1)*(DefiniteIntegral()*f2),atol=atol,rtol=rtol)
+        @test ≈(DefiniteIntegral()*f2f1,(DefiniteIntegral()*f1)*(DefiniteIntegral()*f2),atol=atol,rtol=rtol)
+    end
+    for S=[Laurent(0..2pi),Fourier(0..2pi)]
+        f1=Fun(x->exp(1im*x),S)
+        f2=Fun(x->exp(1im*x),S)
+        f1f2_theory=Fun(x->2pi*exp(1im*x),S)
+        f1f2=Convolution(f1)*f2
+        @test ≈(f1f2,f1f2_theory,atol=atol,rtol=rtol)
+        @test ≈(DefiniteIntegral()*f1f2,(DefiniteIntegral()*f1)*(DefiniteIntegral()*f2),atol=atol,rtol=rtol)
+         
+        @test Convolution(Fun(S,[]))*f1≈0
+    end
+end