diff --git a/docs/src/filters.md b/docs/src/filters.md index 8ae6f1510..b938e7356 100644 --- a/docs/src/filters.md +++ b/docs/src/filters.md @@ -59,6 +59,7 @@ parameters. ```@docs FIRFilter +resample_filter ``` ## Filter application @@ -196,7 +197,7 @@ coefa ## Examples Construct a 4th order elliptic lowpass filter with normalized cutoff -frequency 0.2, 0.5 dB of passband ripple, and 30 dB attentuation in +frequency 0.2, 0.5 dB of passband ripple, and 30 dB attenuation in the stopband and extract the coefficients of the numerator and denominator of the transfer function: diff --git a/src/Filters/design.jl b/src/Filters/design.jl index 576335a2f..06eee7e20 100644 --- a/src/Filters/design.jl +++ b/src/Filters/design.jl @@ -124,7 +124,7 @@ function landen(k::Real) kn = Vector{typeof(k)}(undef, niter) # Eq. (50) for i = 1:niter - kn[i] = k = abs2(k/(1+sqrt(1-abs2(k)))) + kn[i] = k = abs2(k / (1 + sqrt(1 - abs2(k)))) end kn end @@ -223,7 +223,7 @@ end Elliptic(n::Integer, rp::Real, rs::Real) `n` pole elliptic (Cauer) filter with `rp` dB ripple in the -passband and `rs` dB attentuation in the stopband. +passband and `rs` dB attenuation in the stopband. """ Elliptic(n::Integer, rp::Real, rs::Real) = Elliptic(Float64, n, rp, rs) @@ -473,7 +473,7 @@ Here, `m` and `n` are respectively the numbers of zeros and poles in the s-domai then additional `n-m` zeros are added at `z = -1`. """ function bilinear(f::ZeroPoleGain{:s,Z,P,K}, fs::Real) where {Z,P,K} - ztype = typeof(0 + zero(Z)/fs) + ztype = typeof(0 + zero(Z) / fs) z = fill(convert(ztype, -1), max(length(f.p), length(f.z))) ptype = typeof(0 + zero(P) / fs) @@ -481,17 +481,17 @@ function bilinear(f::ZeroPoleGain{:s,Z,P,K}, fs::Real) where {Z,P,K} num = one(one(fs) - one(Z)) for i = 1:length(f.z) - z[i] = (2 + f.z[i] / fs)/(2 - f.z[i] / fs) + z[i] = (2 + f.z[i] / fs) / (2 - f.z[i] / fs) num *= (2 * fs - f.z[i]) end den = one(one(fs) - one(P)) for i = 1:length(f.p) - p[i] = (2 + f.p[i] / fs)/(2 - f.p[i]/fs) + p[i] = (2 + f.p[i] / fs) / (2 - f.p[i] / fs) den *= (2 * fs - f.p[i]) end - ZeroPoleGain{:z}(z, p, f.k * real(num)/real(den)) + ZeroPoleGain{:z}(z, p, f.k * real(num) / real(den)) end # Pre-warp filter frequencies for digital filtering @@ -545,12 +545,12 @@ end # Get length and alpha for Kaiser window filter with specified # transition band width and stopband attenuation in dB function kaiserord(transitionwidth::Real, attenuation::Real=60) - n = ceil(Int, (attenuation - 7.95)/(π*2.285*transitionwidth))+1 + n = ceil(Int, (attenuation - 7.95) / (π * 2.285 * transitionwidth)) + 1 if attenuation > 50 - β = 0.1102*(attenuation - 8.7) + β = 0.1102 * (attenuation - 8.7) elseif attenuation >= 21 - β = 0.5842*(attenuation - 21)^0.4 + 0.07886*(attenuation - 21) + β = 0.5842 * (attenuation - 21)^0.4 + 0.07886 * (attenuation - 21) else β = 0.0 end @@ -670,50 +670,47 @@ function digitalfilter(ftype::FilterType, proto::FIRWindow; fs::Real=2) coefs = firprototype(length(proto.window), ftype, fs) @assert length(proto.window) == length(coefs) out = (coefs .*= proto.window) - proto.scale ? rmul!(out, 1/scalefactor(out, ftype, fs)) : out + proto.scale ? rmul!(out, 1 / scalefactor(out, ftype, fs)) : out end -# Compute FIR coefficients necessary for arbitrary rate resampling -function resample_filter(rate::AbstractFloat, Nϕ = 32, rel_bw = 1.0, attenuation = 60) - f_nyq = rate >= 1.0 ? 1.0/Nϕ : rate/Nϕ - cutoff = f_nyq * rel_bw - trans_width = cutoff * 0.2 - - # Determine resampling filter order - hLen, α = kaiserord(trans_width, attenuation) +""" + resample_filter(rate::AbstractFloat, Nϕ::Integer=32, rel_bw::Real=1.0, att=60) - # Round the number of taps up to a multiple of Nϕ. - # Otherwise the missing taps will be filled with 0. - hLen = Nϕ * ceil(Int, hLen/Nϕ) +Compute FIR coefficients suitable for arbitrary rate resampling +with `Nϕ` filters / phases, stop-band attenuation `att`, and relative bandwidth `rel_bw`. +""" +function resample_filter(rate::AbstractFloat, Nϕ::Integer=32, rel_bw::Real=1.0, attenuation=60) + f_nyq = rate >= 1.0 ? 1.0 / Nϕ : rate / Nϕ + return _resample_filter(f_nyq, Nϕ, rel_bw, attenuation) +end - # Ensure that the filter is an odd length - if (iseven(hLen)) - hLen += 1 - end +""" + resample_filter(rate::Union{Integer,Rational}, rel_bw::Real=1.0, attenuation=60) - # Design filter - h = digitalfilter(Lowpass(cutoff), FIRWindow(kaiser(hLen, α))) - rmul!(h, Nϕ) +Compute FIR coefficients suitable for rational rate resampling, +with stop-band attenuation `att`, and relative bandwidth `rel_bw`. +""" +function resample_filter(rate::Union{Integer,Rational}, rel_bw::Real=1.0, attenuation=60) + Nϕ = numerator(rate) + decimation = denominator(rate) + f_nyq = min(1 / Nϕ, 1 / decimation) + return _resample_filter(f_nyq, Nϕ, rel_bw, attenuation) end -# Compute FIR coefficients necessary for rational rate resampling -function resample_filter(rate::Union{Integer,Rational}, rel_bw = 1.0, attenuation = 60) - Nϕ = numerator(rate) - decimation = denominator(rate) - f_nyq = min(1/Nϕ, 1/decimation) +function _resample_filter(f_nyq, Nϕ, rel_bw, attenuation) cutoff = f_nyq * rel_bw trans_width = cutoff * 0.2 # Determine resampling filter order hLen, α = kaiserord(trans_width, attenuation) - # Round the number of taps up to a multiple of Nϕ (same as interpolation factor). + # Round the number of taps up to a multiple of Nϕ. # Otherwise the missing taps will be filled with 0. - hLen = Nϕ * ceil(Int, hLen/Nϕ) + hLen = Nϕ * ceil(Int, hLen / Nϕ) # Ensure that the filter is an odd length - if (iseven(hLen)) + if iseven(hLen) hLen += 1 end diff --git a/src/Filters/filt.jl b/src/Filters/filt.jl index 0765181c5..aa214650a 100644 --- a/src/Filters/filt.jl +++ b/src/Filters/filt.jl @@ -340,10 +340,10 @@ end # Zero phase digital filtering for second order sections function filtfilt(f::SecondOrderSections{:z,T,G}, x::AbstractArray{S}) where {T,G,S} zi = filt_stepstate(f) - pad_length = 6 * length(f.biquads) + pad_length = min(6 * length(f.biquads), size(x, 1) - 1) t = promote_type(T, G, S) zitmp = similar(zi, t) - extrapolated = Vector{t}(undef, size(x, 1)+pad_length*2) + extrapolated = Vector{t}(undef, size(x, 1) + pad_length * 2) out = similar(x, t) for col in CartesianIndices(axes(x)[2:end]) diff --git a/src/Filters/stream_filt.jl b/src/Filters/stream_filt.jl index 7fb0fecf0..9dee1d20f 100644 --- a/src/Filters/stream_filt.jl +++ b/src/Filters/stream_filt.jl @@ -3,7 +3,7 @@ const PFB{T} = Matrix{T} # polyphase filter bank abstract type FIRKernel{T} end # Single rate FIR kernel -mutable struct FIRStandard{T} <: FIRKernel{T} +struct FIRStandard{T} <: FIRKernel{T} h::Vector{T} hLen::Int end @@ -17,96 +17,92 @@ end # Interpolator FIR kernel mutable struct FIRInterpolator{T} <: FIRKernel{T} - pfb::PFB{T} - interpolation::Int - Nϕ::Int - tapsPerϕ::Int + const pfb::PFB{T} + const interpolation::Int + const Nϕ::Int + const tapsPerϕ::Int + const hLen::Int inputDeficit::Int ϕIdx::Int - hLen::Int end -function FIRInterpolator(h::Vector, interpolation::Integer) - pfb = taps2pfb(h, interpolation) - tapsPerϕ, Nϕ = size(pfb) - interpolation = interpolation - inputDeficit = 1 - ϕIdx = 1 - hLen = length(h) +function FIRInterpolator(h::Vector, interpolation::Int) + pfb = taps2pfb(h, interpolation) + tapsPerϕ, Nϕ = size(pfb) + inputDeficit = 1 + ϕIdx = 1 + hLen = length(h) - FIRInterpolator(pfb, interpolation, Nϕ, tapsPerϕ, inputDeficit, ϕIdx, hLen) + FIRInterpolator(pfb, interpolation, Nϕ, tapsPerϕ, hLen, inputDeficit, ϕIdx) end - +FIRInterpolator(h::Vector, interpolation::Integer) = FIRInterpolator(h, Int(interpolation)) # Decimator FIR kernel mutable struct FIRDecimator{T} <: FIRKernel{T} - h::Vector{T} - hLen::Int - decimation::Int + const h::Vector{T} + const hLen::Int + const decimation::Int inputDeficit::Int end -function FIRDecimator(h::Vector, decimation::Integer) +function FIRDecimator(h::Vector, decimation::Int) h = reverse(h) hLen = length(h) inputDeficit = 1 FIRDecimator(h, hLen, decimation, inputDeficit) end - +FIRDecimator(h::Vector, decimation::Integer) = FIRDecimator(h, Int(decimation)) # Rational resampler FIR kernel mutable struct FIRRational{T} <: FIRKernel{T} - pfb::PFB{T} - ratio::Rational{Int} - Nϕ::Int - ϕIdxStepSize::Int - tapsPerϕ::Int + const pfb::PFB{T} + const ratio::Rational{Int} + const Nϕ::Int + const ϕIdxStepSize::Int + const tapsPerϕ::Int + const hLen::Int ϕIdx::Int inputDeficit::Int - hLen::Int end -function FIRRational(h::Vector, ratio::Rational) +function FIRRational(h::Vector, ratio::Rational{Int}) pfb = taps2pfb(h, numerator(ratio)) tapsPerϕ, Nϕ = size(pfb) ϕIdxStepSize = mod(denominator(ratio), numerator(ratio)) ϕIdx = 1 inputDeficit = 1 hLen = length(h) - FIRRational(pfb, ratio, Nϕ, ϕIdxStepSize, tapsPerϕ, ϕIdx, inputDeficit, hLen) + FIRRational(pfb, ratio, Nϕ, ϕIdxStepSize, tapsPerϕ, hLen, ϕIdx, inputDeficit) end - -FIRRational(h::Vector,ratio::Integer)=FIRRational(h,convert(Rational,ratio)) +FIRRational(h::Vector, ratio::Union{Integer,Rational}) = FIRRational(h, convert(Rational{Int}, ratio)) # # Arbitrary resampler FIR kernel # -# This kernel is different from the others in that it has two polyphase filtlter banks. -# The the second filter bank, dpfb, is the derivative of pfb. The purpose of this is to -# allow us to compute two y values, yLower & yUpper, whitout having to advance the input -# index by 1. It makes the kernel simple by not having to store extra state in the case -# when where's at the last polphase branch and the last available input sample. By using -# a derivitive filter, we can always compute the output in that scenario. +# This kernel is different from the others in that it has two polyphase filter banks. +# The second filter bank, dpfb, is the derivative of pfb. The purpose of this is to +# allow us to compute two y values, yLower & yUpper, without having to advance the input +# index by 1. It makes the kernel simpler by not having to store extra state in the case +# when where's at the last polyphase branch and the last available input sample. By using +# a derivative filter, we can always compute the output in that scenario. # See section 7.6.1 in [1] for a better explanation. mutable struct FIRArbitrary{T} <: FIRKernel{T} - rate::Float64 - pfb::PFB{T} - dpfb::PFB{T} - Nϕ::Int - tapsPerϕ::Int + const rate::Float64 + const pfb::PFB{T} + const dpfb::PFB{T} + const Nϕ::Int + const tapsPerϕ::Int ϕAccumulator::Float64 ϕIdx::Int α::Float64 - Δ::Float64 + const Δ::Float64 inputDeficit::Int xIdx::Int - hLen::Int + const hLen::Int end -function FIRArbitrary(h::Vector, rate_in::Real, Nϕ_in::Integer) - rate = convert(Float64, rate_in) - Nϕ = convert(Int, Nϕ_in) +function FIRArbitrary(h::Vector, rate::Float64, Nϕ::Int) dh = [diff(h); zero(eltype(h))] pfb = taps2pfb(h, Nϕ) dpfb = taps2pfb(dh, Nϕ) @@ -133,7 +129,7 @@ function FIRArbitrary(h::Vector, rate_in::Real, Nϕ_in::Integer) hLen ) 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} @@ -145,12 +141,15 @@ end # Constructor for single-rate, decimating, interpolating, and rational resampling filters """ - FIRFilter(h::Vector[, ratio::Union{Integer,Rational}]) + FIRFilter(h::Vector, ratio::Union{Integer,Rational}=1) + FIRFilter(ratio::Union{Integer,Rational}, args...) Construct a stateful FIRFilter object from the vector of filter taps `h`. `ratio` is an optional rational integer which specifies the input to output sample rate relationship (e.g. `147//160` for converting recorded audio from 48 KHz to 44.1 KHz). + +Constructs `h` with `resample_filter(ratio, args...)` if it is not provided. """ function FIRFilter(h::Vector, resampleRatio::Union{Integer,Rational} = 1) interpolation = numerator(resampleRatio) @@ -160,12 +159,12 @@ function FIRFilter(h::Vector, resampleRatio::Union{Integer,Rational} = 1) if resampleRatio == 1 # single-rate kernel = FIRStandard(h) historyLen = kernel.hLen - 1 - elseif interpolation == 1 # decimate - kernel = FIRDecimator(h, decimation) - historyLen = kernel.hLen - 1 elseif decimation == 1 # interpolate kernel = FIRInterpolator(h, interpolation) historyLen = kernel.tapsPerϕ - 1 + elseif interpolation == 1 # decimate + kernel = FIRDecimator(h, decimation) + historyLen = kernel.hLen - 1 else # rational kernel = FIRRational(h, resampleRatio) historyLen = kernel.tapsPerϕ - 1 @@ -178,13 +177,16 @@ end # Constructor for arbitrary resampling filter (polyphase interpolator w/ intra-phase linear interpolation) """ - FIRFilter(h::Vector, rate::AbstractFloat[, Nϕ::Integer=32]) + FIRFilter(h::Vector, rate::AbstractFloat, Nϕ::Integer=32) + FIRFilter(rate::AbstractFloat, Nϕ::Integer=32, args...) Returns a polyphase FIRFilter object from the vector of filter taps `h`. `rate` is a floating point number that specifies the input to output sample-rate relationship ``\\frac{fs_{out}}{fs_{in}}``. `Nϕ` is an optional parameter which specifies the number of *phases* created from `h`. `Nϕ` defaults to 32. + +Constructs `h` with `resample_filter(rate, Nϕ, args...)` if it is not provided. """ function FIRFilter(h::Vector, rate::AbstractFloat, Nϕ::Integer=32) rate > 0.0 || throw(DomainError(rate, "rate must be greater than 0")) @@ -195,14 +197,14 @@ function FIRFilter(h::Vector, rate::AbstractFloat, Nϕ::Integer=32) end # Constructor for a resampling FIR filter, where the user needs only to set the sampling rate -function FIRFilter(rate::AbstractFloat, Nϕ::Integer=32) - h = resample_filter(rate, Nϕ) - FIRFilter(h, rate) +function FIRFilter(rate::AbstractFloat, Nϕ::Integer=32, args...) + h = resample_filter(rate, Nϕ, args...) + FIRFilter(h, rate, Nϕ) end -function FIRFilter(rate::Union{Integer,Rational}) - h = resample_filter(rate) - FIRFilter(h, rate) +function FIRFilter(ratio::Union{Integer,Rational}, args...) + h = resample_filter(ratio, args...) + FIRFilter(h, ratio) end @@ -290,14 +292,13 @@ end function taps2pfb(h::Vector{T}, Nϕ::Integer) where T hLen = length(h) - tapsPerϕ = ceil(Int, hLen/Nϕ) - pfbSize = tapsPerϕ * Nϕ + tapsPerϕ = ceil(Int, hLen / Nϕ) pfb = Matrix{T}(undef, tapsPerϕ, Nϕ) hIdx = 1 for rowIdx in tapsPerϕ:-1:1, colIdx in 1:Nϕ tap = hIdx > hLen ? zero(T) : h[hIdx] - @inbounds pfb[rowIdx,colIdx] = tap + pfb[rowIdx, colIdx] = tap hIdx += 1 end @@ -414,10 +415,10 @@ function filt!(buffer::AbstractVector{Tb}, self::FIRFilter{FIRStandard{Th}}, x:: h = kernel.h for i = 1:min(kernel.hLen-1, xLen) - @inbounds buffer[i] = unsafe_dot(h, history, x, i) + buffer[i] = unsafe_dot(h, history, x, i) end for i = kernel.hLen:xLen - @inbounds buffer[i] = unsafe_dot(h, x, i) + buffer[i] = unsafe_dot(h, x, i) end self.history = shiftin!(history, x) @@ -545,8 +546,8 @@ function filt!(buffer::AbstractVector{Tb}, self::FIRFilter{FIRDecimator{Th}}, x: accumulator = unsafe_dot(kernel.h, x, inputIdx) end - @inbounds buffer[bufIdx] = accumulator - inputIdx += kernel.decimation + buffer[bufIdx] = accumulator + inputIdx += kernel.decimation end kernel.inputDeficit = inputIdx - xLen @@ -680,9 +681,15 @@ zeros to `x`. The result is that when the input and output signals are plotted on top of each other, they correlate very well, but one signal will have more samples than the other. """ -function resample(x::AbstractVector, rate::Real, h::Vector) - self = FIRFilter(h, rate) +function resample(x::AbstractVector, rate::Union{Integer,Rational}, h::Vector) + _resample!(x, rate, FIRFilter(h, rate)) +end + +function resample(x::AbstractVector, rate::AbstractFloat, h::Vector, Nϕ::Integer=32) + _resample!(x, rate, FIRFilter(h, rate, Nϕ)) +end +function _resample!(x::AbstractVector, rate::Real, self::FIRFilter) # Get delay, in # of samples at the output rate, caused by filtering processes τ = timedelay(self) @@ -692,7 +699,7 @@ function resample(x::AbstractVector, rate::Real, h::Vector) setphase!(self, τ) # Calculate the number of 0's required - outLen = ceil(Int, length(x)*rate) + outLen = ceil(Int, length(x) * rate) reqInlen = inputlength(self, outLen, RoundUp) reqZerosLen = reqInlen - length(x) xPadded = [x; zeros(eltype(x), reqZerosLen)] @@ -703,21 +710,43 @@ function resample(x::AbstractVector, rate::Real, h::Vector) return y end -function resample(x::AbstractVector, rate::Real) - h = resample_filter(rate) - resample(x, rate, h) +""" + resample(x::AbstractVector, rate::Real, args::Real...) + +Constructs a filter with `resample_filter` using the optional arguments `args`, +and resamples the signal `x` with it. +""" +function resample(x::AbstractVector, rate::Real, args::Real...) + sf = FIRFilter(rate, args...) + _resample!(x, rate, sf) end """ resample(x::AbstractArray, rate::Real, h::Vector = resample_filter(rate); dims) + resample(x::AbstractArray, rate::AbstractFloat, h::Vector, Nϕ=32; dims) Resample an array `x` along dimension `dims`. +If `rate` is an `AbstractFloat`, the number of phases `Nϕ` +to be used in constructing `FIRArbitrary` can be supplied +as an optional argument, which defaults to 32. """ -function resample(x::AbstractArray, rate::Real, h::Vector = resample_filter(rate); dims) - mapslices(x; dims) do x - resample(x, rate, h) - end +function resample(x::AbstractArray, rate::Union{Integer,Rational}, h::Vector; dims) + sf = FIRFilter(h, rate) + return _resample!(x, rate, sf; dims) end +function resample(x::AbstractArray, rate::AbstractFloat, h::Vector, Nϕ=32; dims) + sf = FIRFilter(h, rate, Nϕ) + _resample!(x, rate, sf; dims) +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) + end # # References diff --git a/test/filt.jl b/test/filt.jl index 77239272f..f028e4879 100644 --- a/test/filt.jl +++ b/test/filt.jl @@ -301,6 +301,7 @@ end f = digitalfilter(Bandpass(0.1, 0.3), Butterworth(2)) @test filtfilt(convert(SecondOrderSections, f), x) ≈ filtfilt(convert(PolynomialRatio, f), x) + @test_nowarn filtfilt(f, [1.]) # check that pad_length won't throw oob end # diff --git a/test/filt_stream.jl b/test/filt_stream.jl index f54470f60..0d8f44f82 100644 --- a/test/filt_stream.jl +++ b/test/filt_stream.jl @@ -179,9 +179,9 @@ function test_interpolation(h::AbstractVector{T}, x::AbstractVector{V}, interpol x2 = x[pivotPoint+1:end] @printfifinteractive( """\n\n - _ _ _ ___ ____ ____ ___ _ ____ ____ ___ _ ____ _ _ - | |\\ | | |___ |__/ |__] | | | |__| | | | | |\\ | - | | \\| | |___ | \\ | |___ |__| | | | | |__| | \\|\n\n""" ) + _ _ _ ___ ____ ____ ___ ____ _ ____ ___ _ ____ _ _ + | |\\ | | |___ |__/ |__] | | | |__| | | | | |\\ | + | | \\| | |___ | \\ | |__| |__ | | | | |__| | \\|\n\n""") @printfifinteractive( "Testing interpolation, h::%s, x::%s. xLen = %d, hLen = %d, interpolation = %d\n", typeof(h), typeof(x), xLen, hLen, interpolation ) @printfifinteractive( "\n\tNaive interpolation with filt\n\t\t") diff --git a/test/resample.jl b/test/resample.jl index 202496815..e28063cae 100644 --- a/test/resample.jl +++ b/test/resample.jl @@ -2,64 +2,55 @@ using DSP using Test using DelimitedFiles: readdlm +reference_data(s) = vec(readdlm(joinpath(dirname(@__FILE__), "data", s), '\t')) + @testset "rational ratio" begin # AM Modulator # sig(t) = [(1 + sin(2π*0.005*t)) * sin(2π*.05*t) for t in t] - x_ml = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_x.txt"),'\t')) - - # - # [y1,b1] = resample(x, 1, 2) - # - rate = 1//2 - h1_ml = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_taps_1_2.txt"),'\t')) - y1_ml = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_y_1_2.txt"),'\t')) - y1_jl = resample(x_ml, rate, h1_ml) - @test y1_jl ≈ y1_ml - - - # - # [y2,b2] = resample(x, 2, 1) - # - rate = 2//1 - h2_ml = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_taps_2_1.txt"),'\t')) - y2_ml = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_y_2_1.txt"),'\t')) - y2_jl = resample(x_ml, rate, h2_ml) - @test y2_jl ≈ y2_ml - - - # - # [y3,b3] = resample(x, 3, 2) - # - rate = 3//2 - h3_ml = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_taps_3_2.txt"),'\t')) - y3_ml = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_y_3_2.txt"),'\t')) - y3_jl = resample(x_ml, rate, h3_ml) - @test y3_jl ≈ y3_ml - + x_ml = reference_data("resample_x.txt") # - # [y4,b4] = resample(x, 2, 3) + # [y,b] = resample(x, $num, $den) # - rate = 2//3 - h4_ml = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_taps_2_3.txt"),'\t')) - y4_ml = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_y_2_3.txt"),'\t')) - y4_jl = resample(x_ml, rate, h4_ml) - @test y4_jl ≈ y4_ml + for rate in (1//2, 2//1, 3//2, 2//3) + num, den = numerator(rate), denominator(rate) + h1_ml = reference_data("resample_taps_$(num)_$(den).txt") + y1_ml = reference_data("resample_y_$(num)_$(den).txt") + y1_jl = resample(x_ml, rate, h1_ml) + @test y1_ml ≈ y1_jl + @test y1_ml ≈ resample(x_ml, rate) rtol = 0.001 # check default taps are ok + end end @testset "array signal" begin - rate = 1//2 - x_ml = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_x.txt"),'\t')) - h1_ml = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_taps_1_2.txt"),'\t')) - y1_ml = vec(readdlm(joinpath(dirname(@__FILE__), "data", "resample_y_1_2.txt"),'\t')) - - expected_result = [y1_ml 2y1_ml] - X = [x_ml 2x_ml] - y1_jl = resample(X, rate, h1_ml, dims=1) + x, mat, rate = rand(10_000), rand(121, 212), 1.23456789 + h = resample_filter(rate) + # check that filtering with taps from `resample_filter` give equal results + @test resample(x, rate) == resample(x, rate, h) + @test resample(x, rate, 64, 0.8) == resample(x, rate, resample_filter(rate, 64, 0.8), 64) + @test resample(mat, rate; dims=1) == resample(mat, rate, h; dims=1) + @test resample(mat, rate; dims=2) == resample(mat, rate, h; dims=2) + + rate = 1//2 + x_ml = reference_data("resample_x.txt") + h1_ml = reference_data("resample_taps_1_2.txt") + y1_ml = reference_data("resample_y_1_2.txt") + expected_result = [y1_ml ℯ * y1_ml] + X = [x_ml ℯ * x_ml] + + y1_jl = resample(X, rate, h1_ml; dims=1) + arb_y1 = resample(X, float(rate); dims=1) + rat_y1 = resample(X, rate; dims=1) @test y1_jl ≈ expected_result + @test arb_y1 ≈ expected_result rtol = 0.002 # test that default taps are good enough + @test rat_y1 ≈ expected_result rtol = 0.0005 - y1_jl = resample(X', rate, h1_ml, dims=2) + y1_jl = resample(X', rate, h1_ml; dims=2) + arb_y1 = resample(X', float(rate); dims=2) + rat_y1 = resample(X', rate; dims=2) @test y1_jl ≈ expected_result' + @test arb_y1 ≈ expected_result' rtol = 0.002 + @test rat_y1 ≈ expected_result' rtol = 0.0005 expected_result_3d = permutedims(reshape(expected_result, (size(expected_result, 1), size(expected_result, 2), 1)), (3, 1, 2)) X_3d = permutedims(reshape(X, (size(X, 1), size(X, 2), 1)), (3, 1, 2)) @@ -183,3 +174,20 @@ end @test outputlength(H, inputlength(H, yL, RoundUp)-1) < yL <= outputlength(H, inputlength(H, yL, RoundUp)) end end + +@testset "FIRFilter types" begin + using DSP.Filters: FIRStandard, FIRInterpolator + # test FIRRational(::Vector, ::Int(/Int32)) inferred result type + FIRFutyp = Union{FIRFilter{<:FIRInterpolator},FIRFilter{<:FIRStandard}} + @test only(Base.return_types(FIRFilter, (Vector, Int64))) <: FIRFutyp broken=VERSION