Merge 1494d4c into 756f4a4

Project.toml

@@ -4,16 +4,20 @@ authors = ["Hendrik Ranocha <mail@ranocha.de>"]
 version = "0.1.0"
 
 [deps]
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
 
 [compat]
 OffsetArrays = "1.0"
 julia = "1.4"
 
 [extras]
+DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Random", "Test"]
+test = ["DelimitedFiles", "Random", "Test"]

README.md

@@ -18,3 +18,7 @@ This is work in progress. Currently, the following functions are implemented.
 - `group_sparse_total_variation_denoising(y::AbstractVector, λ::Number; group_size::Integer=1, max_iter::Integer=100)`
   Compute `max_iter` iterations of the algorithm described by
   [Selesnick and Chen (2013) Total variation denoising with overlapping group sparsity](https://doi.org/10.1109/ICASSP.2013.6638755).
+- `fourier_pade(x, u, degree_num, degree_den)`
+  Compute the Fourier-Padé reconstruction of `u` with degrees
+  `(degree_num, degree_den)` and evaluate it at the points `x`,
+  cf. [Driscoll and Fornberg (2001) A Padé-based algorithm for overcoming the Gibbs phenomenon](https://doi.org/10.1023/A:1016648530648).

src/Postprocessing.jl

@@ -1,13 +1,24 @@
 module Postprocessing
 
 using LinearAlgebra
+using FFTW
 using OffsetArrays
+using ToeplitzMatrices
 
-include("total_variation.jl")
+"""
+    reference(algorithm)
 
+Returns a citable reference for `algorithm`, e.g. a bibtex entry
+for a scientific article or a book.
+"""
+function reference end
 
-export total_variation_denoising, total_variation_denoising!,
-       group_sparse_total_variation_denoising, group_sparse_total_variation_denoising!
+include("total_variation.jl")
+include("fourier_pade.jl")
 
+export reference
+export total_variation, total_variation_denoising, total_variation_denoising!,
+       group_sparse_total_variation_denoising, group_sparse_total_variation_denoising!
+export fourier_pade
 
 end # module

src/fourier_pade.jl

@@ -0,0 +1,66 @@
+"""
+    fourier_pade(x, u, degree_num, degree_den)
+
+Compute the Fourier-Padé reconstruction of `u` with degrees
+`(degree_num, degree_den)` and evaluate it at the points `x`, cf.
+Driscoll and Fornberg (2001) A Padé-based algorithm for overcoming the Gibbs phenomenon,
+doi: 10.1023/A:1016648530648.
+"""
+function fourier_pade(x, u, degree_num, degree_den)
+  N0 = length(u)
+  if iseven(N0)
+    throw(ArgumentError("u has $N0 coeffficients but is assumed to have an odd number of coefficients."))
+  end
+  N = degree_num + degree_den
+  if N0 <= N
+    throw(ArgumentError("Cannot perform a Fourier-Padé reconstruction with degrees $degree_num and $degree_den given $N0 coefficients."))
+  end
+
+  uh = fft(u) / N0
+  uhat = uh[1:N+1]
+
+  # denominator
+  col_den = uhat[degree_num+2:N+1]
+  row_den = uhat[degree_num+2:-1:max(1, degree_num+2-degree_den)]
+  if degree_den > degree_num
+    row_den = vcat(row_den, fill(0, degree_den-degree_num-1))
+  end
+  Z = nullspace(Matrix(Toeplitz(col_den, row_den)))
+  den_p = Z[:, end]
+  den_p ./= den_p[findfirst(!iszero, den_p)]
+  den_m = conj.(den_p)
+
+  # numerator
+  col_num = uhat[1:degree_num+1]
+  col_num[1] /= 2
+  row_num = zeros(eltype(col_num), degree_den+1)
+  row_num[1] = col_num[1]
+  A = Toeplitz(col_num, row_num)
+  num_p = A * den_p
+  num_m = conj.(num_p)
+
+  # evaluate Fourier-Padé approximation
+  T = real(eltype(uhat))
+  n = length(x)
+  xx = range(zero(T), 2*one(T)*π*(n-1)/n, length=n)
+  x_p = @. exp(im * xx)
+  x_m = @. exp(-im * xx)
+
+  real.(evalpoly.(x_p, Ref(num_p)) ./ evalpoly.(x_p, Ref(den_p)) .+ evalpoly.(x_m, Ref(num_m)) ./ evalpoly.(x_m, Ref(den_m)))
+end
+
+function reference(::typeof(fourier_pade))
+"""
+@article{driscoll2001pade,
+  title={A {P}ad{\'e}-based algorithm for overcoming the {G}ibbs phenomenon},
+  author={Driscoll, Tobin A and Fornberg, Bengt},
+  journal={Numerical Algorithms},
+  volume={26},
+  number={1},
+  pages={77--92},
+  year={2001},
+  publisher={Springer},
+  doi={10.1023/A:1016648530648}
+}
+"""
+end

src/total_variation.jl

@@ -1,3 +1,23 @@
+"""
+    total_variation(u)
+
+Compute the discrete total variation
+```math
+    \\sum_{k} |u_{k+1} - u_{k}|.
+```
+"""
+function total_variation(u)
+  uprev = first(u)
+  res = zero(uprev)
+  @inbounds @simd for i in eachindex(u)[2:end]
+    unext = u[i]
+    res += abs(uprev - unext)
+    uprev = unext
+  end
+  res
+end
+
+
 """
     total_variation_denoising(y::AbstractVector, λ::Number)
 
@@ -110,6 +130,23 @@ function total_variation_denoising!(dest::AbstractVector, source::AbstractVector
   dest
 end
 
+function reference(::Union{typeof(total_variation_denoising),
+                           typeof(total_variation_denoising!)})
+"""
+@article{condat2013direct,
+  title={A direct algorithm for 1-{D} total variation denoising},
+  author={Condat, Laurent},
+  journal={IEEE Signal Processing Letters},
+  volume={20},
+  number={11},
+  pages={1054--1057},
+  year={2013},
+  publisher={IEEE},
+  doi={10.1109/LSP.2013.2278339}
+}
+"""
+end
+
 
 
 """
@@ -189,3 +226,19 @@ function group_sparse_total_variation_denoising!(dest::AbstractVector, source::A
 
   dest
 end
+
+function reference(::Union{typeof(group_sparse_total_variation_denoising),
+                           typeof(group_sparse_total_variation_denoising!)})
+"""
+@inproceedings{selesnick2013total,
+  title={Total variation denoising with overlapping group sparsity},
+  author={Selesnick, Ivan W and Chen, Po-Yu},
+  booktitle={2013 IEEE International Conference on Acoustics,
+             Speech and Signal Processing},
+  pages={5696--5700},
+  year={2013},
+  organization={IEEE},
+  doi={10.1109/ICASSP.2013.6638755}
+}
+"""
+end

test/fourier_pade_test.jl

@@ -0,0 +1,31 @@
+using Postprocessing
+using Test
+using LinearAlgebra, DelimitedFiles
+
+@testset "Fourier Padé reconstruction" begin
+  println(devnull, reference(fourier_pade))
+
+  for T in (Float32, Float64)
+    data = readdlm(joinpath(@__DIR__, "test_u_piecewise_constant.txt"), comments=true)
+    x = T.(data[:, 1])
+    u = T.(data[:, 2])
+    u_true = @. ifelse(abs(x - 0.6) < 0.2, one(T), zero(T))
+
+    u_reconstructed = fourier_pade(x, u, 50, 50)
+    @test norm(u_reconstructed - u_true) < norm(u - u_true)
+    @test total_variation(u_reconstructed) < total_variation(u)
+
+    u_reconstructed = fourier_pade(x, u, 50, 45)
+    @test norm(u_reconstructed - u_true) < norm(u - u_true)
+    @test total_variation(u_reconstructed) < total_variation(u)
+
+    # this choice of the degrees leads to increased oscillations
+    # near the discontinuity but should still be computable
+    u_reconstructed = fourier_pade(x, u, 50, 52)
+
+    @test_throws ArgumentError fourier_pade(x, u[1:end-1], 50, 50)
+
+    N = length(u)
+    @test_throws ArgumentError fourier_pade(x, u[1:end-1], N-(N÷2)+1, N÷2)
+  end
+end

test/runtests.jl

@@ -1,2 +1,3 @@
 
 @time include("total_variation_test.jl")
+@time include("fourier_pade_test.jl")