remez: New API (WIP)

src/Filters/remez_fir.jl

@@ -118,16 +118,10 @@ end
 
 
 """
-build_grid(numtaps, bands, desired, weight, grid_density)
+build_grid(numtaps, band_defs, Hz, grid_density, neg)
 return "grid" and "des" and "wt" arrays
 """
-function build_grid(numtaps, bands, eff, wate, grid_density, neg)
-    # translate from scipy remez argument names
-    lgrid = grid_density
-
-    nbands = length(eff)
-    nfilt = numtaps
-
+function build_grid(nfilt, band_defs, Hz, grid_density, neg)
     nodd = isodd(nfilt)
     nfcns = nfilt ÷ 2
     if nodd && !neg
@@ -138,20 +132,40 @@ function build_grid(numtaps, bands, eff, wate, grid_density, neg)
     # SET UP THE DENSE GRID. THE NUMBER OF POINTS IN THE GRID
     # IS (FILTER LENGTH + 1)*GRID DENSITY/2
     #
-    delf = lgrid * nfcns
-    delf = 0.5 / delf
-
-    # calculate clamped band-edges
-    edges = reshape(bands, 2, nbands)
-    if neg || !nodd
-        flimlow = neg ? delf : 0.0
-        flimhigh = (neg == nodd) ? 0.5 - delf : 0.5
-        edges = map(f -> clamp(f, flimlow, flimhigh), edges)
+    delf = 0.5 / (grid_density * nfcns)
+
+    flimlow = neg ? delf : 0.0
+    flimhigh = (neg == nodd) ? 0.5 - delf : 0.5
+    function normalize_banddef_entry(b)
+        # normalize and clamp band-edges
+        fl = convert(Float64, clamp(b.first[1] / Hz, flimlow, flimhigh))
+        fu = convert(Float64, clamp(b.first[2] / Hz, flimlow, flimhigh))
+        # make sure desired and weight are functions
+        if !isa(b.second, Tuple{Any, Any})
+            desired = b.second
+            weight = 1.0
+        else
+            desired = b.second[1]
+            weight = b.second[2]
+        end
+        if isa(desired, Real)
+            let d = desired
+                desired = _ -> d
+            end
+        end
+        if isa(weight, Real)
+            let w = weight
+                weight = _ -> w
+            end
+        end
+        return Pair{Tuple{Float64,Float64},Tuple{Any,Any}}((fl, fu), (desired, weight))
     end
+    normalized_band_defs = normalize_banddef_entry.(band_defs)
+
     local ngrid
     # work around JuliaLang/julia#15276
-    let delf=delf, edges=edges
-        ngrid = sum(max(length(edges[1,lband]:delf:edges[2,lband]), 1) for lband in 1:nbands)#::Int
+    let delf=delf
+        ngrid = sum(max(length(band_def.first[1]:delf:band_def.first[2]), 1) for band_def in normalized_band_defs)#::Int
     end
 
     grid = zeros(Float64, ngrid)  # the array of frequencies, between 0 and 0.5
@@ -169,22 +183,29 @@ function build_grid(numtaps, bands, eff, wate, grid_density, neg)
     # SET UP A NEW APPROXIMATION PROBLEM WHICH IS EQUIVALENT
     # TO THE ORIGINAL PROBLEM
     #
-    for lband in 1:nbands
-        flow = edges[1, lband]
-        fup = edges[2, lband]
-        for f in [(flow:delf:fup)[1:end-1]; fup]
-            change = neg ? (nodd ? sinpi(2f) : sinpi(f)) : (nodd ? 1.0 : cospi(f))
-            grid[j] = cospi(2f)
-            des[j] = eff[lband](f) / change
-            wt[j] = wate[lband](f) * change
-            j += 1
-        end
+    for band_def in normalized_band_defs
+        flow = band_def.first[1]
+        fup = band_def.first[2]
+        # outline inner loop to have it type-stable (band_def.second is Tuple{Any,Any})
+        j = _buildgrid!(grid, des, wt, j, [(flow:delf:fup)[1:end-1]; fup], neg, nodd, Hz,
+                        band_def.second)
     end
     @assert ngrid == j - 1
 
     return grid, des, wt
 end
 
+function _buildgrid!(grid, des, wt, j, fs, neg, nodd, Hz, des_wt)
+    for f in fs
+        change = neg ? (nodd ? sinpi(2f) : sinpi(f)) : (nodd ? 1.0 : cospi(f))
+        grid[j] = cospi(2f)
+        des[j] = des_wt[1](f * Hz) / change
+        wt[j] = des_wt[2](f * Hz) * change
+        j += 1
+    end
+    return j
+end
+
 """
 function freq_eval(xf, x::AbstractVector, y::AbstractVector, ad::AbstractVector)
 -----------------------------------------------------------------------
@@ -332,38 +353,39 @@ grid()
 function remez(numtaps::Integer, bands::Vector, desired::Vector;
                weight::Vector=fill(1.0, length(desired)),
                Hz::Real=1.0,
-               filter_type::Union{RemezFilterType,Nothing}=nothing,
-               neg::Bool=filter_type in (filter_type_hilbert, filter_type_differentiator),
-               maxiter::Integer=25,
-               grid_density::Integer=16)
-    bands = bands/Hz
-
-    # Sanity checks on arguments
+               filter_type::RemezFilterType=filter_type_bandpass,
+               kwargs...)
     issorted(bands) || throw(ArgumentError("`bands` is not monotonically increasing"))
-    (0 <= bands[1]) && (bands[end] <= 0.5) ||
-      throw(ArgumentError("`bands` values must be between 0 and `Hz`/2"))
-    length(bands) == 2length(desired) || 
-      throw(ArgumentError("`desired` must be half the length of `bands`."))
-    length(bands) == 2length(weight) || 
-      throw(ArgumentError("`weight` must be half the length of `bands`."))
-    filter_type === nothing || all(d -> isa(d, Number), desired) ||
-      throw(ArgumentError("`filter_type` given but non-numeric entry in `desired`"))
-    filter_type === nothing || all(w -> isa(w, Number), weight) ||
-      throw(ArgumentError("`filter_type` given but non-numeric entry in `weight`"))
-    filter_type === nothing || neg == (filter_type !== filter_type_bandpass) ||
-      throw(ArgumentError("`neg` does not match given `filter_type`"))
-
-    bands = convert(Vector{Float64}, bands)   # in C, known as "edge"
-
+    length(bands) == 2length(desired) ||
+        throw(ArgumentError("`desired` must be half the length of `bands`."))
+    length(bands) == 2length(weight) ||
+        throw(ArgumentError("`weight` must be half the length of `bands`."))
+    band_ranges = [(bands[i], bands[i+1]) for i in 1:2:length(bands)]
     if filter_type == filter_type_differentiator
-        eff = [f -> d*f for d in desired]
-        wate = [d > 0.0001 ? f -> w/f : f -> w for (w, d) in zip(weight, desired)]
+        eff = [f -> d*f/Hz for d in desired]
+        wate = [d > 0.0001 ? f -> w/(f/Hz) : f -> w for (w, d) in zip(weight, desired)]
     else
-        eff = [isa(d, Number) ? f->d : d for d in desired]
-        wate = [isa(w, Number) ? f->w : w for w in weight]
+        eff = desired
+        wate = weight
     end
+    band_defs = [r => (d, w) for (r, d, w) in zip(band_ranges, eff, wate)]
+    neg = filter_type in (filter_type_hilbert, filter_type_differentiator)
+    return remez(numtaps, band_defs; Hz=Hz, neg=neg, kwargs...)
+end
+
+function remez(numtaps::Integer, band_defs;
+                Hz::Real=1.0,
+                neg::Bool=false,
+                maxiter::Integer=25,
+                grid_density::Integer=16)
+    all(b.first[1] <= b.first[2] for b in band_defs) ||
+        throw(ArgumentError("lower band edge higher then upper band edge"))
+    all(band_defs[i].first[2] < band_defs[i+1].first[1] for i in 1:length(band_defs)-1) ||
+        throw(ArgumentError("band edges is not monotonically increasing"))
+    (0 <= band_defs[1].first[1]) && (band_defs[end].first[2] <= 0.5*Hz) ||
+        throw(ArgumentError("band edges must be between 0 and `Hz`/2"))
 
-    grid, des, wt = build_grid(numtaps, bands, eff, wate, grid_density, neg)
+    grid, des, wt = build_grid(numtaps, band_defs, Hz, grid_density, neg)
 
     nfilt = numtaps
     ngrid = length(grid)
@@ -592,7 +614,7 @@ function remez(numtaps::Integer, bands::Vector, desired::Vector;
     l = 1
 
     # Boolean for "kkk" in C code.
-    full_grid = (bands[1] == 0.0 && bands[end] == 0.5) || (nfcns <= 3)
+    full_grid = (band_defs[1].first[1] == 0.0 && band_defs[end].first[2] == 0.5*Hz) || (nfcns <= 3)
     if !full_grid
         aa    = 2.0/(grid[1]-grid[ngrid])
         bb    = -(grid[1]+grid[ngrid])/(grid[1]-grid[ngrid])

test/remez_fir.jl

@@ -3,7 +3,8 @@ using DSP, Compat, Compat.Test, Compat.DelimitedFiles, FilterTestHelpers
 
 @testset "remez_argument_check1" begin
     # bands not monotonically increasing
-    @test_throws ArgumentError remez(151, [0, 0.25, 0.23, 0.5], [1.0, 0.0])
+    @test_throws ArgumentError remez(151, [0, 0.25, 0.25, 0.5], [1.0, 0.0])
+    @test_throws ArgumentError remez(151, [0.2, 0.1, 0.25, 0.5], [1.0, 0.0])
 end
 
 @testset "remez_argument_check2" begin
@@ -30,21 +31,25 @@ end
 # Length 151 LPF (Low Pass Filter).
 #
 @testset "remez_151_lpf" begin
-    h = remez(151, [0, 0.475, 0.5, 1.0], [1.0, 0.0]; Hz=2.0);
     h_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "remez_151_lpf.txt"),'\t')
+    h = remez(151, [0, 0.475, 0.5, 1.0], [1.0, 0.0]; Hz=2.0);
+    @test h ≈ h_scipy
+    h = remez(151, [(0, 0.475) => 1, (0.5, 1.0) => 0]; Hz=2.0);
     @test h ≈ h_scipy
 end
 
 #
 # Length 152 LPF. Non-default "weight" input.
-# 
+#
 #    from scipy.signal import remez
 #    lpf = remez(152, [0, 0.475, 0.5, 1.0], [1.0, 0.0], weight=[1,2], Hz=2.0)
 #    lpf.tofile('remez_152_lpf.txt', sep='\n')
 #
 @testset "remez_152_lpf" begin
-    h = remez(152, [0, 0.475, 0.5, 1.0], [1.0, 0.0]; weight=[1,2], Hz=2.0);
     h_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "remez_152_lpf.txt"),'\t')
+    h = remez(152, [0, 0.475, 0.5, 1.0], [1.0, 0.0]; weight=[1,2], Hz=2.0);
+    @test h ≈ h_scipy
+    h = remez(152, [(0, 0.475) => (1, 1), (0.5, 1.0) => (0, 2)]; Hz=2.0);
     @test h ≈ h_scipy
 end
 
@@ -56,8 +61,10 @@ end
 #    hpf.tofile('remez_51_hpf.txt', sep='\n')
 #
 @testset "remez_51_hpf" begin
-    h = remez(51, [0, 0.75, 0.8, 1.0], [0.0, 1.0]; Hz=2.0);
     h_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "remez_51_hpf.txt"),'\t')
+    h = remez(51, [0, 0.75, 0.8, 1.0], [0.0, 1.0]; Hz=2.0);
+    @test h ≈ h_scipy
+    h = remez(51, [(0, 0.75) => 0, (0.8, 1.0) => 1]; Hz=2.0);
     @test h ≈ h_scipy
 end
 
@@ -69,27 +76,32 @@ end
 #    bpf.tofile('remez_180_bpf.txt', sep='\n')
 #
 @testset "remez_180_bpf" begin
+h_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "remez_180_bpf.txt"),'\t')
     h = remez(180, [0, 0.375, 0.4, 0.5, 0.525, 1.0], [0.0, 1.0, 0.0]; Hz=2.0, maxiter=30);
-    h_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "remez_180_bpf.txt"),'\t')
+    @test h ≈ h_scipy
+    h = remez(180, [(0, 0.375) => 0, (0.4, 0.5) => 1, (0.525, 1.0) => 0]; Hz=2.0, maxiter=30);
     @test h ≈ h_scipy
 end
 
 @testset "remez_warn_no_converge_after_maxiter_iterations" begin
     @static if VERSION ≥ v"0.7.0-DEV.2988"
         @test_logs (:warn, r"filter is not converged") remez(180, [0, 0.375, 0.4, 0.5, 0.525, 1.0], [0.0, 1.0, 0.0]; Hz=2.0)
+        @test_logs (:warn, r"filter is not converged") remez(180, [(0, 0.375) => 0, (0.4, 0.5) => 1, (0.525, 1.0) => 0]; Hz=2.0)
     else
         warn_check(msg) = contains(msg, "filter is not converged")
         @test_warn warn_check remez(180, [0, 0.375, 0.4, 0.5, 0.525, 1.0], [0.0, 1.0, 0.0]; Hz=2.0)
+        @test_warn warn_check remez(180, [(0, 0.375) => 0, (0.4, 0.5) => 1, (0.525, 1.0) => 0]; Hz=2.0)
     end
 end
 
 @testset "remez_error_no_converge_transition_band_too_wide" begin
     @test_throws ErrorException remez(151, [0, 0.1, 0.4, 0.5], [1.0, 0.0])
+    @test_throws ErrorException remez(151, [(0, 0.1) => 1, (0.4, 0.5) => 0])
 end
 
 #
 #  Odd-symmetric filters - hilbert and differentiators type.
-#  Even length - much better approximation since it is not constrained to 0 at 
+#  Even length - much better approximation since it is not constrained to 0 at
 #  the nyquist frequency
 #
 # Length 20 hilbert
@@ -99,8 +111,10 @@ end
 #    h.tofile('remez_20_hilbert.txt', sep='\n')
 #
 @testset "remez_20_hilbert" begin
-    h = remez(20, [0.1, 0.95],[1]; filter_type=filter_type_hilbert, Hz=2.0);
     h_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "remez_20_hilbert.txt"),'\t')
+    h = remez(20, [0.1, 0.95], [1]; filter_type=filter_type_hilbert, Hz=2.0);
+    @test h ≈ h_scipy
+    h = remez(20, [(0.1, 0.95) => 1]; neg=true, Hz=2.0);
     @test h ≈ h_scipy
 end
 
@@ -112,8 +126,10 @@ end
 #    h.tofile('remez_21_hilbert.txt', sep='\n')
 #
 @testset "remez_21_hilbert" begin
-    h = remez(21, [0.1, 0.95],[1]; filter_type=filter_type_hilbert, Hz=2.0);
     h_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "remez_21_hilbert.txt"),'\t')
+    h = remez(21, [0.1, 0.95], [1]; filter_type=filter_type_hilbert, Hz=2.0);
+    @test h ≈ h_scipy
+    h = remez(21, [(0.1, 0.95) => 1]; neg=true, Hz=2.0);
     @test h ≈ h_scipy
 end
 
@@ -125,8 +141,10 @@ end
 #    h.tofile('remez_200_differentiator.txt', sep='\n')
 #
 @testset "remez_200_differentiator" begin
-    h = remez(200, [0.01, 0.99],[1]; filter_type=filter_type_differentiator, Hz=2.0);
     h_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "remez_200_differentiator.txt"),'\t')
+    h = remez(200, [0.01, 0.99], [1]; filter_type=filter_type_differentiator, Hz=2.0);
+    @test h ≈ h_scipy
+    h = remez(200, [(0.01, 0.99) => (f -> f/2, f -> 1/f)]; neg=true, Hz=2.0);
     @test h ≈ h_scipy
 end
 
@@ -138,7 +156,9 @@ end
 #    h.tofile('remez_201_differentiator.txt', sep='\n')
 #
 @testset "remez_201_differentiator" begin
-    h = remez(201, [0.05, 0.95],[1]; filter_type=filter_type_differentiator, Hz=2.0);
     h_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "remez_201_differentiator.txt"),'\t')
+    h = remez(201, [0.05, 0.95], [1]; filter_type=filter_type_differentiator, Hz=2.0);
+    @test h ≈ h_scipy
+    h = remez(201, [(0.05, 0.95) => (f -> f/2, f -> 1/f)]; neg=true, Hz=2.0);
     @test h ≈ h_scipy
 end