# Not sure how to apply a filter

I’ve constructed a filter as follows:

``````filter1 = let
f1 = 6.31
f2 = 1258.9
Q1 = 0.71
f3 = 15.915
f4 = 15.915
Q2 = 0.64
SecondOrderSections([
Biquad(0.0, 2π*f4^2/f3, 4π^2*f4^2, 2π*f4/Q2, 4π^2*f4^2)], 1.0)
end
``````

and I’ve plotted the frequency response:

`````` let
fr=1.0:1.0:1200.0
plot(fr, abs.(freqs(filter1, fr.*2π)),
yaxis=:log,xaxis=:log, title = "filter response", leg=false)
end
``````

assuming `using DSP, Plots`.

I would like to apply this filter to sampled data. I have successfully created digital filters following the examples in the docs, starting with `responsetype`s and `designmethod`s, but I’m not sure how to create a `digitalfilter` from a `SecondOrderSections`.

You can just use `filt(filter1, sampled_data)`, I believe.

That works for a digital filter where you have designed it for a given sample rate, but I don’t know to use an analog filter. The docs say that it will “Apply the filter … along the first dimension of array x”, but nowhere have I set the sample rate and when I tried anyway just to see what I would get, I got extreme values and NaNs. I think there’s a distinction between analog filters and digital filters in DSP.jl that I don’t understand.

The distinction between analog and digital filter in this context is better expressed as continuous-time (Laplace transform) or discrete-time (Z-Transform). If you intend to apply this filter to a sampled signal, you should be operating in (or convert your coefficients to) a discrete-time representation. In the discrete-time case, frequencies are represented in radians/sample, so the values given for `fn` in your filter coefficients and response plotting are suspiciously large. I suspect you meant to divide each of those by your sampling frequency `Fs`.

1 Like

I’ve been reading through the docs, some of the source, and some Matlab docs for `sosfilt`, and I think you are right. The cutoff frequencies I used above are indeed in Hz (this is a filter specified in a standard by the polynomials, and the above definition gives the right bode plot), where the Matlab example defines them in normalized frequency in units of π rd / sample or half-cycles per sample. I’m at the wrong computer right now but in the morning will try multiplying the `fn` values by `2/fsample`, and then applying the filter.

Thanks!

I still don’t know what I’m doing. I have a alternate solution using a digital filter constructed differently but I would still like to know what I’m doing wrong.

Here’s my example:

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

using Markdown
using InteractiveUtils

# ╔═╡ 84c78010-284a-11eb-1385-cd8bd55f042d
begin
using DSP
using Plots
end

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

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

# ╔═╡ 0938a0a0-28e5-11eb-0e2f-9d28fb0ab8b9
let
f=10.0
tr = 0.0:Tsample:1.0
x=sin.(2π*f.*tr)
p1=plot(tr, x, label="input")
fx = filt(filter1, x)
p2=plot(tr,fx, label="output")
end

# ╔═╡ d6b631f0-2863-11eb-1f59-bb236634efcc
let
fr=1.0:1.0:1200.0
plot(fr, abs.(freqs(filter1, fr.*2π*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
# ╠═84c78010-284a-11eb-1385-cd8bd55f042d
``````

The bode plot that is produced is correct, but the output waveform in the plot has a magnitude of about 0.004, regardless of if I set `f` to 1.0, 10.0, 100.0 or 1000.0 Hz.

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())
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}(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)
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

``````

Bilinear transform, bilinear. Now that’s a term I’ve not heard in a long, long time. A long time.

I did study this in one of my EE classes, but back where there were only three Star Wars movies.

Thank you, it’s making sense now.