Not sure how to apply a filter

Specific Domains Signal and Image Processing
7

Ah, I think I see now. The fact that freqs gives the correct plot (as opposed to freqz) means that the Biquad coefficients you have are in the Lapace domain (continuous time, usually represented as a polynomial in s). filt assumes it’s input is in the Z domain (discrete time, usually represented as a polynomial in z). The Bilinear Transform can convert coefficients between the two representations.

The (unreleased) latest version of DSP.jl has tools to make this difference explicit, try this out

### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# ╔═╡ 84c78010-284a-11eb-1385-cd8bd55f042d
begin
        using Pkg
        Pkg.activate(mktempdir())
        Pkg.add(PackageSpec(name="DSP", rev="master"))
        Pkg.add("Plots")
end

# ╔═╡ 958dba22-2949-11eb-15b6-7f836480cc87
begin
        using DSP
        using Plots
end

# ╔═╡ 620b8630-2871-11eb-06f5-c7d982acd7bd
Tsample=50e-6

# ╔═╡ 83664af0-2861-11eb-084e-6de55511c1dc
filter_s = let
        f1 = 6.31*Tsample*2
        f2 = 1258.9*Tsample*2
        Q1 = 0.71
        f3 = 15.915*Tsample*2
        f4 = 15.915*Tsample*2
        Q2 = 0.64
        SecondOrderSections([
                        Biquad{:s}(4π^2*f2^2,0.0,0.0,2π*f1/Q1, 4π^2*f1^2),
                        Biquad{:s}(0.0,0.0,1.0,2π*f2/Q1, 4π^2*f2^2),
                        Biquad{:s}(0.0, 2π*f4^2/f3, 4π^2*f4^2, 2π*f4/Q2, 4π^2*f4^2)], 1.0)
end

# ╔═╡ b504b86a-2958-11eb-28b5-0df48a71b1b2
 # since frequencies are normalized already, fs is just 2
filter1 = DSP.Filters.bilinear(filter_s, 2)

# ╔═╡ 0938a0a0-28e5-11eb-0e2f-9d28fb0ab8b9
let
        f=10.0
        tr = 0.0:Tsample:10.0
        # x=sin.(2π*(f/4)*tr) + sin.(2π*f*tr) + sin.(2π*4*f*tr)
        x = rand(length(tr))
        perin = periodogram(x; fs=1. / Tsample)
        imax=10000
        imin=10
        p1=plot(freq(perin)[imin:imax], power(perin)[imin:imax], yaxis=:log, xaxis=:log, label="input", leg=false)
        fx = filt(filter1, x)
        perout = periodogram(fx; fs=1. / Tsample)
        p2=plot(freq(perout)[imin:imax], power(perout)[imin:imax], yaxis=:log, xaxis=:log, label="output", leg=false)
        plot(p1, p2, layout=(2,1), link=:x)
end

# ╔═╡ d6b631f0-2863-11eb-1f59-bb236634efcc
let
        fr=1.0:1.0:1200.0
        plot(fr, abs.(freqz(filter1, fr, 1. / Tsample)),
                yaxis=:log,xaxis=:log, title = "filter response", leg=false)
end

# ╔═╡ Cell order:
# ╠═0938a0a0-28e5-11eb-0e2f-9d28fb0ab8b9
# ╠═d6b631f0-2863-11eb-1f59-bb236634efcc
# ╟─620b8630-2871-11eb-06f5-c7d982acd7bd
# ╠═83664af0-2861-11eb-084e-6de55511c1dc
# ╠═b504b86a-2958-11eb-28b5-0df48a71b1b2
# ╠═84c78010-284a-11eb-1385-cd8bd55f042d
# ╠═958dba22-2949-11eb-15b6-7f836480cc87