Not sure how to apply a filter

klaff · November 16, 2020, 11:56pm

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(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:

 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 responsetypes and designmethods, but I’m not sure how to create a digitalfilter from a SecondOrderSections.

mbaz · November 17, 2020, 2:01am

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

klaff · November 17, 2020, 5:06am

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.

contradict · November 17, 2020, 6:26am

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.

klaff · November 17, 2020, 6:50am

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!

klaff · November 17, 2020, 11:04pm

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

contradict · November 18, 2020, 5:36am

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

klaff · November 18, 2020, 4:11pm

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.