From 4ffeb623a4e46bf186b824475cdc92390d529471 Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Wed, 9 Oct 2024 15:03:13 +0800 Subject: [PATCH 1/8] iterate over x one_buffer --- test/filt_stream.jl | 35 ++++++++++++++++++++--------------- 1 file changed, 20 insertions(+), 15 deletions(-) diff --git a/test/filt_stream.jl b/test/filt_stream.jl index 8aa0d17b3..03fc991a5 100644 --- a/test/filt_stream.jl +++ b/test/filt_stream.jl @@ -13,7 +13,7 @@ function naivefilt(h::Vector, x::Vector{T}, resamplerate::Union{Integer, Rationa end y = filt(h, one(T), xZeroStuffed) - y = [y[n] for n = 1:downfactor:length(y)] + return y[1:downfactor:length(y)] end @@ -87,7 +87,7 @@ function test_singlerate(h::AbstractVector{T}, x::AbstractVector) where T @timeifinteractive naiveResult = filt(h, one(T), x) @printfifinteractive( "\n\tfilt( h, x, 1//1 )\n\t\t" ) - @timeifinteractive statelessResult = filt( h, x ) + @timeifinteractive statelessResult = filt(h, x) @test naiveResult ≈ statelessResult @printfifinteractive( "\n\tfilt. length(x1) = %d, length(x2) = %d\n\t\t", length(x1), length(x2) ) @@ -102,8 +102,9 @@ function test_singlerate(h::AbstractVector{T}, x::AbstractVector) where T @printfifinteractive( "\n\tfilt filt. Piecewise for first %d inputs\n\t\t", length(x1) ) reset!(myfilt) @timeifinteractive begin - for i in 1:length(x1) - y1[i] = filt(myfilt, x1[i:i])[1] + one_buffer = similar(x1, 1) + for (i, p) in enumerate(x1) + y1[i] = filt(myfilt, setindex!(one_buffer, p, 1))[1] end y2 = filt(myfilt, x2) end @@ -150,10 +151,11 @@ function test_decimation(h, x, decimation) @printfifinteractive( "\n\tfilt decimation. Piecewise for first %d inputs.\n\t\t", length(x1) ) reset!(myfilt) - y1 = similar( x, 0 ) + y1 = similar(x, 0) @timeifinteractive begin - for i in 1:length(x1) - append!(y1, filt(myfilt, x1[i:i])) + one_buffer = similar(x1, 1) + for p in x1 + append!(y1, filt(myfilt, setindex!(one_buffer, p, 1))) end y2 = filt(myfilt, x2) end @@ -192,11 +194,11 @@ function test_interpolation(h::AbstractVector{T}, x::AbstractVector{V}, interpol end @printfifinteractive( "\n\tfilt( h, x, %d//1 )\n\t\t", interpolation ) - @timeifinteractive statelessResult = filt( h, x, interpolation//1 ) + @timeifinteractive statelessResult = filt(h, x, interpolation//1) @test naiveResult ≈ statelessResult @printfifinteractive( "\n\tfilt interpolation. length(x1) = %d, length(x2) = %d\n\t\t", length(x1), length(x2) ) - myfilt = FIRFilter( h, interpolation//1 ) + myfilt = FIRFilter(h, interpolation//1) @timeifinteractive begin y1 = filt(myfilt, x1) y2 = filt(myfilt, x2) @@ -211,8 +213,9 @@ function test_interpolation(h::AbstractVector{T}, x::AbstractVector{V}, interpol reset!(myfilt) y1 = similar(x, 0) @timeifinteractive begin - for i in 1:length(x1) - append!(y1, filt(myfilt, x1[i:i])) + one_buffer = similar(x1, 1) + for p in x1 + append!(y1, filt(myfilt, setindex!(one_buffer, p, 1))) end y2 = filt(myfilt, x2) end @@ -265,8 +268,9 @@ function test_rational(h, x, ratio) reset!(myfilt) y1 = similar(x, 0) @timeifinteractive begin - for i in 1:length(x) - append!(y1, filt(myfilt, x[i:i])) + one_buffer = similar(x1, 1) + for p in x + append!(y1, filt(myfilt, setindex!(one_buffer, p, 1))) end end piecewiseResult = y1 @@ -312,8 +316,9 @@ function test_arbitrary(Th, x, resampleRate, numFilters) reset!(myfilt) piecewiseResult = similar(x, 0) sizehint!(piecewiseResult, ceil(Int, length(x)*resampleRate)) - @timeifinteractive for i in 1:length(x) - thisY = filt(myfilt, x[i:i]) + one_buffer = similar(x, 1) + @timeifinteractive for p in x + thisY = filt(myfilt, setindex!(one_buffer, p, 1)) append!(piecewiseResult, thisY) end From 4a2852f8136085b9f9b80903adf6e7973b408945 Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Wed, 9 Oct 2024 20:24:07 +0800 Subject: [PATCH 2/8] improve tests --- test/filt_stream.jl | 32 ++++++++++++++------------------ 1 file changed, 14 insertions(+), 18 deletions(-) diff --git a/test/filt_stream.jl b/test/filt_stream.jl index 03fc991a5..38a617a2c 100644 --- a/test/filt_stream.jl +++ b/test/filt_stream.jl @@ -120,7 +120,7 @@ end # Decimation # -function test_decimation(h, x, decimation) +function test_decimation(h, x, decimation::Integer) xLen = length(x) hLen = length(h) pivotPoint = min(rand(50:150), xLen ÷ 4) @@ -171,7 +171,7 @@ end # Interpolation # -function test_interpolation(h::AbstractVector{T}, x::AbstractVector{V}, interpolation) where {T,V} +function test_interpolation(h::AbstractVector{T}, x::AbstractVector{V}, interpolation::Integer) where {T,V} xLen = length(x) hLen = length(h) pivotPoint = min(rand(50:150), xLen ÷ 4) @@ -285,7 +285,7 @@ end # Arbitrary resampling # -function test_arbitrary(Th, x, resampleRate, numFilters) +function test_arbitrary(Th, x, resampleRate, numFilters::Integer) cutoffFreq = 0.45 transitionWidth = 0.05 h = digitalfilter(Lowpass(cutoffFreq), FIRWindow(transitionwidth=transitionWidth/numFilters); fs=numFilters) .* numFilters @@ -322,11 +322,9 @@ function test_arbitrary(Th, x, resampleRate, numFilters) append!(piecewiseResult, thisY) end - commonLen = minimum(length.((naiveResult, statelessResult, statefulResult, piecewiseResult))) - resize!(naiveResult, commonLen) - resize!(statelessResult, commonLen) - resize!(statefulResult, commonLen) - resize!(piecewiseResult, commonLen) + results = (naiveResult, statelessResult, statefulResult, piecewiseResult) + commonLen = minimum(length, results) + resize!.(results, commonLen) @test naiveResult ≈ statelessResult @test naiveResult ≈ statefulResult @@ -347,24 +345,22 @@ end xLen = xLen-mod(xLen, decimation) x = rand(Tx, xLen) ratio = interpolation//decimation - if ratio == 1 + if isone(ratio) test_singlerate(h, x) + elseif numerator(ratio) == interpolation + test_rational(h, x, ratio) + if Tx in (Float32, ComplexF32) + test_arbitrary(Th, x, convert(Float64, ratio) + rand(), 32) + end end if decimation != 1 test_decimation(h, x, decimation) + else + test_rational(h, x, interpolation) end if interpolation != 1 test_interpolation(h, x, interpolation) end - if numerator(ratio) == interpolation && denominator(ratio) == decimation && numerator(ratio) != 1 && denominator(ratio) != 1 - test_rational(h, x, ratio) - if Tx in [Float32, ComplexF32] - test_arbitrary(Th, x, convert(Float64, ratio)+rand(), 32) - end - end - if decimation == 1 - test_rational(h, x, interpolation) - end end @test resample(1:2, 3, [zeros(2); 1; zeros(3)]) == [1, 0, 0, 2, 0, 0] From f299a1f66d4911e7b365c4aea4de1d943b125336 Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Sat, 12 Oct 2024 12:33:44 +0800 Subject: [PATCH 3/8] jacobsen, sinpi, cospi --- test/estimation.jl | 81 ++++++++++++++++++---------------------------- test/multitaper.jl | 14 ++++---- test/remez_fir.jl | 8 ++--- test/util.jl | 2 +- 4 files changed, 43 insertions(+), 62 deletions(-) diff --git a/test/estimation.jl b/test/estimation.jl index 867e277ae..50cac86d2 100644 --- a/test/estimation.jl +++ b/test/estimation.jl @@ -22,77 +22,58 @@ end @testset "jacobsen" begin fs = 100 - t = range(0, 5, step = 1/fs) + t = range(0, 5; step=1/fs) + function test_complex_jacobsen(fs, fc, f=0, t=t) + sc = cispi.(2 * fc * t .+ f) + f_est_complex = jacobsen(sc, fs) + isapprox(f_est_complex, fc; atol=1e-5) + end # test at two arbitrary frequencies - fc = -40.3 - sc = cis.(2π*fc*t .+ π/1.4) - f_est_complex = jacobsen(sc, fs) - @test isapprox(f_est_complex, fc, atol = 1e-5) - fc = 14.3 - sc = cis.(2π*fc*t .+ π/3) - f_est_complex = jacobsen(sc, fs) - @test isapprox(f_est_complex, fc, atol = 1e-5) + @test test_complex_jacobsen(fs, -40.3, 1 / 1.4) + @test test_complex_jacobsen(fs, 14.3, 1 / 3) # test near fs/2 - fc = 49.90019 - sc = cis.(2π*fc*t) - f_est_complex = jacobsen(sc, fs) - @test isapprox(f_est_complex, fc, atol = 1e-5) + @test test_complex_jacobsen(fs, 49.90019) # test near -fs/2 - fc = -49.90019 - sc = cis.(2π*fc*t) - f_est_complex = jacobsen(sc, fs) - @test isapprox(f_est_complex, fc, atol = 1e-5) + @test test_complex_jacobsen(fs, -49.90019) # test near +zero - fc = 0.04 - sc = cis.(2π*fc*t) - f_est_complex = jacobsen(sc, fs) - @test isapprox(f_est_complex, fc, atol = 1e-5) + @test test_complex_jacobsen(fs, 0.04) # test near -zero - fc = -0.1 - sc = cis.(2π*fc*t) - f_est_complex = jacobsen(sc, fs) - @test isapprox(f_est_complex, fc, atol = 1e-5) + @test test_complex_jacobsen(fs, -0.1) # tests for real signals: test only around fs/4, where the # expected error is small. fr = 28.3 - sr = cos.(2π*fr*t .+ π/4.2) + sr = cospi.(2 * fr * t .+ 1 / 4.2) f_est_real = jacobsen(sr, fs) - @test isapprox(f_est_real, fr, atol = 1e-5) + @test isapprox(f_est_real, fr; atol=1e-5) fr = 23.45 - sr = sin.(2π*fr*t .+ 3π/2.2) + sr = sinpi.(2 * fr * t .+ 3 / 2.2) f_est_real = jacobsen(sr, fs) - @test isapprox(f_est_real, fr, atol = 1e-5) + @test isapprox(f_est_real, fr; atol=1e-5) end @testset "quinn" begin ### real input fs = 100 - t = range(0, 5, step = 1/fs) + t = range(0, 5; step=1/fs) fr = 28.3 - sr = cos.(2π*fr*t .+ π/4.2) - (f_est_real, maxiter) = quinn(sr, 50, fs) - @test maxiter == false - @test isapprox(f_est_real, fr, atol = 1e-3) + sr = cospi.(2 * fr * t .+ 1 / 4.2) + function test_quinn(f, s, args...) + (f_est_real, maxiter) = quinn(s, args...) + @test maxiter == false + @test isapprox(f_est_real, f; atol=1e-3) + return nothing + end + test_quinn(fr, sr, 50, fs) # use default initial guess - (f_est_real, maxiter) = quinn(sr, fs) # initial guess given by Jacobsen - @test maxiter == false - @test isapprox(f_est_real, fr, atol = 1e-3) + test_quinn(fr, sr, fs) # initial guess given by Jacobsen # use default fs - (f_est_real, maxiter) = quinn(sr) # fs = 1.0, initial guess given by Jacobsen - @test maxiter == false - @test isapprox(f_est_real, fr/fs, atol = 1e-3) + test_quinn(fr / fs, sr) # fs = 1.0, initial guess given by Jacobsen ### complex input fc = -40.3 - sc = cis.(2π*fc*t .+ π/1.4) - (f_est_real, maxiter) = quinn(sc, -20, fs) - @test maxiter == false - @test isapprox(f_est_real, fc, atol = 1e-3) + sc = cispi.(2 * fc * t .+ 1 / 1.4) + test_quinn(fc, sc, -20, fs) # use default initial guess - (f_est_real, maxiter) = quinn(sc, fs) # initial guess given by Jacobsen - @test maxiter == false - @test isapprox(f_est_real, fc, atol = 1e-3) + test_quinn(fc, sc, fs) # initial guess given by Jacobsen # use default fs - (f_est_real, maxiter) = quinn(sc) # fs = 1.0, initial guess by Jacobsen - @test maxiter == false - @test isapprox(f_est_real, fc/fs, atol = 1e-3) + test_quinn(fc / fs, sc) # fs = 1.0, initial guess by Jacobsen end diff --git a/test/multitaper.jl b/test/multitaper.jl index 07b348398..dfe042378 100644 --- a/test/multitaper.jl +++ b/test/multitaper.jl @@ -91,8 +91,8 @@ avg_coh(x) = dropdims(mean(coherence(x); dims=3); dims=3) n_samples = 1024 t = (0:1023) ./ fs freq_range = (10, 15) - sin_1 = sin.(2π * 12.0 * t) # 12 Hz sinusoid signal - sin_2 = sin.(2π * 12.0 * t .+ π) + sin_1 = sinpi.(2 * 12.0 * t) # 12 Hz sinusoid signal + sin_2 = sinpi.(2 * 12.0 * t .+ 1) noise = rand(1024) * 2 .- 1 T = Float64 @@ -248,8 +248,8 @@ end fs = 1000.0 n_samples = 1024 t = (0:1023) ./ fs - sin_1 = sin.(2π * 12.0 * t) # 12 Hz sinusoid signal - sin_2 = sin.(2π * 12.0 * t .+ π) + sin_1 = sinpi.(2 * 12.0 * t) # 12 Hz sinusoid signal + sin_2 = sinpi.(2 * 12.0 * t .+ 1) noise = vec(readdlm(joinpath(@__DIR__, "data", "noise.txt"))) # generated by `rand(1024) * 2 .- 1` more_noisy = Array{Float64,3}(undef, 1, 2, n_samples) more_noisy[1, 1, :] = sin_1 @@ -272,8 +272,8 @@ end n_samples = 1024 t = (0:1023) ./ fs data = Array{Float64, 3}(undef, 1, 2, n_samples) - data[1, 1, :] = sin.(2π * 12.0 * t) # 12 Hz sinusoid signal - data[1, 2, :] = sin.(2π * 12.0 * t .+ π) + data[1, 1, :] = sinpi.(2 * 12.0 * t) # 12 Hz sinusoid signal + data[1, 2, :] = sinpi.(2 * 12.0 * t .+ 1) ## Script to generate reference data: # using PyMNE, DelimitedFiles # result = PyMNE.time_frequency.csd_array_multitaper(data, fs, n_fft = nextpow(2, n_samples), low_bias=true, adaptive=false) @@ -328,7 +328,7 @@ end n_samples = 1024 t = (0:1023) ./ fs noise = vec(readdlm(joinpath(@__DIR__, "data", "noise.txt"))) - signal = sin.(2π * 12.0 * t) .+ 3 * noise + signal = sinpi.(2 * 12.0 * t) .+ 3 * noise cs = mt_cross_power_spectra(reshape(signal, 1, :); fs) p = mt_pgram(signal; fs) diff --git a/test/remez_fir.jl b/test/remez_fir.jl index 9e8bfa98f..de34a779f 100644 --- a/test/remez_fir.jl +++ b/test/remez_fir.jl @@ -167,13 +167,13 @@ end Fs = 4800*L f = range(0, stop=0.5, length=10000) - P = (π*f*Fs/4800) ./ sin.(π*f*Fs/4800) - Pdb = 20*log10.(abs.(P)) + P = (π * f * Fs / 4800) ./ sinpi.(f * Fs / 4800) + Pdb = 20 * log10.(abs.(P)) Pdb[1] = 0.0 g_vec = remez(201, [ - ( 0.0, 2880.0) => (f -> (f==0) ? 1.0 : abs.((π*f/4800) ./ sin.(π*f/4800)), 1.0), - (10000.0, Fs/2) => (0.0, 100.0) + ( 0.0, 2880.0) => (f -> (f==0) ? 1.0 : abs.((π*f/4800) ./ sinpi.(f/4800)), 1.0), + (10000.0, Fs/2) => (0.0, 100.0) ]; Hz=Fs) g = PolynomialRatio(g_vec, [1.0]) Gdb = 20*log10.(abs.(freqresp(g, 2π*f))) diff --git a/test/util.jl b/test/util.jl index f613e891e..6611fc448 100644 --- a/test/util.jl +++ b/test/util.jl @@ -50,7 +50,7 @@ using Statistics: mean end @testset "fft helpers" begin - @test meanfreq(sin.(2*π*10*(0:1e-3:10*π)),1e3) ≈ 10.0 rtol=1e-3 + @test meanfreq(sinpi.(2 * 10 * (0:1e-3:10*π)), 1e3) ≈ 10.0 rtol = 1e-3 # nextfastfft @test nextfastfft(64) == 64 From e99810a7165c92fb87cd998e68316f80745de584 Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Mon, 21 Oct 2024 17:51:32 +0800 Subject: [PATCH 4/8] triple quotes restore double newline --- test/filt_stream.jl | 58 ++++++++++++++++++++++----------------------- 1 file changed, 29 insertions(+), 29 deletions(-) diff --git a/test/filt_stream.jl b/test/filt_stream.jl index 38a617a2c..ce7d086af 100644 --- a/test/filt_stream.jl +++ b/test/filt_stream.jl @@ -77,11 +77,11 @@ function test_singlerate(h::AbstractVector{T}, x::AbstractVector) where T x1 = x[1:pivotPoint] x2 = x[pivotPoint+1:end] - @printfifinteractive( "\n\n" ) - @printfifinteractive( "____ _ _ _ ____ _ ____ ____ ____ ___ ____\n" ) - @printfifinteractive( "[__ | |\\ | | __ | |___ |__/ |__| | |___\n" ) - @printfifinteractive( "___] | | \\| |__] |___ |___ | \\ | | | |___\n" ) - @printfifinteractive( "\nTesting single-rate filtering, h::%s, x::%s. xLen = %d, hLen = %d\n", typeof(h), typeof(x), xLen, hLen ) + @printfifinteractive( """\n\n + ____ _ _ _ ____ _ ____ ____ ____ ___ ____ + [__ | |\\ | | __ | |___ |__/ |__| | |___ + ___] | | \\| |__] |___ |___ | \\ | | | |___\n\n""" ) + @printfifinteractive( "Testing single-rate filtering, h::%s, x::%s. xLen = %d, hLen = %d\n", typeof(h), typeof(x), xLen, hLen ) @printfifinteractive( "\n\tfilt\n\t\t") @timeifinteractive naiveResult = filt(h, one(T), x) @@ -127,11 +127,11 @@ function test_decimation(h, x, decimation::Integer) x1 = x[1:pivotPoint] x2 = x[pivotPoint+1:end] - @printfifinteractive( "\n\n" ) - @printfifinteractive( "___ ____ ____ _ _ _ ____ ___ _ ____ _ _ \n" ) - @printfifinteractive( "| \\ |___ | | |\\/| |__| | | | | |\\ | \n" ) - @printfifinteractive( "|__/ |___ |___ | | | | | | | |__| | \\| \n" ) - @printfifinteractive( "\nTesting decimation. h::%s, x::%s. xLen = %d, hLen = %d, decimation = %d\n", typeof(h), typeof(x), xLen, hLen, decimation ) + @printfifinteractive( """\n\n + ___ ____ ____ _ _ _ ____ ___ _ ____ _ _ + | \\ |___ | | |\\/| |__| | | | | |\\ | + |__/ |___ |___ | | | | | | | |__| | \\|\n\n""" ) + @printfifinteractive( "Testing decimation. h::%s, x::%s. xLen = %d, hLen = %d, decimation = %d\n", typeof(h), typeof(x), xLen, hLen, decimation ) @printfifinteractive( "\n\tNaive decimation\n\t\t") @timeifinteractive naiveResult = naivefilt(h, x, 1//decimation) @@ -178,11 +178,11 @@ function test_interpolation(h::AbstractVector{T}, x::AbstractVector{V}, interpol x1 = x[1:pivotPoint] x2 = x[pivotPoint+1:end] - @printfifinteractive( "\n\n" ) - @printfifinteractive( "_ _ _ ___ ____ ____ ___ _ ____ ____ ___ _ ____ _ _ \n" ) - @printfifinteractive( "| |\\ | | |___ |__/ |__] | | | |__| | | | | |\\ | \n" ) - @printfifinteractive( "| | \\| | |___ | \\ | |___ |__| | | | | |__| | \\| \n" ) - @printfifinteractive( "\nTesting interpolation, h::%s, x::%s. xLen = %d, hLen = %d, interpolation = %d\n", typeof(h), typeof(x), xLen, hLen, interpolation ) + @printfifinteractive( """\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") @timeifinteractive begin @@ -238,15 +238,15 @@ function test_rational(h, x, ratio) downfactor = denominator(ratio) resultType = promote_type(eltype(h), eltype(x)) - @printfifinteractive( "\n\n" ) - @printfifinteractive( " ____ ____ ___ _ ____ _ _ ____ _ \n" ) - @printfifinteractive( " |__/ |__| | | | | |\\ | |__| | \n" ) - @printfifinteractive( " | \\ | | | | |__| | \\| | | |___ \n" ) - @printfifinteractive( " \n" ) - @printfifinteractive( "____ ____ ____ ____ _ _ ___ _ _ _ _ ____\n" ) - @printfifinteractive( "|__/ |___ [__ |__| |\\/| |__] | | |\\ | | __\n" ) - @printfifinteractive( "| \\ |___ ___] | | | | | |___ | | \\| |__]\n" ) - @printfifinteractive( "\n\nTesting rational resampling, h::%s, x::%s. xLen = %d, hLen = %d, ratio = %d//%d\n", typeof(h), typeof(x), xLen, hLen, upfactor, downfactor ) + @printfifinteractive( """\n\n + ____ ____ ___ _ ____ _ _ ____ _ + |__/ |__| | | | | |\\ | |__| | + | \\ | | | | |__| | \\| | | |___ + + ____ ____ ____ ____ _ _ ___ _ _ _ _ ____ + |__/ |___ [__ |__| |\\/| |__] | | |\\ | | __ + | \\ |___ ___] | | | | | |___ | | \\| |__]\n\n""") + @printfifinteractive( "Testing rational resampling, h::%s, x::%s. xLen = %d, hLen = %d, ratio = %d//%d\n", typeof(h), typeof(x), xLen, hLen, upfactor, downfactor ) @printfifinteractive( "\n\tNaive rational resampling\n\t\t") @timeifinteractive naiveResult = naivefilt(h, x, ratio) @@ -293,11 +293,11 @@ function test_arbitrary(Th, x, resampleRate, numFilters::Integer) myfilt = FIRFilter(h, resampleRate, numFilters) xLen = length(x) - @printfifinteractive("\n\n") - @printfifinteractive( "____ ____ ___ ____ ____ ____ ____ _ _ ___ _ _ _ _ ____\n" ) - @printfifinteractive( "|__| |__/ |__] |__/ |___ [__ |__| |\\/| |__] | | |\\ | | __\n" ) - @printfifinteractive( "| | | \\ |__] . | \\ |___ ___] | | | | | |___ | | \\| |__]\n" ) - @printfifinteractive( "\nh::%s, x::%s, rate = %f, Nϕ = %d, xLen = %d\n", typeof(h), typeof(x), resampleRate, numFilters, xLen ) + @printfifinteractive("""\n\n + ____ ____ ___ ____ ____ ____ ____ _ _ ___ _ _ _ _ ____ + |__| |__/ |__] |__/ |___ [__ |__| |\\/| |__] | | |\\ | | __ + | | | \\ |__] . | \\ |___ ___] | | | | | |___ | | \\| |__]\n\n""") + @printfifinteractive( "h::%s, x::%s, rate = %f, Nϕ = %d, xLen = %d\n", typeof(h), typeof(x), resampleRate, numFilters, xLen ) @printfifinteractive( "\n\tNaive arbitrary resampling\n\t\t" ) @timeifinteractive naiveResult = naivefilt(h, x, resampleRate, numFilters) From 46037642d68e0683f1aca5c06196f567e187c8cf Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Sat, 26 Oct 2024 12:11:27 +0800 Subject: [PATCH 5/8] test_mtspec, atol=eps --- test/FilterTestHelpers.jl | 37 +++++++++++++------------------------ test/periodograms.jl | 39 +++++++++++++++------------------------ 2 files changed, 28 insertions(+), 48 deletions(-) diff --git a/test/FilterTestHelpers.jl b/test/FilterTestHelpers.jl index 81bda3906..3041be557 100644 --- a/test/FilterTestHelpers.jl +++ b/test/FilterTestHelpers.jl @@ -27,12 +27,12 @@ function zpkfilter_eq(f1, f2) @test map(Float64, f1.k) ≈ map(Float64, f2.k) end -function zpkfilter_eq(f1, f2, eps) +function zpkfilter_eq(f1, f2, atol) if !isempty(f1.z) || !isempty(f2.z) - @test ≈(map(ComplexF64, sort(f1.z; lt)), map(ComplexF64, sort(f2.z; lt)); atol=eps) + @test ≈(map(ComplexF64, sort(f1.z; lt)), map(ComplexF64, sort(f2.z; lt)); atol) end - @test ≈(map(ComplexF64, sort(f1.p; lt)), map(ComplexF64, sort(f2.p; lt)); atol=eps) - @test ≈(map(Float64, f1.k), map(Float64, f2.k); atol=eps) + @test ≈(map(ComplexF64, sort(f1.p; lt)), map(ComplexF64, sort(f2.p; lt)); atol) + @test ≈(map(Float64, f1.k), map(Float64, f2.k); atol) end loss(x::Real, y::Real) = abs(float(x) - float(y))/eps(float(x)) @@ -58,33 +58,22 @@ function tffilter_accuracy(f1, f2, accurate_f) accuracy_check(loss(a1, accurate_a), loss(a2, accurate_a), "a") end -function zpkfilter_accuracy(f1, f2, accurate_f; relerr=1, compare_gain_at=nothing, eps=nothing) +function zpkfilter_accuracy(f1, f2, accurate_f; relerr=1, compare_gain_at=nothing, eps=0.0) z1, p1 = sort(f1.z; lt), sort(f1.p; lt) z2, p2 = sort(f2.z; lt), sort(f2.p; lt) accurate_z, accurate_p = sort(accurate_f.z; lt), sort(accurate_f.p; lt) + atol = eps if !isempty(z1) || !isempty(z2) || !isempty(accurate_z) - if eps !== nothing - @test ≈(z1, accurate_z; atol=eps) - @test ≈(z2, accurate_z; atol=eps) - else - @test z1 ≈ accurate_z - @test z2 ≈ accurate_z - end + @test ≈(z1, accurate_z; atol) + @test ≈(z2, accurate_z; atol) accuracy_check(loss(z1, accurate_z), loss(z2, accurate_z), "z", relerr) end - if eps !== nothing - @test ≈(p1, accurate_p; atol=eps) - @test ≈(p2, accurate_p; atol=eps) - @test ≈(f1.k, accurate_f.k; atol=eps) - @test ≈(f2.k, accurate_f.k; atol=eps) - else - @test p1 ≈ accurate_p - @test p2 ≈ accurate_p - @test f1.k ≈ accurate_f.k - @test f2.k ≈ accurate_f.k - end + @test ≈(p1, accurate_p; atol) + @test ≈(p2, accurate_p; atol) + @test ≈(f1.k, accurate_f.k; atol) + @test ≈(f2.k, accurate_f.k; atol) accuracy_check(loss(p1, accurate_p), loss(p2, accurate_p), "p", relerr) - if compare_gain_at === nothing + if isnothing(compare_gain_at) accuracy_check(loss(f1.k, accurate_f.k), loss(f2.k, accurate_f.k), "k", relerr) else jω = compare_gain_at*im diff --git a/test/periodograms.jl b/test/periodograms.jl index d25d93767..9db260eae 100644 --- a/test/periodograms.jl +++ b/test/periodograms.jl @@ -37,58 +37,49 @@ using FFTW: fftfreq @test time(mt_spec) == t @test power(mt_spec)[:,1] ≈ power(mt_pgram(x0[1:256]; fs=10)) + function test_mtspec(m, mt_spec) + @test power(m) ≈ power(mt_spec) + @test freq(m) == freq(mt_spec) + @test time(m) == time(mt_spec) + nothing + end + # in-place full precision version: spec_config = MTSpectrogramConfig{Float64}(length(x0), 256, 128; fs=10) out = allocate_output(spec_config) mt_spec2 = mt_spectrogram!(out, x0, spec_config) - @test power(mt_spec2) ≈ power(mt_spec) - @test freq(mt_spec2) == freq(mt_spec) - @test time(mt_spec2) == time(mt_spec) + test_mtspec(mt_spec2, mt_spec) # out-of-place with config: r = mt_spectrogram(x0, spec_config) - @test power(r) ≈ power(mt_spec) - @test freq(r) == freq(mt_spec) - @test time(r) == time(mt_spec) + test_mtspec(r, mt_spec) # in-place without config: r = mt_spectrogram!(out, x0, 256, 128; fs=10) - @test power(r) ≈ power(mt_spec) - @test freq(r) == freq(mt_spec) - @test time(r) == time(mt_spec) + test_mtspec(r, mt_spec) # in-place half precision version: spec_config = MTSpectrogramConfig{Float32}(length(x0), 256, 128; fs=10) out = allocate_output(spec_config) mt_spec3 = mt_spectrogram!(out, x0, spec_config) - @test power(mt_spec3) ≈ power(mt_spec) - @test freq(mt_spec3) == freq(mt_spec) - @test time(mt_spec3) == time(mt_spec) + test_mtspec(mt_spec3, mt_spec) mt_spec4 = mt_spectrogram!(out, Float32.(x0), spec_config) - @test power(mt_spec4) ≈ power(mt_spec) - @test freq(mt_spec4) == freq(mt_spec) - @test time(mt_spec4) == time(mt_spec) + test_mtspec(mt_spec4, mt_spec) # We can also only pass the window config. Full precision: config = MTConfig{Float64}(256; fs=10) mt_spec5 = mt_spectrogram(x0, config, 128) - @test power(mt_spec5) ≈ power(mt_spec) - @test freq(mt_spec5) == freq(mt_spec) - @test time(mt_spec5) == time(mt_spec) + test_mtspec(mt_spec5, mt_spec) # We can also only pass the window config. Half precision: config = MTConfig{Float32}(256; fs=10) mt_spec6 = mt_spectrogram(x0, config, 128) - @test power(mt_spec6) ≈ power(mt_spec) - @test freq(mt_spec6) == freq(mt_spec) - @test time(mt_spec6) == time(mt_spec) + test_mtspec(mt_spec6, mt_spec) # with Float32 input: mt_spec7 = mt_spectrogram(Float32.(x0), config, 128) - @test power(mt_spec7) ≈ power(mt_spec) - @test freq(mt_spec7) == freq(mt_spec) - @test time(mt_spec7) == time(mt_spec) + test_mtspec(mt_spec7, mt_spec) @test_throws DimensionMismatch mt_spectrogram!(similar(out, size(out, 1), size(out, 2)+1), x0, spec_config) @test_throws DimensionMismatch mt_spectrogram!(out, vcat(x0, x0), spec_config) From 9a6c8480196d6d5b8ce2317df9183f407f0d48de Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Sun, 27 Oct 2024 14:52:38 +0800 Subject: [PATCH 6/8] remove more DSP. --- test/filt.jl | 6 +++--- test/filter_conversion.jl | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/test/filt.jl b/test/filt.jl index 570a4ec42..77239272f 100644 --- a/test/filt.jl +++ b/test/filt.jl @@ -293,13 +293,13 @@ end @testset "filtfilt SOS" begin x = readdlm(joinpath(dirname(@__FILE__), "data", "spectrogram_x.txt"),'\t') - f = DSP.digitalfilter(DSP.Lowpass(0.2), DSP.Butterworth(4)) + f = digitalfilter(Lowpass(0.2), Butterworth(4)) @test filtfilt(convert(SecondOrderSections, f), x) ≈ filtfilt(convert(PolynomialRatio, f), x) - f = DSP.digitalfilter(DSP.Highpass(0.1), DSP.Butterworth(6)) + f = digitalfilter(Highpass(0.1), Butterworth(6)) @test filtfilt(convert(SecondOrderSections, f), x) ≈ filtfilt(convert(PolynomialRatio, f), x) - f = DSP.digitalfilter(DSP.Bandpass(0.1, 0.3), DSP.Butterworth(2)) + f = digitalfilter(Bandpass(0.1, 0.3), Butterworth(2)) @test filtfilt(convert(SecondOrderSections, f), x) ≈ filtfilt(convert(PolynomialRatio, f), x) end diff --git a/test/filter_conversion.jl b/test/filter_conversion.jl index 85f3d8a88..3471c6b73 100644 --- a/test/filter_conversion.jl +++ b/test/filter_conversion.jl @@ -133,7 +133,7 @@ using Polynomials.PolyCompat 1.0000000000000000e+00 -1.9021224191804869e+00 1.0000000000000000e+00 1.0000000000000000e+00 -1.8964983429993663e+00 9.9553672990017417e-01 1.0000000000000000e+00 -1.9021224191804869e+00 1.0000000000000000e+00 1.0000000000000000e+00 -1.8992956433548462e+00 9.9559721515078736e-01 ] - f = convert(SecondOrderSections, digitalfilter(Bandstop(49.5, 50.5), DSP.Butterworth(2); fs=1000)) + f = convert(SecondOrderSections, digitalfilter(Bandstop(49.5, 50.5), Butterworth(2); fs=1000)) @test m_sos_butterworth_bs ≈ sosfilter_to_matrix(f) @test f.g ≈ 0.995566972017647 From c4787f28ba609a11f9492b7913afc445ad45640c Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Thu, 12 Dec 2024 13:52:16 +0800 Subject: [PATCH 7/8] diric DomainError --- src/diric.jl | 2 +- test/diric.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/diric.jl b/src/diric.jl index 6580ed8e0..2b760a4d4 100644 --- a/src/diric.jl +++ b/src/diric.jl @@ -36,7 +36,7 @@ julia> diric(0, 4) ``` """ function diric(Ω::T, n::Integer) where T <: AbstractFloat - n > 0 || throw(ArgumentError("n=$n not positive")) + n > 0 || throw(DomainError(n, "n not positive")) sign = one(T) if isodd(n) Ω = rem2pi(Ω, RoundNearest) # [-π,π) diff --git a/test/diric.jl b/test/diric.jl index d6dc5b6b9..6d3efaaaa 100644 --- a/test/diric.jl +++ b/test/diric.jl @@ -1,7 +1,7 @@ # diric.jl Dirichlet kernel tests @testset "diric" begin - @test_throws ArgumentError diric(0, -2) + @test_throws DomainError diric(0, -2) @test @inferred(diric(0, 4)) ≈ 1 @test @inferred(diric(0 // 1, 5)) == 1 @test @inferred(diric(4π, 4)) ≈ 1 From ef474bbf00fb954853375ea9fb95394fb26dedcc Mon Sep 17 00:00:00 2001 From: wheeheee <104880306+wheeheee@users.noreply.github.com> Date: Fri, 13 Dec 2024 18:16:27 +0800 Subject: [PATCH 8/8] Apply suggestions from code review - foreach - move test_quinn - eps -> atol --- test/FilterTestHelpers.jl | 3 +-- test/estimation.jl | 11 ++++++----- test/filt_stream.jl | 2 +- test/filter_design.jl | 4 ++-- 4 files changed, 10 insertions(+), 10 deletions(-) diff --git a/test/FilterTestHelpers.jl b/test/FilterTestHelpers.jl index 3041be557..dbf908149 100644 --- a/test/FilterTestHelpers.jl +++ b/test/FilterTestHelpers.jl @@ -58,11 +58,10 @@ function tffilter_accuracy(f1, f2, accurate_f) accuracy_check(loss(a1, accurate_a), loss(a2, accurate_a), "a") end -function zpkfilter_accuracy(f1, f2, accurate_f; relerr=1, compare_gain_at=nothing, eps=0.0) +function zpkfilter_accuracy(f1, f2, accurate_f; relerr=1, compare_gain_at=nothing, atol=0.0) z1, p1 = sort(f1.z; lt), sort(f1.p; lt) z2, p2 = sort(f2.z; lt), sort(f2.p; lt) accurate_z, accurate_p = sort(accurate_f.z; lt), sort(accurate_f.p; lt) - atol = eps if !isempty(z1) || !isempty(z2) || !isempty(accurate_z) @test ≈(z1, accurate_z; atol) @test ≈(z2, accurate_z; atol) diff --git a/test/estimation.jl b/test/estimation.jl index 50cac86d2..b68dbe870 100644 --- a/test/estimation.jl +++ b/test/estimation.jl @@ -52,22 +52,23 @@ end end @testset "quinn" begin - ### real input - fs = 100 - t = range(0, 5; step=1/fs) - fr = 28.3 - sr = cospi.(2 * fr * t .+ 1 / 4.2) function test_quinn(f, s, args...) (f_est_real, maxiter) = quinn(s, args...) @test maxiter == false @test isapprox(f_est_real, f; atol=1e-3) return nothing end + ### real input + fs = 100 + t = range(0, 5; step=1/fs) + fr = 28.3 + sr = cospi.(2 * fr * t .+ 1 / 4.2) test_quinn(fr, sr, 50, fs) # use default initial guess test_quinn(fr, sr, fs) # initial guess given by Jacobsen # use default fs test_quinn(fr / fs, sr) # fs = 1.0, initial guess given by Jacobsen + ### complex input fc = -40.3 sc = cispi.(2 * fc * t .+ 1 / 1.4) diff --git a/test/filt_stream.jl b/test/filt_stream.jl index ce7d086af..f54470f60 100644 --- a/test/filt_stream.jl +++ b/test/filt_stream.jl @@ -324,7 +324,7 @@ function test_arbitrary(Th, x, resampleRate, numFilters::Integer) results = (naiveResult, statelessResult, statefulResult, piecewiseResult) commonLen = minimum(length, results) - resize!.(results, commonLen) + foreach(x -> resize!(x, commonLen), results) @test naiveResult ≈ statelessResult @test naiveResult ≈ statefulResult diff --git a/test/filter_design.jl b/test/filter_design.jl index e07d697c6..ed10d7102 100644 --- a/test/filter_design.jl +++ b/test/filter_design.jl @@ -266,7 +266,7 @@ end f = Elliptic(20, 0.1, 10) zpkfilter_eq(f, matlab_f) - zpkfilter_accuracy(f, matlab_f, Elliptic(BigFloat, 20, 0.1, 10), eps=1e-9) + zpkfilter_accuracy(f, matlab_f, Elliptic(BigFloat, 20, 0.1, 10), atol=1e-9) # 19 pole elliptic filter prototype with 0.1 dB passband ripple and 10 # dB stopband ripple from MATLAB 2013b: @@ -309,7 +309,7 @@ end f = Elliptic(19, 0.1, 10) zpkfilter_eq(f, matlab_f) - zpkfilter_accuracy(f, matlab_f, Elliptic(BigFloat, 19, 0.1, 10), eps=4e-9) + zpkfilter_accuracy(f, matlab_f, Elliptic(BigFloat, 19, 0.1, 10), atol=4e-9) end #