|
| 1 | +export diric |
| 2 | + |
| 3 | +""" |
| 4 | + kernel = diric(Ω::Real, n::Integer) |
| 5 | +
|
| 6 | +Dirichlet kernel, also known as periodic sinc function, |
| 7 | +where `n` should be a positive integer. |
| 8 | +Returns `sin(Ω * n/2) / (n * sin(Ω / 2))` typically, |
| 9 | +but `±1` when `Ω` is a multiple of 2π. |
| 10 | +
|
| 11 | +In the usual case where 'n' is odd, the output is equivalent to |
| 12 | +``1/n \\sum_{k=-(n-1)/2}^{(n-1)/2} e^{i k Ω}``, |
| 13 | +which is the discrete-time Fourier transform (DTFT) |
| 14 | +of a `n`-point moving average filter. |
| 15 | +
|
| 16 | +When `n` is odd or even, the function is 2π or 4π periodic, respectively. |
| 17 | +The formula for general `n` is |
| 18 | +`diric(Ω,n) = ``e^{-i (n-1) Ω/2}/n \\sum_{k=0}^{n-1} e^{i k Ω}``, |
| 19 | +which is a real and symmetric function of `Ω`. |
| 20 | +
|
| 21 | +As of 2021-03-19, the Wikipedia definition has different factors. |
| 22 | +The definition here is consistent with scipy and other software frameworks. |
| 23 | +
|
| 24 | +# Examples |
| 25 | +
|
| 26 | +```jldoctest |
| 27 | +julia> round.(diric.((-2:0.5:2)*π, 5), digits=9)' |
| 28 | +1×9 adjoint(::Vector{Float64}) with eltype Float64: |
| 29 | + 1.0 -0.2 0.2 -0.2 1.0 -0.2 0.2 -0.2 1.0 |
| 30 | +
|
| 31 | +julia> diric(0, 4) |
| 32 | +1.0 |
| 33 | +``` |
| 34 | +""" |
| 35 | +function diric(Ω::T, n::Integer) where T <: AbstractFloat |
| 36 | + n > 0 || throw(ArgumentError("n=$n not positive")) |
| 37 | + sign = one(T) |
| 38 | + if isodd(n) |
| 39 | + Ω = rem2pi(Ω, RoundNearest) # [-π,π) |
| 40 | + else # wrap to interval [-π,π) to improve precision near ±2π |
| 41 | + Ω = 2 * rem2pi(Ω/2, RoundNearest) # [-2π,2π) |
| 42 | + if Ω > π # [π,2π) |
| 43 | + sign = -one(T) |
| 44 | + Ω -= T(2π) # [-π,0) |
| 45 | + elseif Ω < -π # [-2π,-π) |
| 46 | + sign = -one(T) |
| 47 | + Ω += T(2π) # [0,π) |
| 48 | + end |
| 49 | + end |
| 50 | + |
| 51 | + denom = sin(Ω / 2) |
| 52 | + atol = eps(T) |
| 53 | + if abs(denom) ≤ atol # denom ≈ 0 ? |
| 54 | + return sign |
| 55 | + end |
| 56 | + |
| 57 | + return sign * sin(Ω * n/2) / (n * denom) # typical case |
| 58 | +end |
| 59 | + |
| 60 | +# handle non AbstractFloat types (e.g., Int, Rational, π) |
| 61 | +diric(Ω::Real, n::Integer) = diric(float(Ω), n) |
0 commit comments