Merge pull request #11 from ranocha/fourier_pade

improve Fourier Padé

README.md

@@ -18,7 +18,8 @@ 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)`
+- `fourier_pade(u, degree_num, degree_den, num_output=length(u))`
   Compute the Fourier-Padé reconstruction of `u` with degrees
-  `(degree_num, degree_den)` and evaluate it at the points `x`,
+  `(degree_num, degree_den)` and evaluate it at `num_output`
+  equispaced points,
   cf. [Driscoll and Fornberg (2001) A Padé-based algorithm for overcoming the Gibbs phenomenon](https://doi.org/10.1023/A:1016648530648).

src/fourier_pade.jl

@@ -1,54 +1,46 @@
 """
-    fourier_pade(x, u, degree_num, degree_den)
+    fourier_pade(u, degree_num, degree_den, num_output=length(u))
 
 Compute the Fourier-Padé reconstruction of `u` with degrees
-`(degree_num, degree_den)` and evaluate it at the points `x`, cf.
+`(degree_num, degree_den)` and evaluate it at `num_output`
+equispaced points, 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
+function fourier_pade(u, degree_num, degree_den, num_output=length(u))
   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."))
+  modes = length(u) ÷ 2 + 1
+  if modes <= N
+    throw(ArgumentError("Cannot perform a Fourier-Padé reconstruction with degrees ($degree_num, $degree_den) given $(length(u)) real coefficients corresponding to $modes complex modes."))
   end
 
-  uh = fft(u) / N0
-  uhat = uh[1:N+1]
+  uhat = rfft(u)
+  uhat ./= length(u)
 
   # 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
+  row_den = zeros(eltype(col_den), degree_den+1)
+  row_den[1:min(degree_num+2,degree_den+1)] = uhat[degree_num+2:-1:max(1, degree_num+2-degree_den)]
   Z = nullspace(Matrix(Toeplitz(col_den, row_den)))
-  den_p = Z[:, end]
+  den_p = zeros(eltype(uhat), num_output)
+  den_p[1:degree_den+1] = 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)
+  num_p = zeros(eltype(uhat), num_output)
+  mul!(view(num_p, 1:degree_num+1), A, view(den_p, 1:degree_den+1))
 
   # 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)))
+  bfft_plan = plan_bfft(num_p)
+  2 .* real.((bfft_plan * num_p) ./ (bfft_plan * den_p))
 end
 
+
 function reference(::typeof(fourier_pade))
 """
 @article{driscoll2001pade,

test/fourier_pade_test.jl

@@ -11,21 +11,19 @@ using LinearAlgebra, DelimitedFiles
     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)
+    u_reconstructed = fourier_pade(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)
+    u_reconstructed = fourier_pade(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)
+    u_reconstructed = fourier_pade(u, 50, 52)
 
     N = length(u)
-    @test_throws ArgumentError fourier_pade(x, u[1:end-1], N-(N÷2)+1, N÷2)
+    @test_throws ArgumentError fourier_pade(u[1:end-1], N-(N÷2)+1, N÷2)
   end
 end