# Not sure how to apply a filter

**URL:** <https://discourse.julialang.org/t/not-sure-how-to-apply-a-filter/50265>\
**Category:** Signal and Image Processing\
**Tags:** dsp\
**Created:** [November 16, 2020, 11:56pm UTC](https://discourse.julialang.org/t/not-sure-how-to-apply-a-filter/50265 "2020-11-16T23:56:29Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [November 16, 2020, 11:56pm UTC](https://discourse.julialang.org/t/not-sure-how-to-apply-a-filter/50265/1 "2020-11-16T23:56:29Z")

</div>

I’ve constructed a filter as follows:

```julia
filter1 = let
	f1 = 6.31
	f2 = 1258.9
	Q1 = 0.71
	f3 = 15.915
	f4 = 15.915
	Q2 = 0.64
	SecondOrderSections([
			Biquad(4π^2*f2^2,0.0,0.0,2π*f1/Q1, 4π^2*f1^2),
			Biquad(0.0,0.0,1.0,2π*f2/Q1, 4π^2*f2^2),
			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:

```julia
 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`.

---

<div class="post-metadata">

**Author:** ![mbaz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbaz/32/17295_2.png) [@mbaz](https://discourse.julialang.org/u/mbaz)\
**Post date:** [November 17, 2020, 2:01am UTC](https://discourse.julialang.org/t/not-sure-how-to-apply-a-filter/50265/2 "2020-11-17T02:01:35Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [November 17, 2020, 5:06am UTC](https://discourse.julialang.org/t/not-sure-how-to-apply-a-filter/50265/3 "2020-11-17T05:06:40Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [November 17, 2020, 6:26am UTC](https://discourse.julialang.org/t/not-sure-how-to-apply-a-filter/50265/4 "2020-11-17T06:26:15Z")

</div>

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`.

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [November 17, 2020, 6:50am UTC](https://discourse.julialang.org/t/not-sure-how-to-apply-a-filter/50265/5 "2020-11-17T06:50:08Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [November 17, 2020, 11:04pm UTC](https://discourse.julialang.org/t/not-sure-how-to-apply-a-filter/50265/6 "2020-11-17T23:04:31Z")

</div>

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:

```julia
### 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(4π^2*f2^2,0.0,0.0,2π*f1/Q1, 4π^2*f1^2),
			Biquad(0.0,0.0,1.0,2π*f2/Q1, 4π^2*f2^2),
			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")
	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.(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.

---

<div class="post-metadata">

**Author:** ![contradict](https://avatars.discourse-cdn.com/v4/letter/c/ac91a4/32.png) [@contradict](https://discourse.julialang.org/u/contradict)\
**Post date:** [November 18, 2020, 5:36am UTC](https://discourse.julialang.org/t/not-sure-how-to-apply-a-filter/50265/7 "2020-11-18T05:36:29Z")

</div>

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](https://en.wikipedia.org/wiki/Bilinear_transform#Transformation_of_a_second-order_biquad) can convert coefficients between the two representations.

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

```julia
### 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

```

---

<div class="post-metadata">

**Author:** ![klaff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/klaff/32/7637_2.png) [@klaff](https://discourse.julialang.org/u/klaff)\
**Post date:** [November 18, 2020, 4:11pm UTC](https://discourse.julialang.org/t/not-sure-how-to-apply-a-filter/50265/8 "2020-11-18T16:11:53Z")

</div>

_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.