Skip to content

Commit

Permalink
Buffer array resample and reference_data (#623)
Browse files Browse the repository at this point in the history
  • Loading branch information
@wheeheee
wheeheee authored Feb 6, 2025
1 parent cd2912c commit 2257010
Show file tree
Hide file tree
Showing 12 changed files with 203 additions and 153 deletions.
2 changes: 1 addition & 1 deletion Project.toml
View file
Original file line number Diff line number Diff line change
@@ -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"
Expand Down
99 changes: 63 additions & 36 deletions src/Filters/stream_filt.jl
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,5 @@
using ..DSP: _zeropad

const PFB{T} = Matrix{T} # polyphase filter bank

abstract type FIRKernel{T} end
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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


Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down
7 changes: 5 additions & 2 deletions test/FilterTestHelpers.jl
View file
Original file line number Diff line number Diff line change
@@ -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
Expand Down
24 changes: 12 additions & 12 deletions test/filt.jl
View file
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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))

Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion test/filter_conversion.jl
View file
Original file line number Diff line number Diff line change
@@ -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

Expand Down
28 changes: 14 additions & 14 deletions test/filter_design.jl
View file
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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

Expand Down
Loading

2 comments on commit 2257010

@martinholters
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@JuliaRegistrator register

@JuliaRegistrator
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Registration pull request created: JuliaRegistries/General/124439

Tip: Release Notes

Did you know you can add release notes too? Just add markdown formatted text underneath the comment after the text
"Release notes:" and it will be added to the registry PR, and if TagBot is installed it will also be added to the
release that TagBot creates. i.e.

@JuliaRegistrator register

Release notes:

## Breaking changes

- blah

To add them here just re-invoke and the PR will be updated.

Tagging

After the above pull request is merged, it is recommended that a tag is created on this repository for the registered package version.

This will be done automatically if the Julia TagBot GitHub Action is installed, or can be done manually through the github interface, or via:

git tag -a v0.8.1 -m "<description of version>" 225701045a806ff4df1a292ab2728ca7580a4eee
git push origin v0.8.1

Please sign in to comment.