From 0759d66be6ab6690d1e3cd918c37888056101b8d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Thimot=C3=A9=20Dupuch?= <50801309+thimotedupuch@users.noreply.github.com> Date: Sat, 5 Sep 2026 22:38:20 +0200 Subject: [PATCH] Refactor Remez exchange into structured Julia and reduce allocations --- benchmark/README.md | 60 +++ benchmark/remez.jl | 90 ++++ benchmark/remez_compatibility.jl | 36 ++ src/Filters/remez_fir.jl | 777 ++++++++++++++----------------- test/remez_fir.jl | 97 +++- 5 files changed, 634 insertions(+), 426 deletions(-) create mode 100644 benchmark/README.md create mode 100644 benchmark/remez.jl create mode 100644 benchmark/remez_compatibility.jl diff --git a/benchmark/README.md b/benchmark/README.md new file mode 100644 index 00000000..e47cc493 --- /dev/null +++ b/benchmark/README.md @@ -0,0 +1,60 @@ +# Remez rewrite validation + +Run from a checkout containing the original revision: + +```sh +julia --startup-file=no benchmark/remez.jl +julia --startup-file=no benchmark/remez_compatibility.jl +julia --project -e 'using Pkg; Pkg.test(; julia_args=["--check-bounds=yes"])' +``` + +The first two scripts use only standard libraries. Their optional first argument +selects a baseline git revision; the default is +`7c798756cca39251da14bb5d38146feb2cdd1717`. They load the original source and the +working-tree source in separate modules in the same Julia process. They neither +change the checkout nor install benchmark dependencies. + +The timing script checks exact coefficient equality for each case, warms up both +implementations, and reports the median of 31 samples of five calls. Sample order +alternates between old/new and new/old. Allocated bytes include the returned +coefficient vector and are measured after compilation. Compilation time and +process resident memory are not measured. Timing results depend on the CPU and +Julia version; they are not assertions in the package tests. + +The differential script compares coefficients bit for bit, warning messages, and +exception types/messages on 600 seeded random designs. It varies tap count, +symmetry, number of bands, weights, grid density, partial/full frequency coverage, +and iteration limits. Both scripts require the baseline revision in git history; +the ordinary package tests do not. + +## Measured results + +Julia 1.12.7, `znver3`, one Julia thread. Microseconds per call: + +| Design | Original | Rewrite | Original bytes | Rewrite bytes | +| --- | ---: | ---: | ---: | ---: | +| 15-tap low-pass | 13.76 | 11.45 | 5,152 | 3,776 | +| 151-tap low-pass | 863.71 | 620.49 | 45,656 | 34,328 | +| 152-tap weighted low-pass | 1,297.99 | 922.68 | 45,592 | 34,328 | +| 51-tap high-pass | 59.44 | 42.47 | 16,080 | 11,976 | +| 180-tap band-pass | 2,430.06 | 1,771.71 | 53,064 | 39,848 | +| 20-tap Hilbert | 10.71 | 9.20 | 6,304 | 4,736 | +| 21-tap Hilbert | 9.35 | 8.53 | 6,304 | 4,736 | +| 200-tap differentiator | 954.59 | 715.02 | 59,808 | 45,272 | +| 201-tap differentiator | 1,001.03 | 770.01 | 55,712 | 42,200 | +| 63-tap partial-band design | 78.45 | 56.03 | 15,432 | 11,400 | +| 63-tap low-pass, density 64 | 270.58 | 186.17 | 65,512 | 49,032 | +| 511-tap low-pass | 5,852.15 | 4,957.51 | 151,888 | 114,768 | +| 201-tap inverse sinc | 1,240.34 | 929.87 | 59,024 | 44,808 | +| Three-argument weighted low-pass | 1,311.01 | 934.69 | 45,944 | 34,680 | +| Three-argument Hilbert, 20 taps | 11.42 | 9.76 | 6,672 | 5,104 | +| Three-argument Hilbert, 21 taps | 9.75 | 8.51 | 6,672 | 5,104 | +| Three-argument differentiator, 200 taps | 964.16 | 729.10 | 60,400 | 45,816 | +| Three-argument differentiator, 201 taps | 1,041.57 | 782.92 | 56,304 | 42,744 | + +These cases take approximately 9–31% less time and allocate 24–27% fewer bytes. +The rewrite fills the frequency grid directly and reuses interpolation workspace +for the partial-band basis conversion. Its barycentric evaluator groups four +independent divisions while preserving accumulation order with fused operations. +The exchange step keeps the original search order, endpoint margin, and stopping +rules, including the returned approximation when `maxiter` is reached. diff --git a/benchmark/remez.jl b/benchmark/remez.jl new file mode 100644 index 00000000..bfd0a7c9 --- /dev/null +++ b/benchmark/remez.jl @@ -0,0 +1,90 @@ +# julia --startup-file=no benchmark/remez.jl [baseline git revision] +# No non-stdlib dependencies; both implementations run in the same process. +using Printf + +module BaselineRemez + revision = isempty(ARGS) ? "7c798756cca39251da14bb5d38146feb2cdd1717" : ARGS[1] + source = read(`git -C $(@__DIR__) show $revision:src/Filters/remez_fir.jl`, String) + include_string(@__MODULE__, source, "baseline_remez.jl") +end + +module CurrentRemez + include("../src/Filters/remez_fir.jl") +end + +const REMEZ_CASES = [ + ("small LP", 15, [(0., .2) => 1., (.3, .5) => 0.], (;)), + ("LP odd", 151, [(0., .475) => 1., (.5, 1.) => 0.], (; Hz=2.)), + ("LP even weighted", 152, [(0., .475) => (1., 1.), (.5, 1.) => (0., 2.)], (; Hz=2.)), + ("HP", 51, [(0., .75) => 0., (.8, 1.) => 1.], (; Hz=2.)), + ("BP", 180, [(0., .375) => 0., (.4, .5) => 1., (.525, 1.) => 0.], (; Hz=2., maxiter=30)), + ("Hilbert even", 20, [(.1, .95) => 1.], (; Hz=2., neg=true)), + ("Hilbert odd", 21, [(.1, .95) => 1.], (; Hz=2., neg=true)), + ("diff even", 200, [(.01, .99) => (f -> f / 2, f -> 1 / f)], (; Hz=2., neg=true)), + ("diff odd", 201, [(.05, .95) => (f -> f / 2, f -> 1 / f)], (; Hz=2., neg=true)), + ("partial band", 63, [(.05, .2) => 1., (.25, .45) => 0.], (;)), + ("dense grid", 63, [(0., .2) => 1., (.23, .5) => 0.], (; grid_density=64)), + ("long LP", 511, [(0., .2) => 1., (.205, .5) => 0.], (; maxiter=100)), + ("inverse sinc", 201, [(0., 2880.) => (f -> inv(sinc(f / 4800)), 1.), + (10000., 153600.) => (0., 100.)], (; Hz=307200.)), +] + +function elapsed_batch(f, args, kwargs, repetitions) + start = time_ns() + for _ in 1:repetitions + f(args...; kwargs...) + end + return (time_ns() - start) / repetitions +end + +function compare_case(name, args, kwargs; newkwargs=kwargs, samples=31, repetitions=5) + old, new = BaselineRemez.remez, CurrentRemez.remez + reference = old(args...; kwargs...) + @assert isequal(reference, new(args...; newkwargs...)) name + # Compile both the filters and the timing harness before collecting samples. + elapsed_batch(old, args, kwargs, 2) + elapsed_batch(new, args, newkwargs, 2) + oldbytes = @allocated old(args...; kwargs...) + newbytes = @allocated new(args...; newkwargs...) + oldtimes, newtimes = Float64[], Float64[] + GC.gc() + for i in 1:samples + # Alternate order to reduce systematic thermal / frequency bias. + if isodd(i) + push!(oldtimes, elapsed_batch(old, args, kwargs, repetitions)) + push!(newtimes, elapsed_batch(new, args, newkwargs, repetitions)) + else + push!(newtimes, elapsed_batch(new, args, newkwargs, repetitions)) + push!(oldtimes, elapsed_batch(old, args, kwargs, repetitions)) + end + end + oldtime = sort!(oldtimes)[cld(samples, 2)] / 1000 + newtime = sort!(newtimes)[cld(samples, 2)] / 1000 + @printf("%-24s %10.2f %10.2f %7.3f %10d %10d\n", + name, oldtime, newtime, newtime / oldtime, oldbytes, newbytes) + return (; oldtime, newtime, oldbytes, newbytes) +end + +function main() + println("Julia ", VERSION, " / ", Sys.CPU_NAME, " / threads=", Threads.nthreads()) + println("Times: median microseconds per call; bytes: allocations per call, after warmup.") + @printf("%-24s %10s %10s %7s %10s %10s\n", + "case", "old μs", "new μs", "new/old", "old bytes", "new bytes") + for (name, taps, bands, kwargs) in REMEZ_CASES + compare_case(name, (taps, bands), kwargs) + end + compare_case("three-argument weighted", (152, [0., .475, .5, 1.], [1., 0.]), + (; weight=[1., 2.], Hz=2.)) + for n in (20, 21) + compare_case("three-arg Hilbert $n", (n, [.1, .95], [1.]), + (; Hz=2., filter_type=BaselineRemez.filter_type_hilbert); + newkwargs=(; Hz=2., filter_type=CurrentRemez.filter_type_hilbert)) + end + for (n, edges) in ((200, [.01, .99]), (201, [.05, .95])) + compare_case("three-arg diff $n", (n, edges, [1.]), + (; Hz=2., filter_type=BaselineRemez.filter_type_differentiator); + newkwargs=(; Hz=2., filter_type=CurrentRemez.filter_type_differentiator)) + end +end + +abspath(PROGRAM_FILE) == (@__FILE__) && main() diff --git a/benchmark/remez_compatibility.jl b/benchmark/remez_compatibility.jl new file mode 100644 index 00000000..74ae137e --- /dev/null +++ b/benchmark/remez_compatibility.jl @@ -0,0 +1,36 @@ +# julia --startup-file=no benchmark/remez_compatibility.jl [baseline git revision] +# Keep this differential check outside the package tests: it requires git history. +using Test, Random, Logging +include("remez.jl") + +function outcome(f, taps, bands, kwargs) + logger = Test.TestLogger() + result = with_logger(logger) do + try + f(taps, bands; kwargs...) + catch e + (typeof(e), sprint(showerror, e)) + end + end + return result, [(entry.level, entry.message) for entry in logger.logs] +end + +@testset "Remez differential compatibility" begin + rng = MersenneTwister(9138) + for trial in 1:600 + taps = rand(rng, 6:201) + nbands = rand(rng, 1:4) + edges = sort!(rand(rng, 2nbands) ./ 2) + if isodd(trial) + edges[1], edges[end] = 0., .5 + end + bands = [(edges[2i-1], edges[2i]) => (rand(rng), 10.0^rand(rng, -1:1)) + for i in 1:nbands] + kwargs = (; neg=rand(rng, Bool), grid_density=rand(rng, (4, 16, 32)), + maxiter=rand(rng, (1, 2, 10, 25, 40))) + @testset "case $trial: taps=$taps, $kwargs" begin + @test isequal(outcome(BaselineRemez.remez, taps, bands, kwargs), + outcome(CurrentRemez.remez, taps, bands, kwargs)) + end + end +end diff --git a/src/Filters/remez_fir.jl b/src/Filters/remez_fir.jl index e1908f75..d611613f 100644 --- a/src/Filters/remez_fir.jl +++ b/src/Filters/remez_fir.jl @@ -91,161 +91,135 @@ C CODE BANNER @enum RemezFilterType filter_type_bandpass filter_type_differentiator filter_type_hilbert - -""" - lagrange_interp(k::Integer, n::Integer, m::Integer, x::AbstractVector) - -CALCULATE THE LAGRANGE INTERPOLATION COEFFICIENTS -""" -function lagrange_interp(k::Integer, n::Integer, m::Integer, x::AbstractVector) - retval = 1.0 - q = x[k] - for l = 1:m, j = l:m:n - if j != k - retval *= 2.0 * (q - x[j]) - end +# The strided product order limits roundoff in high-order interpolation. +function _remez_barycentric_weight(k, count, stride, nodes) + product = 1.0 + node = nodes[k] + @inbounds for offset in 1:stride, j in offset:stride:count + j == k && continue + product *= 2.0 * (node - nodes[j]) end - 1.0 / retval + return inv(product) end +_remez_callable(value::Real) = Returns(value) +_remez_callable(value) = value -""" - build_grid(numtaps, band_defs, Hz, grid_density, neg) +function _remez_band_edges(band, Hz, lower, upper) + low = Float64(clamp(band.first[1] / Hz, lower, upper)) + high = Float64(clamp(band.first[2] / Hz, lower, upper)) + return low, high +end -Returns `grid`, `des`, and `wt` arrays -""" -function build_grid(nfilt, band_defs, Hz, grid_density, neg::Bool) - nodd = isodd(nfilt) - nfcns = nfilt ÷ 2 - nfcns += nodd & !neg - - # - # SET UP THE DENSE GRID. THE NUMBER OF POINTS IN THE GRID - # IS (FILTER LENGTH + 1)*GRID DENSITY/2 - # - 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[1] - weight = b.second[2] - else - desired = b.second - weight = 1.0 - end - if isa(desired, Real) - desired = Returns(desired) - end - if isa(weight, Real) - weight = Returns(weight) - end - return Pair{Tuple{Float64,Float64},Tuple{Any,Any}}((fl, fu), (desired, weight)) +function _remez_grid(numtaps, band_defs, Hz, density, neg) + odd = isodd(numtaps) + ncosines = numtaps ÷ 2 + (odd && !neg) + step = 0.5 / (density * ncosines) + lower = neg ? step : 0.0 + upper = neg == odd ? 0.5 - step : 0.5 + + # Count first, then fill directly; no normalized band or frequency arrays. + ngrid = 0 + for band in band_defs + low, high = _remez_band_edges(band, Hz, lower, upper) + ngrid += max(length(low:step:high), 1) end - normalized_band_defs = normalize_banddef_entry.(band_defs) - - # work around JuliaLang/julia#15276 - ngrid = let delf=delf - sum(max(length(band_def.first[1]:delf:band_def.first[2]), 1) for band_def in normalized_band_defs)#::Int + grid = Vector{Float64}(undef, ngrid) + desired = similar(grid) + weight = similar(grid) + + index = 1 + for band in band_defs + low, high = _remez_band_edges(band, Hz, lower, upper) + response, importance = band.second isa Tuple{Any,Any} ? + band.second : (band.second, 1.0) + index = _remez_grid_band!(grid, desired, weight, index, low:step:high, + high, neg, odd, Hz, + _remez_callable(response), _remez_callable(importance)) end + @assert index == ngrid + 1 + return grid, desired, weight +end - grid = zeros(Float64, ngrid) # the array of frequencies, between 0 and 0.5 - des = zeros(Float64, ngrid) # the desired function on the grid - wt = zeros(Float64, ngrid) # array of weights - - j = 1 - - # - # CALCULATE THE DESIRED MAGNITUDE RESPONSE AND THE WEIGHT - # FUNCTION ON THE GRID - # - # and - # - # SET UP A NEW APPROXIMATION PROBLEM WHICH IS EQUIVALENT - # TO THE ORIGINAL PROBLEM - # - for band_def in normalized_band_defs - flow, fup = band_def.first - # 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) +# Specialize on each band's callables, including heterogeneous band definitions. +function _remez_grid_band!(grid, desired, weight, index, frequencies, high, neg, odd, Hz, + response::D, importance::W) where {D,W} + count = max(length(frequencies), 1) + for i in 1:count + # The final point is exactly the upper edge, including zero-width bands. + frequency = i == count ? high : frequencies[i] + change = neg ? sinpi(odd ? 2frequency : frequency) : + (odd ? 1.0 : cospi(frequency)) + grid[index] = cospi(2frequency) + desired[index] = response(frequency * Hz) / change + weight[index] = importance(frequency * Hz) * change + index += 1 end - @assert ngrid == j - 1 - - return grid, des, wt + return index end -function _buildgrid!(grid, des, wt, j, fs, neg, nodd, Hz, des_wt) - desired, weight = des_wt - for f in fs - change = neg ? sinpi(nodd ? 2f : f) : (nodd ? 1.0 : cospi(f)) - grid[j] = cospi(2f) - des[j] = desired(f * Hz) / change - wt[j] = weight(f * Hz) * change - j += 1 +# Barycentric evaluation away from the interpolation nodes. Callers use the +# stored values at nodes, so the hot loop needs no equality check or allocation. +@inline function _remez_eval(frequency, nodes, values, weights) + denominator = 0.0 + numerator = 0.0 + # Independent divisions can run in parallel. Keep both accumulations in + # their original order: reassociating these sums changes Remez convergence. + j = 1 + @inbounds while j + 3 <= length(weights) + c1 = weights[j] / (frequency - nodes[j]) + c2 = weights[j+1] / (frequency - nodes[j+1]) + c3 = weights[j+2] / (frequency - nodes[j+2]) + c4 = weights[j+3] / (frequency - nodes[j+3]) + denominator += c1 + numerator = fma(c1, values[j], numerator) + denominator += c2 + numerator = fma(c2, values[j+1], numerator) + denominator += c3 + numerator = fma(c3, values[j+2], numerator) + denominator += c4 + numerator = fma(c4, values[j+3], numerator) + j += 4 + end + @inbounds for j in j:length(weights) + coefficient = weights[j] / (frequency - nodes[j]) + denominator += coefficient + numerator = muladd(coefficient, values[j], numerator) end - return j + return numerator / denominator end -""" - function freq_eval(xf, x::AbstractVector, y::AbstractVector, ad::AbstractVector) ------------------------------------------------------------------------ - FUNCTION: freq_eval (gee) - FUNCTION TO EVALUATE THE FREQUENCY RESPONSE USING THE - LAGRANGE INTERPOLATION FORMULA IN THE BARYCENTRIC FORM - ------------------------------------------------------------------------ -""" -function freq_eval(xf, x::AbstractVector, y::AbstractVector, ad::AbstractVector) - d = 0.0 - p = 0.0 - - for j in eachindex(ad) - c = ad[j] / (xf - x[j]) - d += c - p = muladd(c, y[j], p) +# Fit the alternating weighted error on the current extremal set. +function _remez_interpolate!(nodes, values, weights, extrema, grid, desired, importance) + count = length(values) + stride = (count - 2) ÷ 15 + 1 + for j in 1:count + nodes[j] = grid[extrema[j]] + end + for j in 1:count + weights[j] = _remez_barycentric_weight(j, count, stride, nodes) end - p/d -end - -function initialize_y(dev::Float64, nz::Integer, iext::AbstractArray, des::AbstractArray, wt::AbstractArray, y::AbstractArray) - nu = 1 - if dev > 0.0 - nu = -1 + numerator = denominator = 0.0 + sign = 1 + for j in 1:count + index = extrema[j] + numerator = muladd(weights[j], desired[index], numerator) + denominator = muladd(sign, weights[j] / importance[index], denominator) + sign = -sign end - dev = -nu * dev - k = nu - for j = 1 : nz - l = iext[j] - y[j] = des[l] + k * dev / wt[l] - k = -k + deviation = numerator / denominator + first_sign = deviation > 0.0 ? -1 : 1 + ripple = -first_sign * deviation + sign = first_sign + for j in 1:count + index = extrema[j] + values[j] = desired[index] + sign * ripple / importance[index] + sign = -sign end - nu, dev + return first_sign, ripple end -#========= -Banner from C code - ------------------------------------------------------------------------ - SUBROUTINE: remez - THIS SUBROUTINE IMPLEMENTS THE REMEZ EXCHANGE ALGORITHM - FOR THE WEIGHTED CHEBYSHEV APPROXIMATION OF A CONTINUOUS - FUNCTION WITH A SUM OF COSINES. INPUTS TO THE SUBROUTINE - ARE A DENSE GRID WHICH REPLACES THE FREQUENCY AXIS, THE - DESIRED FUNCTION ON THIS GRID, THE WEIGHT FUNCTION ON THE - GRID, THE NUMBER OF COSINES, AND AN INITIAL GUESS OF THE - EXTREMAL FREQUENCIES. THE PROGRAM MINIMIZES THE CHEBYSHEV - ERROR BY DETERMINING THE BSMINEST LOCATION OF THE EXTREMAL - FREQUENCIES (POINTS OF MAXIMUM ERROR) AND THEN CALCULATES - THE COEFFICIENTS OF THE BEST APPROXIMATION. ------------------------------------------------------------------------ -=========# """ remez(numtaps::Integer, band_defs; Hz::Real=1.0, @@ -396,366 +370,320 @@ function remez(numtaps::Integer, band_defs; (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, band_defs, Hz, grid_density, neg) - - nfilt = numtaps - ngrid = length(grid) - - nodd = isodd(nfilt) # boolean: "nodd" means filter length is odd - nfcns = numtaps ÷ 2 # integer divide - if nodd && !neg - nfcns = nfcns + 1 - end - - nz = nfcns+1 - nzz = nfcns+2 - iext = zeros(Int64, nzz) # indices of extremals - x = zeros(Float64, nzz) - y = zeros(Float64, nz) + grid, desired, weight = _remez_grid(numtaps, band_defs, Hz, grid_density, neg) + full_grid = band_defs[1].first[1] == 0.0 && band_defs[end].first[2] == 0.5Hz + return _remez_solve(numtaps, neg, maxiter, full_grid, grid, desired, weight) +end - for j = 1:nz - iext[j] = (j-1)*(ngrid-1) ÷ nfcns + 1 +function _remez_solve(numtaps, neg, maxiter, full_grid, grid, desired, weight) + ncosines = numtaps ÷ 2 + (isodd(numtaps) && !neg) + count = ncosines + 1 + extrema = Vector{Int}(undef, count + 1) + nodes = zeros(count + 1) # Extra node is a sentinel during coefficient recovery. + values = zeros(count) + weights = zeros(count) + for j in 1:count + extrema[j] = (j - 1) * (length(grid) - 1) ÷ ncosines + 1 end + extrema[end] = length(grid) + 1 - dev = 0.0 # deviation from the desired function, - # that is, the amount of "ripple" on the extremal set - devl = -1.0 # deviation on last iteration - niter = 0 - ad = zeros(Float64, nz) - - jet = ((nfcns-1) ÷ 15) + 1 - + previous_ripple = -1.0 + iteration = 0 while true - - # - # Start next iteration - # - # @label L100 - iext[nzz] = ngrid + 1 - niter += 1 - if niter > maxiter + iteration += 1 + if iteration > maxiter @warn("remez() iteration count exceeds maxiter = $maxiter, filter is not converged; try increasing maxiter") - # the filter is returned in its current, unconverged state. break end - for j = 1:nz - x[j] = grid[iext[j]] - end - - for j = 1 : nz - ad[j] = lagrange_interp(j, nz, jet, x) + first_sign, ripple = _remez_interpolate!(nodes, values, weights, extrema, + grid, desired, weight) + if ripple <= previous_ripple + throw(ErrorException("remez() - failure to converge at iteration $iteration, try reducing transition band width")) end + previous_ripple = ripple + _remez_exchange!(extrema, grid, desired, weight, nodes, values, weights, + first_sign, ripple) || break + end - dnum = 0.0 - dden = 0.0 - k = 1 - for j = 1 : nz - l = iext[j] - dnum = muladd(ad[j], des[l], dnum) - dden = muladd(k, ad[j] / wt[l], dden) - k = -k - end - dev = dnum / dden + return _remez_coefficients(numtaps, neg, full_grid, grid, nodes, values, weights) +end - fill!(y, 0.0) - nu, dev = initialize_y(dev, nz, iext, des, wt, y) +# Search a monotone run of signed error, stopping before the exclusive bound. +# The direction is a value parameter so the two hot loops specialize separately. +@inline function _remez_peak(error, index, bound, amplitude, ::Val{step}) where {step} + while index + step != bound + candidate = error(index + step) + candidate <= amplitude && break + index += step + amplitude = candidate + end + return index, amplitude +end - if dev <= devl - # finished - throw(ErrorException("remez() - failure to converge at iteration $niter, try reducing transition band width")) +# Exchange extrema in place. Searching each interval in the same order preserves +# the original algorithm's tie breaking and its behavior at the iteration limit. +function _remez_exchange!(extrema, grid, desired, weight, nodes, values, baryweights, + first_sign, ripple) + count = length(values) + gridlength = length(grid) + oldfirst, oldlast = extrema[1], extrema[count] + lower = 0 + changed = false + first_amplitude = last_amplitude = ripple + sign = first_sign + + for j in 1:count + current, upper = extrema[j], extrema[j + 1] + error = let sign = sign + i -> sign * ((_remez_eval(grid[i], nodes, values, baryweights) - + desired[i]) * weight[i]) end - devl = dev - - # - # SEARCH FOR THE EXTREMAL FREQUENCIES OF THE BEST APPROXIMATION - # - - # Between here and L370, the extremal index set is updated in a loop - # roughly over the index "j" - although the logic is complicated as - # the extremal set may grow or shrink in an iteration. - # j - the index of the current extremal being updated - # nz - the number of cosines in the approximation (including the constant term). - # nz = nfcns + 1 where nfcns = nfilt / 2, and - # nfilt is the filter length or number of taps. - # For example, for a length 15 filter, nfcns = 7 and nz = 8. - # jchgne - the number of extremal indices that changed this iteration - jchnge = 0 - k1 = iext[1] - knz = iext[nz] - klow = 0 - nut = -nu - j = 1 - - local comp - - @label L200 - j == nzz && (ynz = comp) # equivalent to "if (j == nzz) ynz = comp; end" - j >= nzz && @goto L300 - kup = iext[j+1] - l = iext[j]+1 - nut = -nut - j == 2 && (y1 = comp) - comp = dev - l >= kup && @goto L220 - err = (freq_eval(grid[l], x, y, ad) - des[l]) * wt[l] - nut*err <= comp && @goto L220 - comp = nut * err - @label L210 - l += 1; l >= kup && @goto L215 - err = (freq_eval(grid[l], x, y, ad) - des[l]) * wt[l] - nut*err <= comp && @goto L215 - comp = nut * err - @goto L210 - - @label L215 - iext[j] = l - 1; j += 1 - klow = l - 1 - jchnge += 1 - @goto L200 - - @label L220 - l -= 1 - @label L225 - l -= 1; l <= klow && @goto L250 - err = (freq_eval(grid[l], x, y, ad) - des[l]) * wt[l] - nut*err > comp && @goto L230 - jchnge <= 0 && @goto L225 - @goto L260 - - @label L230 - comp = nut * err - @label L235 - l -= 1; l <= klow && @goto L240 - err = (freq_eval(grid[l], x, y, ad) - des[l]) * wt[l] - nut*err <= comp && @goto L240 - comp = nut * err - @goto L235 - @label L240 - klow = iext[j] - iext[j] = l+1 - j += 1 - jchnge += 1 - @goto L200 - - @label L250 - l = iext[j]+1 - jchnge > 0 && @goto L215 - - @label L255 - l += 1; l >= kup && @goto L260 - err = (freq_eval(grid[l], x, y, ad) - des[l]) * wt[l] - nut*err <= comp && @goto L255 - comp = nut * err - - @goto L210 - @label L260 - klow = iext[j]; j += 1 - @goto L200 - - @label L300 - j > nzz && @goto L320 - k1 > iext[1] && (k1 = iext[1]) - knz < iext[nz] && (knz = iext[nz]) - nut1 = nut - nut = -nu - l = 0 - kup = k1 - comp = ynz*(1.00001) - luck = 1 - @label L310 - l += 1; l >= kup && @goto L315 - err = (freq_eval(grid[l], x, y, ad) - des[l]) * wt[l] - nut*err <= comp && @goto L310 - comp = nut * err - j = nzz - @goto L210 - - @label L315 - luck = 6 - @goto L325 - - @label L320 - luck > 9 && @goto L350 - comp > y1 && (y1 = comp) - k1 = iext[nzz] - @label L325 - l = ngrid+1 - klow = knz - nut = -nut1 - comp = y1*(1.00001) - @label L330 - l -= 1; l <= klow && @goto L340 - err = (freq_eval(grid[l], x, y, ad) - des[l]) * wt[l] - nut*err <= comp && @goto L330 - j = nzz - comp = nut * err - luck = luck + 10 - @goto L235 - @label L340 - luck == 6 && @goto L370 - for j = 1 : nfcns - iext[nzz-j] = iext[nz-j] - end - iext[1] = k1 - continue # @goto L100 - @label L350 - for j = 1:nz - iext[j] = iext[j+1] + amplitude = ripple + candidate = current + 1 + + value = candidate < upper ? error(candidate) : amplitude + if value > amplitude + candidate, amplitude = _remez_peak(error, candidate, upper, + value, Val(1)) + lower = candidate + changed = true + else + candidate = current - 1 + while candidate > lower + value = error(candidate) + if value > amplitude + candidate, amplitude = _remez_peak(error, candidate, lower, + value, Val(-1)) + break + end + changed && break + candidate -= 1 + end + + if candidate > lower && amplitude > ripple + # Keep the old extremum as the lower bound when moving left. + lower = current + changed = true + elseif candidate <= lower && changed + candidate = current + lower = current + elseif !changed + candidate = current + 2 + while candidate < upper + value = error(candidate) + value > amplitude && break + candidate += 1 + end + if candidate < upper + candidate, amplitude = _remez_peak(error, candidate, upper, + value, Val(1)) + changed = true + else + candidate = current + end + lower = candidate + else + candidate = current + lower = current + end end - continue # @goto L100 - @label L370 + extrema[j] = candidate + j == 1 && (first_amplitude = amplitude) + last_amplitude = amplitude + sign = -sign + end + # An extra extremum at either end may displace the opposite endpoint. + # The small margin avoids cycling on nearly equal endpoint errors. + left_bound = min(oldfirst, extrema[1]) + right_bound = max(oldlast, extrema[count]) + left_error = i -> -first_sign * ((_remez_eval(grid[i], nodes, values, baryweights) - + desired[i]) * weight[i]) + left = 1 + threshold = last_amplitude * 1.00001 + amplitude = threshold + while left < left_bound + amplitude = left_error(left) + amplitude > threshold && break + left += 1 + end + have_left = left < left_bound + if have_left + left, amplitude = _remez_peak(left_error, left, left_bound, + amplitude, Val(1)) + first_amplitude = max(first_amplitude, amplitude) + end - if jchnge <= 0 # we are done if none of the extremal indices changed - break + right_error = let sign = sign + i -> sign * ((_remez_eval(grid[i], nodes, values, baryweights) - + desired[i]) * weight[i]) + end + right = gridlength + threshold = first_amplitude * 1.00001 + while right > right_bound + amplitude = right_error(right) + amplitude > threshold && break + right -= 1 + end + if right > right_bound + right, _ = _remez_peak(right_error, right, right_bound, + amplitude, Val(-1)) + for j in 1:count-1 + extrema[j] = extrema[j + 1] + end + extrema[count] = right + return true + elseif have_left + for j in count:-1:2 + extrema[j] = extrema[j - 1] end - end # while + extrema[1] = left + return true + end + return changed +end - # - # CALCULATION OF THE COEFFICIENTS OF THE BEST APPROXIMATION - # USING THE INVERSE DISCRETE FOURIER TRANSFORM - # +function _remez_coefficients(numtaps, neg, full_grid, grid, nodes, values, weights) + ncosines = length(values) - 1 + count, workspace_length = ncosines + 1, ncosines + 2 + gridlength = length(grid) + # Sample the interpolant, then recover its cosine coefficients by an IDFT. - a = zeros(Float64, nfcns) # frequency response on evenly spaced grid - p = zeros(Float64, nfcns) - q = zeros(Float64, nfcns-2) - alpha = zeros(Float64, nzz) # return vector + samples = Vector{Float64}(undef, ncosines) # frequency response on evenly spaced grid + coefficients = zeros(Float64, workspace_length) # return vector - fsh = 1.0e-06 - x[nzz] = -2.0 - delf = 1 / (2*nfcns - 1) - l = 1 + node_tolerance = 1.0e-06 + nodes[workspace_length] = -2.0 + step = 1 / (2*ncosines - 1) + node_index = 1 - # Boolean for "kkk" in C code. - full_grid = (band_defs[1].first[1] == 0.0 && band_defs[end].first[2] == 0.5*Hz) || (nfcns <= 3) + # Partial grids use an affine change of the Chebyshev variable. + full_grid |= ncosines <= 3 if !full_grid - aa = 2.0/(grid[1]-grid[ngrid]) - bb = -(grid[1]+grid[ngrid])/(grid[1]-grid[ngrid]) + scale = 2.0/(grid[1]-grid[gridlength]) + shift = -(grid[1]+grid[gridlength])/(grid[1]-grid[gridlength]) end - # Fill in "a" array with the frequency response on an evenly - # spaced set of frequencies. Care is taken to use "y[l]" - - # an already computed response on one of the extremals "l" - - # if the extremal is equal to the frequency ft. If no y[l] - # matches, a[j] is computed using freq_eval. - for j = 1 : nfcns - xt = cospi(2 * (j - 1) * delf) + # Reuse values near interpolation nodes to avoid division by zero. + for j = 1:ncosines + frequency = cospi(2 * (j - 1) * step) if !full_grid - xt = (xt-bb)/aa + frequency = (frequency-shift)/scale end - if (l > 1) - l = l-1 + if (node_index > 1) + node_index = node_index-1 end - while x[l]-xt >= fsh - l += 1 + while nodes[node_index]-frequency >= node_tolerance + node_index += 1 end - if xt-x[l] < fsh - a[j] = y[l] + if frequency-nodes[node_index] < node_tolerance + samples[j] = values[node_index] else - a[j] = freq_eval(xt, x, y, ad) + samples[j] = _remez_eval(frequency, nodes, values, weights) end end - nm1 = nfcns - 1 # nm1 => "nfcns minus 1" + degree = ncosines - 1 - for j = 1 : nfcns - dtemp = 0.0 - for k = 1 : nm1 - dtemp = muladd(a[k+1], cospi(2 * (j-1) * delf * k), dtemp) + for j in 1:ncosines + total = 0.0 + for k in 1:degree + total = muladd(samples[k+1], cospi(2 * (j-1) * step * k), total) end - alpha[j] = 2dtemp + a[1] + coefficients[j] = 2total + samples[1] end - for j = 2 : nfcns - alpha[j] *= 2.0 * delf + for j in 2:ncosines + coefficients[j] *= 2.0 * step end - alpha[1] *= delf + coefficients[1] *= step if !full_grid - p[1] = muladd(2alpha[nfcns], bb, alpha[nm1]) - p[2] = 2.0*aa*alpha[nfcns] - q[1] = alpha[nfcns-2]-alpha[nfcns] - for j = 2 : nm1 - if j >= nm1 - aa *= 0.5 - bb *= 0.5 + # Interpolation is finished; reuse its workspace for basis conversion. + p = weights + q = values + p[1] = muladd(2coefficients[ncosines], shift, coefficients[degree]) + p[2] = 2.0*scale*coefficients[ncosines] + q[1] = coefficients[ncosines-2]-coefficients[ncosines] + for j in 2:degree + if j >= degree + scale *= 0.5 + shift *= 0.5 end p[j+1] = 0.0 - for k = 1 : j - a[k] = p[k] - p[k] = 2.0 * bb * a[k] + for k in 1:j + samples[k] = p[k] + p[k] = 2.0 * shift * samples[k] end - p[2] = muladd(a[1], 2aa, p[2]) - for k = 1 : j-1 - p[k] += muladd(aa, a[k+1], q[k]) + p[2] = muladd(samples[1], 2scale, p[2]) + for k in 1:j-1 + p[k] += muladd(scale, samples[k+1], q[k]) end - for k = 3 : j+1 - p[k] = muladd(aa, a[k-1], p[k]) + for k in 3:j+1 + p[k] = muladd(scale, samples[k-1], p[k]) end - if j != nm1 - for k = 1 : j - q[k] = -a[k] + if j != degree + for k in 1:j + q[k] = -samples[k] end - q[1] += alpha[nfcns - 1 - j] + q[1] += coefficients[ncosines - 1 - j] end end - for j = 1 : nfcns - alpha[j] = p[j] + for j in 1:ncosines + coefficients[j] = p[j] end end - if nfcns <= 3 - alpha[nfcns+1] = alpha[nfcns+2] = 0.0 - end + return _remez_impulse_response(numtaps, neg, coefficients, ncosines) +end - # - # CALCULATE THE IMPULSE RESPONSE. - # - h = zeros(Float64, nfilt) +# Undo the linear-phase transformation for FIR types I, II, III and IV. +function _remez_impulse_response(numtaps, neg, coefficients, ncosines) + odd = isodd(numtaps) + count = ncosines + 1 + degree = ncosines - 1 + h = Vector{Float64}(undef, numtaps) if !neg - if nodd - for j = 1 : nm1 - h[j] = 0.5 * alpha[nz-j] + if odd + for j in 1:degree + h[j] = 0.5 * coefficients[count-j] end - h[nfcns] = alpha[1] + h[ncosines] = coefficients[1] else - h[1] = 0.25 * alpha[nfcns] - for j = 2 : nm1 - h[j] = 0.25 * (alpha[nz-j] + alpha[nfcns+2-j]) + h[1] = 0.25 * coefficients[ncosines] + for j in 2:degree + h[j] = 0.25 * (coefficients[count-j] + coefficients[ncosines+2-j]) end - h[nfcns] = muladd(0.5, alpha[1], 0.25 * alpha[2]) + h[ncosines] = muladd(0.5, coefficients[1], 0.25 * coefficients[2]) end else - if nodd - h[1] = 0.25 * alpha[nfcns] - h[2] = 0.25 * alpha[nm1] - for j = 1 : nm1 - h[j] = 0.25 * (alpha[nz-j] - alpha[nfcns+3-j]) + if odd + for j in 1:degree + h[j] = 0.25 * (coefficients[count-j] - coefficients[ncosines+3-j]) end - h[nfcns] = muladd(0.5, alpha[1], -0.25 * alpha[3]) - h[nz] = 0.0 + h[ncosines] = muladd(0.5, coefficients[1], -0.25 * coefficients[3]) else - h[1] = 0.25 * alpha[nfcns] - for j = 2 : nm1 - h[j] = 0.25 * (alpha[nz-j] - alpha[nfcns+2-j]) + h[1] = 0.25 * coefficients[ncosines] + for j in 2:degree + h[j] = 0.25 * (coefficients[count-j] - coefficients[ncosines+2-j]) end - h[nfcns] = muladd(0.5, alpha[1], -0.25 * alpha[2]) + h[ncosines] = muladd(0.5, coefficients[1], -0.25 * coefficients[2]) end end - for j = 1 : nfcns - k = nfilt + 1 - j + for j in 1:ncosines + k = numtaps + 1 - j if !neg h[k] = h[j] else h[k] = -h[j] end end - if neg && nodd - h[nz] = 0.0 + if neg && odd + h[count] = 0.0 end return h @@ -853,4 +781,3 @@ function remez(numtaps::Integer, bands::Vector, desired::Vector; neg = filter_type in (filter_type_hilbert, filter_type_differentiator) return remez(numtaps, band_defs; Hz, neg, kwargs...) end - diff --git a/test/remez_fir.jl b/test/remez_fir.jl index d93a4a57..3c52599b 100644 --- a/test/remez_fir.jl +++ b/test/remez_fir.jl @@ -7,6 +7,102 @@ using DSP, Test, DelimitedFiles, FilterTestHelpers @test_throws ArgumentError remez(151, [0.2, 0.1, 0.25, 0.5], [1.0, 0.0]) end +@testset "remez API equivalence and symmetry" begin + for n in (20, 21), density in (8, 32) + bands = [.05, .2, .3, .45] + desired, weight = [1., 0.], [2., 3.] + definitions = [(.05, .2) => (1., 2.), (.3, .45) => (0., 3.)] + h = remez(n, bands, desired; weight, grid_density=density) + @test h == remez(n, definitions; grid_density=density) + @test h == reverse(h) + @test eltype(h) === Float64 + @test length(h) == n + # Changing frequency units must leave the design unchanged. + @test h ≈ remez(n, 1000bands, desired; weight, Hz=1000, grid_density=density) + if iseven(n) + @test abs(sum(h .* (-1) .^ (0:n-1))) < 1e-14 + end + + h = remez(n, [.05, .45], [1.]; filter_type=filter_type_hilbert, + grid_density=density) + @test h == remez(n, [(.05, .45) => 1.]; neg=true, grid_density=density) + @test h == -reverse(h) + @test abs(sum(h)) < 1e-14 + if isodd(n) + @test h[cld(n, 2)] == 0 + @test abs(sum(h .* (-1) .^ (0:n-1))) < 1e-14 + end + + # Both branches of the differentiator's relative weighting rule. + h = remez(n, bands, desired; weight, filter_type=filter_type_differentiator, + grid_density=density) + functional = [(.05, .2) => (f -> f, f -> 2 / f), + (.3, .45) => (f -> 0f, f -> 3.)] + @test h == remez(n, functional; neg=true, grid_density=density) + @test bands == [.05, .2, .3, .45] + @test desired == [1., 0.] + @test weight == [2., 3.] + end +end + +@testset "remez heterogeneous responses" begin + definitions = ((.02, .2) => (f -> 1 + f, f -> 1 + 2f), (.3, .48) => (0., 3.)) + h = remez(32, definitions) + @test h == remez(32, Any[definitions...]) + @test h == remez(32, [(.02, .2) => (f -> 1 + f, f -> 1 + 2f), + (.3, .48) => (Returns(0.), Returns(3.))]) + @test all(isfinite, h) +end + +@testset "remez exact grid endpoints" begin + # For 15 taps and density 4 the step is 1/64. The upper edge replaces + # the last range point; a band narrower than a step still contributes one. + seen = Float64[] + response = f -> (push!(seen, f); 1.) + bands = [(0., .1) => response, (.2, .201) => response, (.3, .3) => response] + grid, desired, weight = DSP.Filters._remez_grid(15, bands, 1., 4, false) + expected = [0., 1/64, 2/64, 3/64, 4/64, 5/64, .1, .201, .3] + @test seen == expected + @test grid == cospi.(2expected) + @test desired == ones(length(expected)) + @test weight == ones(length(expected)) +end + +@testset "remez iteration limit coefficients" begin + # Snapshot of the original implementation, not just a warning assertion. + half = [-0.0057324947447139635, -0.012308085906187437, 0.026588205554288386, + 0.026835890017840345, -0.07953783263891252, -0.03818487882449576, + 0.3086821218293382, 0.5424953165331308] + expected = [half; reverse(half[1:end-1])] + h = @test_logs (:warn, r"filter is not converged") remez(15, + [(0., .2) => 1., (.3, .5) => 0.]; maxiter=1) + @test h ≈ expected + h = @test_logs (:warn, r"filter is not converged") remez(15, + [0., .2, .3, .5], [1., 0.]; maxiter=1) + @test h ≈ expected +end + +@testset "remez weighted equiripple response" begin + h = remez(31, [(0., .2) => (1., 1.), (.25, .5) => (0., 3.)]) + errors = Float64[] + for (band, desired, weight) in ((range(0., .2; length=2001), 1., 1.), + (range(.25, .5; length=2501), 0., 3.)) + # Evaluate the zero-phase response directly, independently of the + # barycentric interpolation used by the implementation. + response = [sum(h[k] * cospi(2f * (k - 16)) for k in eachindex(h)) for f in band] + error = weight .* (response .- desired) + for i in eachindex(error) + left = i == firstindex(error) || abs(error[i]) >= abs(error[i-1]) + right = i == lastindex(error) || abs(error[i]) >= abs(error[i+1]) + left && right && push!(errors, error[i]) + end + end + @test length(errors) == 17 + @test all(errors[i] * errors[i+1] < 0 for i in 1:length(errors)-1) + # Allow dense-grid discretization error relative to the continuous extrema. + @test maximum(abs, errors) / minimum(abs, errors) < 1.02 +end + @testset "remez_argument_check2" begin # bands values out of range @test_throws ArgumentError remez(151, [0, 0.23, 0.25, 0.6], [1.0, 0.0]) @@ -182,4 +278,3 @@ end # Test that maximum passband error is less than 1/4 dB. @test maximum(abs.(Pdb[passband_indices] - Gdb[passband_indices])) < 0.25 end -