Buffer array resample and reference_data (#623)

Author: wheeheee

Project.toml

@@ -1,6 +1,6 @@
 name = "DSP"
 uuid = "717857b8-e6f2-59f4-9121-6e50c889abd2"
-version = "0.8.0"
+version = "0.8.1"
 
 [deps]
 Bessels = "0e736298-9ec6-45e8-9647-e4fc86a2fe38"

src/Filters/stream_filt.jl

@@ -1,3 +1,5 @@
+using ..DSP: _zeropad
+
 const PFB{T} = Matrix{T}          # polyphase filter bank
 
 abstract type FIRKernel{T} end
@@ -132,11 +134,11 @@ end
 FIRArbitrary(h::Vector, rate::Real, Nϕ::Integer) = FIRArbitrary(h, convert(Float64, rate), convert(Int, Nϕ))
 
 # FIRFilter - the kernel does the heavy lifting
-mutable struct FIRFilter{Tk<:FIRKernel}
-    kernel::Tk
+mutable struct FIRFilter{Tk<:FIRKernel,T}
+    const kernel::Tk
+    const h::Vector{T}
+    const historyLen::Int
     history::Vector
-    historyLen::Int
-    h::Vector
 end
 
 # Constructor for single-rate, decimating, interpolating, and rational resampling filters
@@ -172,7 +174,7 @@ function FIRFilter(h::Vector, resampleRatio::Union{Integer,Rational} = 1)
 
     history = zeros(historyLen)
 
-    FIRFilter(kernel, history, historyLen, h)
+    FIRFilter(kernel, h, historyLen, history)
 end
 
 # Constructor for arbitrary resampling filter (polyphase interpolator w/ intra-phase linear interpolation)
@@ -193,7 +195,7 @@ function FIRFilter(h::Vector, rate::AbstractFloat, Nϕ::Integer=32)
     kernel     = FIRArbitrary(h, rate, Nϕ)
     historyLen = kernel.tapsPerϕ - 1
     history    = zeros(historyLen)
-    FIRFilter(kernel, history, historyLen, h)
+    FIRFilter(kernel, h, historyLen, history)
 end
 
 # Constructor for a resampling FIR filter, where the user needs only to set the sampling rate
@@ -623,31 +625,34 @@ function filt!(
     return bufIdx
 end
 
-function filt(self::FIRFilter{Tk}, x::AbstractVector{Tx}) where {Th,Tx,Tk<:FIRKernel{Th}}
-    bufLen         = outputlength(self, length(x))
+function filt(self::FIRFilter{Tk}, x::AbstractVector) where Tk<:FIRKernel
+    buffer = allocate_output(self, x)
+    bufLen = length(buffer)
+    samplesWritten = filt!(buffer, self, x)
+    if Tk <: FIRArbitrary
+        samplesWritten == bufLen || resize!(buffer, samplesWritten)
+    else
+        samplesWritten == bufLen || throw(AssertionError("Length of resampled output different from expectation."))
+    end
+    return buffer
+end
+
+function allocate_output(sf::FIRFilter{Tk}, x::AbstractVector{Tx}) where {Th,Tx,Tk<:FIRKernel{Th}}
     # In some cases when `filt(::FIRFilter{FIRArbitrary}, x)` is called
     # with certain values of `x`, `filt!(buffer, ::FIRFilter{FIRArbitrary}, x)`
     # tries to write one sample too many to the buffer and a `BoundsError`
-    # is thrown.  Add one extra sample to catch these exceptional cases.
+    # is thrown. Add one extra sample to catch these exceptional cases.
     #
     # See https://github.com/JuliaDSP/DSP.jl/issues/317
     #
     # FIXME: Remove this if and when the code in
     #        `filt!(buffer, ::FIRFilter{FIRArbitrary}, x)`
     #        is updated to properly account for pathological arbitrary rates.
+    outLen = outputlength(sf, length(x))
     if Tk <: FIRArbitrary
-        bufLen += 1
+        outLen += 1
     end
-    buffer         = Vector{promote_type(Th,Tx)}(undef, bufLen)
-    samplesWritten = filt!(buffer, self, x)
-
-    if Tk <: FIRArbitrary
-        samplesWritten == bufLen || resize!(buffer, samplesWritten)
-    else
-        @assert samplesWritten == bufLen
-    end
-
-    return buffer
+    return Vector{promote_type(Th, Tx)}(undef, outLen)
 end
 
 
@@ -689,24 +694,34 @@ function resample(x::AbstractVector, rate::AbstractFloat, h::Vector, Nϕ::Intege
     _resample!(x, rate, FIRFilter(h, rate, Nϕ))
 end
 
-function _resample!(x::AbstractVector, rate::Real, self::FIRFilter)
+function _resample!(x::AbstractVector, rate::Real, sf::FIRFilter)
+    undelay!(sf)
+    outLen  = ceil(Int, length(x) * rate)
+    xPadded = _zeropad(x, inputlength(sf, outLen, RoundUp))
+
+    buffer = allocate_output(sf, xPadded)
+    samplesWritten = filt!(buffer, sf, xPadded)
+    return checked_resample_output!(buffer, outLen, samplesWritten, sf)
+end
+
+function undelay!(sf::FIRFilter)
     # Get delay, in # of samples at the output rate, caused by filtering processes
-    τ = timedelay(self)
+    τ = timedelay(sf)
 
     # Use setphase! to
     #   a) adjust the input samples to skip over before producing and output (integer part of τ)
     #   b) set the ϕ index of the PFB (fractional part of τ)
-    setphase!(self, τ)
-
-    # Calculate the number of 0's required
-    outLen      = ceil(Int, length(x) * rate)
-    reqInlen    = inputlength(self, outLen, RoundUp)
-    reqZerosLen = reqInlen - length(x)
-    xPadded     = [x; zeros(eltype(x), reqZerosLen)]
+    setphase!(sf, τ)
+end
 
-    y = filt(self, xPadded)
-    @assert length(y) >= outLen
-    length(y) > outLen && resize!(y, outLen)
+function checked_resample_output!(y::AbstractVector, outLen, samplesWritten, ::FIRFilter{Tk}) where Tk<:FIRKernel
+    if !(Tk <: FIRArbitrary)
+        samplesWritten == length(y) || throw(AssertionError("Length of resampled output different from expectation."))
+    end
+    # outLen: the desired output length ceil(Int, rate * length(input)), but we can overshoot
+    # samplesWritten: number of samples actually written to y; if longer, y[samplesWritten+1:end] contains invalid data
+    samplesWritten >= outLen || throw(AssertionError("Resample output shorter than expected."))
+    length(y) == outLen || resize!(y, outLen)
     return y
 end
 
@@ -742,11 +757,23 @@ end
 resample(x::AbstractArray, rate::Real, args::Real...; dims) =
     _resample!(x, rate, FIRFilter(rate, args...); dims)
 
-_resample!(x::AbstractArray, rate::Real, sf::FIRFilter; dims) =
-    mapslices(x; dims) do v
-        reset!(sf)
-        _resample!(v, rate, sf)
+function _resample!(x::AbstractArray{T}, rate::Real, sf::FIRFilter; dims::Int) where T
+    undelay!(sf)
+    size_v  = size(x, dims)
+    outLen  = ceil(Int, size_v * rate)
+    xPadded = Vector{T}(undef, inputlength(sf, outLen, RoundUp))
+    xPadded[size_v+1:end] .= zero(T)
+    buffer  = allocate_output(sf, xPadded)
+    bufLen  = length(buffer)
+
+    mapslices(x; dims) do v::AbstractVector
+        undelay!(reset!(sf))
+        length(buffer) == bufLen || resize!(buffer, bufLen)
+        copyto!(xPadded, v)
+        samplesWritten = filt!(buffer, sf, xPadded)
+        return checked_resample_output!(buffer, outLen, samplesWritten, sf)
     end
+end
 
 #
 # References

test/FilterTestHelpers.jl

@@ -1,8 +1,11 @@
 module FilterTestHelpers
-using DSP, Test
+using DSP, Test, DelimitedFiles
 
 export tffilter_eq, zpkfilter_eq, tffilter_accuracy, zpkfilter_accuracy,
-       matrix_to_sosfilter, sosfilter_to_matrix
+       matrix_to_sosfilter, sosfilter_to_matrix,
+       read_reference_data
+
+read_reference_data(s, delim='\t') = readdlm(joinpath(@__DIR__, "data", s), delim)
 
 function lt(a, b)
     if abs(real(a) - real(b)) > 1e-10

test/filt.jl

@@ -1,4 +1,4 @@
-!(dirname(@__FILE__) in LOAD_PATH) && push!(LOAD_PATH, dirname(@__FILE__))
+!(@__DIR__() in LOAD_PATH) && push!(LOAD_PATH, @__DIR__)
 using DSP, Test, Random, FilterTestHelpers
 #
 # filt with different filter forms
@@ -181,17 +181,17 @@ end
 #
 ##############
 @testset "filt! with init. cond." begin
-    matlab_filt  = readdlm(joinpath(dirname(@__FILE__), "data", "filt_check.txt"),'\t')
+    matlab_filt = read_reference_data("filt_check.txt")
 
     a = [0.9, 0.6]
     b = [0.4, 1]
     z = [0.4750]
-    x  = vec(readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_x.txt"),'\t'))
+    x = vec(read_reference_data("spectrogram_x.txt"))
     @test_deprecated(filt!(x, b, a, x, z))
 
     @test matlab_filt ≈ x
 
-    x = vec(readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_x.txt"),'\t'))
+    x = vec(read_reference_data("spectrogram_x.txt"))
     filt!(x, DF2TFilter(PolynomialRatio(b, a), z), x)
 
     @test matlab_filt ≈ x
@@ -215,22 +215,22 @@ end
 #
 #######################################
 @testset "filtfilt 1D" begin
-    x2_matlab = readdlm(joinpath(dirname(@__FILE__), "data", "filtfilt_output.txt"),'\t')
+    x2_matlab = read_reference_data("filtfilt_output.txt")
 
     b = [ 0.00327922,  0.01639608,  0.03279216,  0.03279216,  0.01639608,  0.00327922]
     a = [ 1.        , -2.47441617,  2.81100631, -1.70377224,  0.54443269, -0.07231567]
-    x  = readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_x.txt"),'\t')
+    x = read_reference_data("spectrogram_x.txt")
 
     @test x2_matlab ≈ filtfilt(b, a, x)
 end
 
 # Make sure above doesn't crash for real coeffs & complex data.
 @testset "filtfilt 1D Complex" begin
-    x2_matlab = readdlm(joinpath(dirname(@__FILE__), "data", "filtfilt_output.txt"),'\t')
+    x2_matlab = read_reference_data("filtfilt_output.txt")
 
     b = [ 0.00327922,  0.01639608,  0.03279216,  0.03279216,  0.01639608,  0.00327922]
     a = [ 1.        , -2.47441617,  2.81100631, -1.70377224,  0.54443269, -0.07231567]
-    x  = readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_x.txt"),'\t')
+    x = read_reference_data("spectrogram_x.txt")
 
     y = x .+ 1im .* randn(size(x, 1))
 
@@ -260,10 +260,10 @@ end
     b = [ 0.00327922,  0.01639608,  0.03279216,  0.03279216,  0.01639608,  0.00327922]
     a = [ 1.        , -2.47441617,  2.81100631, -1.70377224,  0.54443269, -0.07231567]
 
-    x2_output = readdlm(joinpath(dirname(@__FILE__), "data", "filtfilt_output_2d.txt"),'\t')
+    x2_output = read_reference_data("filtfilt_output_2d.txt")
 
-    x  = readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_x.txt"),'\t')
-    x  = repeat(x, outer=(1, 3))
+    x = read_reference_data("spectrogram_x.txt")
+    x = repeat(x, outer=(1, 3))
     x[:,2] = circshift(x[:,2], 64)
     x[:,3] = circshift(x[:,3], 128)
 
@@ -291,7 +291,7 @@ end
 # the extrapolation will differ slightly.)
 #######################################
 @testset "filtfilt SOS" begin
-    x  = readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_x.txt"),'\t')
+    x = read_reference_data("spectrogram_x.txt")
 
     f = digitalfilter(Lowpass(0.2), Butterworth(4))
     @test filtfilt(convert(SecondOrderSections, f), x) ≈ filtfilt(convert(PolynomialRatio, f), x)

test/filter_conversion.jl

@@ -1,4 +1,4 @@
-!(dirname(@__FILE__) in LOAD_PATH) && push!(LOAD_PATH, dirname(@__FILE__))
+!(@__DIR__() in LOAD_PATH) && push!(LOAD_PATH, @__DIR__)
 using DSP, Test, FilterTestHelpers, Polynomials
 using Polynomials.PolyCompat
 

test/filter_design.jl

@@ -1,7 +1,7 @@
-!(dirname(@__FILE__) in LOAD_PATH) && push!(LOAD_PATH, dirname(@__FILE__))
+!(@__DIR__() in LOAD_PATH) && push!(LOAD_PATH, @__DIR__)
 using DSP, Test, FilterTestHelpers
 using LinearAlgebra: norm
-using DelimitedFiles: readdlm
+using FFTW: fft
 
 #
 # Butterworth filter prototype
@@ -987,29 +987,29 @@ end
 @testset "window FIR" begin
     winfirtaps_jl    = digitalfilter(Lowpass(0.25),FIRWindow(hamming(128), scale=false); fs=1)
     # firwin(128, 0.25, nyq=.5, scale=False)
-    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_128_lowpass_fc0.25_fs1.0.txt"),'\t')
+    winfirtaps_scipy = read_reference_data("digitalfilter_hamming_128_lowpass_fc0.25_fs1.0.txt")
     @test winfirtaps_jl ≈ winfirtaps_scipy
 
     winfirtaps_jl    = digitalfilter(Lowpass(0.25),FIRWindow(hamming(129), scale=false); fs=1)
     # firwin(129, 0.25, nyq=.5, scale=False)
-    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_lowpass_fc0.25_fs1.0.txt"),'\t')
+    winfirtaps_scipy = read_reference_data("digitalfilter_hamming_129_lowpass_fc0.25_fs1.0.txt")
     @test winfirtaps_jl ≈ winfirtaps_scipy
 
     @test_throws ArgumentError digitalfilter(Highpass(0.25),FIRWindow(hamming(128), scale=false))
 
     winfirtaps_jl    = digitalfilter(Highpass(0.25),FIRWindow(hamming(129), scale=false); fs=1)
     # firwin(129, 0.25, nyq=.5, scale=False, pass_zero=False)
-    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_highpass_fc0.25_fs1.0.txt"),'\t')
+    winfirtaps_scipy = read_reference_data("digitalfilter_hamming_129_highpass_fc0.25_fs1.0.txt")
     @test winfirtaps_jl ≈ winfirtaps_scipy
 
     winfirtaps_jl    = digitalfilter(Bandpass(0.1, 0.2),FIRWindow(hamming(128), scale=false); fs=1)
     # firwin(128, [0.1, 0.2], nyq=.5, scale=False, pass_zero=False)
-    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_128_bandpass_fc0.1_0.2_fs1.0.txt"),'\t')
+    winfirtaps_scipy = read_reference_data("digitalfilter_hamming_128_bandpass_fc0.1_0.2_fs1.0.txt")
     @test winfirtaps_jl ≈ winfirtaps_scipy
 
     winfirtaps_jl    = digitalfilter(Bandpass(0.1, 0.2),FIRWindow(hamming(129), scale=false); fs=1)
     # firwin(129, [0.1, 0.2], nyq=.5, scale=False, pass_zero=False)
-    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_bandpass_fc0.1_0.2_fs1.0.txt"),'\t')
+    winfirtaps_scipy = read_reference_data("digitalfilter_hamming_129_bandpass_fc0.1_0.2_fs1.0.txt")
     @test winfirtaps_jl ≈ winfirtaps_scipy
 
     @test_throws ArgumentError ComplexBandpass(2, 1)
@@ -1046,37 +1046,37 @@ end
 
     winfirtaps_jl    = digitalfilter(Bandstop(0.1, 0.2),FIRWindow(hamming(129), scale=false); fs=1)
     # firwin(129, [0.1, 0.2], nyq=.5, scale=False)
-    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_bandstop_fc0.1_0.2_fs1.0.txt"),'\t')
+    winfirtaps_scipy = read_reference_data("digitalfilter_hamming_129_bandstop_fc0.1_0.2_fs1.0.txt")
     @test winfirtaps_jl ≈ vec(winfirtaps_scipy)
 
     winfirtaps_jl    = digitalfilter(Lowpass(0.25),FIRWindow(hamming(128), scale=true); fs=1)
     # firwin(128, 0.25, nyq=.5)
-    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_128_lowpass_scaled_fc0.25_fs1.0.txt"),'\t')
+    winfirtaps_scipy = read_reference_data("digitalfilter_hamming_128_lowpass_scaled_fc0.25_fs1.0.txt")
     @test winfirtaps_jl ≈ winfirtaps_scipy
 
     winfirtaps_jl    = digitalfilter(Lowpass(0.25),FIRWindow(hamming(129), scale=true); fs=1)
     # firwin(129, 0.25, nyq=.5)
-    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_lowpass_scaled_fc0.25_fs1.0.txt"),'\t')
+    winfirtaps_scipy = read_reference_data("digitalfilter_hamming_129_lowpass_scaled_fc0.25_fs1.0.txt")
     @test winfirtaps_jl ≈ winfirtaps_scipy
 
     winfirtaps_jl    = digitalfilter(Highpass(0.25),FIRWindow(hamming(129), scale=true); fs=1)
     # firwin(129, 0.25, nyq=.5, pass_zero=False)
-    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_highpass_scaled_fc0.25_fs1.0.txt"),'\t')
+    winfirtaps_scipy = read_reference_data("digitalfilter_hamming_129_highpass_scaled_fc0.25_fs1.0.txt")
     @test winfirtaps_jl ≈ winfirtaps_scipy
 
     winfirtaps_jl    = digitalfilter(Bandpass(0.1, 0.2),FIRWindow(hamming(128), scale=true); fs=1)
     # firwin(128, [0.1, 0.2], nyq=.5, pass_zero=False)
-    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_128_bandpass_scaled_fc0.1_0.2_fs1.0.txt"),'\t')
+    winfirtaps_scipy = read_reference_data("digitalfilter_hamming_128_bandpass_scaled_fc0.1_0.2_fs1.0.txt")
     @test winfirtaps_jl ≈ winfirtaps_scipy
 
     winfirtaps_jl    = digitalfilter(Bandpass(0.1, 0.2),FIRWindow(hamming(129), scale=true); fs=1)
     # firwin(129, [0.1, 0.2], nyq=.5, scale=False, pass_zero=False)
-    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_bandpass_scaled_fc0.1_0.2_fs1.0.txt"),'\t')
+    winfirtaps_scipy = read_reference_data("digitalfilter_hamming_129_bandpass_scaled_fc0.1_0.2_fs1.0.txt")
     @test winfirtaps_jl ≈ winfirtaps_scipy
 
     winfirtaps_jl    = digitalfilter(Bandstop(0.1, 0.2),FIRWindow(hamming(129), scale=true); fs=1)
     # firwin(129, [0.1, 0.2], nyq=.5, scale=False)
-    winfirtaps_scipy = readdlm(joinpath(dirname(@__FILE__), "data", "digitalfilter_hamming_129_bandstop_scaled_fc0.1_0.2_fs1.0.txt"),'\t')
+    winfirtaps_scipy = read_reference_data("digitalfilter_hamming_129_bandstop_scaled_fc0.1_0.2_fs1.0.txt")
     @test winfirtaps_jl ≈ vec(winfirtaps_scipy)
 end