add Dirichlet kernel diric (#412)

* add diric

Co-authored-by: Martin Holters <[email protected]>

Project.toml

@@ -1,6 +1,6 @@
 name = "DSP"
 uuid = "717857b8-e6f2-59f4-9121-6e50c889abd2"
-version = "0.7.2"
+version = "0.7.3"
 
 [deps]
 Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"

docs/src/util.md

@@ -23,4 +23,5 @@ shiftsignal
 shiftsignal!
 alignsignals
 alignsignals!
+diric
 ```

src/DSP.jl

@@ -20,6 +20,7 @@ include("periodograms.jl")
 include("Filters/Filters.jl")
 include("lpc.jl")
 include("estimation.jl")
+include("diric.jl")
 
 using Reexport
 @reexport using .Util, .Windows, .Periodograms, .Filters, .LPC, .Unwrap, .Estimation

src/diric.jl

@@ -0,0 +1,61 @@
+export diric
+
+"""
+    kernel = diric(Ω::Real, n::Integer)
+
+Dirichlet kernel, also known as periodic sinc function,
+where `n` should be a positive integer.
+Returns `sin(Ω * n/2) / (n * sin(Ω / 2))` typically,
+but `±1` when `Ω` is a multiple of 2π.
+
+In the usual case where 'n' is odd, the output is equivalent to
+``1/n \\sum_{k=-(n-1)/2}^{(n-1)/2} e^{i k Ω}``,
+which is the discrete-time Fourier transform (DTFT)
+of a `n`-point moving average filter.
+
+When `n` is odd or even, the function is 2π or 4π periodic, respectively.
+The formula for general `n` is
+`diric(Ω,n) = ``e^{-i (n-1) Ω/2}/n \\sum_{k=0}^{n-1} e^{i k Ω}``,
+which is a real and symmetric function of `Ω`.
+
+As of 2021-03-19, the Wikipedia definition has different factors.
+The definition here is consistent with scipy and other software frameworks.
+
+# Examples
+
+```jldoctest
+julia> round.(diric.((-2:0.5:2)*π, 5), digits=9)'
+1×9 adjoint(::Vector{Float64}) with eltype Float64:
+ 1.0  -0.2  0.2  -0.2  1.0  -0.2  0.2  -0.2  1.0
+
+julia> diric(0, 4)
+1.0
+```
+"""
+function diric(Ω::T, n::Integer) where T <: AbstractFloat
+    n > 0 || throw(ArgumentError("n=$n not positive"))
+    sign = one(T)
+    if isodd(n)
+        Ω = rem2pi(Ω, RoundNearest) # [-π,π)
+    else # wrap to interval [-π,π) to improve precision near ±2π
+        Ω = 2 * rem2pi(Ω/2, RoundNearest) # [-2π,2π)
+        if Ω > π # [π,2π)
+            sign = -one(T)
+            Ω -= T(2π) # [-π,0)
+        elseif Ω < -π # [-2π,-π)
+            sign = -one(T)
+            Ω += T(2π) # [0,π)
+        end
+    end
+
+    denom = sin(Ω / 2)
+    atol = eps(T)
+    if abs(denom) ≤ atol # denom ≈ 0 ?
+        return sign
+    end
+
+    return sign * sin(Ω * n/2) / (n * denom) # typical case
+end
+
+# handle non AbstractFloat types (e.g., Int, Rational, π)
+diric(Ω::Real, n::Integer) = diric(float(Ω), n)

test/diric.jl

@@ -0,0 +1,28 @@
+# diric.jl Dirichlet kernel tests
+
+@testset "diric" begin
+    @test_throws ArgumentError diric(0, -2)
+    @test @inferred(diric(0, 4)) ≈ 1
+    @test @inferred(diric(0 // 1, 5)) == 1
+    @test @inferred(diric(4π, 4)) ≈ 1
+    @test @inferred(diric(2π, 4)) ≈ -1
+    @test @inferred(diric(0, 5)) ≈ 1
+    @test @inferred(diric(2π, 5)) ≈ 1
+    @test @inferred(diric(π, 5)) ≈ 1 // 5
+    @test @inferred(diric(π/2, 5)) ≈ diric(-π/2, 5) ≈ -0.2
+    @test abs(@inferred(diric(π/2, 4))) < eps()
+    @test abs(@inferred(diric(3π/2, 4))) < eps()
+
+    # check type inference
+    @inferred diric(2, 4)
+    @inferred diric(Float16(1), 4)
+    @inferred diric(0.1f0, 4)
+    @inferred diric(0., 4)
+    @inferred diric(BigFloat(1), 4)
+
+    # accuracy test near 2π for even n
+    dirics(Ω,n) = real(exp(-1im * (n-1) * Ω/2)/n * sum(exp.(1im*(0:(n-1))*Ω)))
+    Ω = 2π .+ 4π * (-101:2:101)/100 * 1e-9
+    n = 400
+    @test dirics.(Ω,n) ≈ diric.(Ω,n)
+end

test/runtests.jl

@@ -1,8 +1,10 @@
 using DSP, FFTW, Test
+
 using DSP: allocate_output
 using Random: seed!
 
 testfiles = ["dsp.jl", "util.jl", "windows.jl", "filter_conversion.jl",
+    "diric.jl",
     "filter_design.jl", "filter_response.jl", "filt.jl", "filt_stream.jl",
     "periodograms.jl", "multitaper.jl", "resample.jl", "lpc.jl", "estimation.jl", "unwrap.jl",
     "remez_fir.jl" ]