How to apply a butterworth filter in a loop?

Specific Domains Signal and Image Processing
, ,
1

I constructed a Butterworth filter, which works fine:

using ControlPlots, DSP

# Define the cut-off frequency in Hz
cut_off_freq = 2.0
Wn = cut_off_freq

# Define the sampling time in seconds
dt = 0.05
fs = 1/dt
sim_time = 4.0
N = Int(sim_time / dt)

# Design the filter
butter = Filters.digitalfilter(Filters.Lowpass(Wn; fs), Filters.Butterworth(4))

# Create an array of measurements
measurements = zeros(N)
for i in Int(N/2):N
    measurements[i] = 1.0
end

results = filt(butter, measurements)

# Plot the results
plot((1:N)*dt, [measurements, results]; xlabel="Time (s)", ylabel="Amplitude", 
     fig="Forth order Butterworth Filter")

image
image624×504 17.9 KB

But now I want to apply it online in a loop, similar to the way I am applying a first order filter in the following example:

results = zeros(N)

# Apply the EMA filter
last_measurement::Float64 = 0.0
for i in 1:N
    global last_measurement
    results[i] = ema_filter(measurements[i], last_measurement, cut_off_freq, dt)
    last_measurement = results[i]
end

function ema_filter(measurement, last_measurement, cut_off_freq, dt)
    if cut_off_freq > 0.0
        alpha = dt / (dt + one(measurement) / (2π * cut_off_freq))
        filtered_value = alpha * measurement + (one(measurement) - alpha) * last_measurement
    else
        filtered_value = measurement
    end
    return filtered_value
end

The point is, the filter will become part of a larger control system.

How can I apply the Butterworth filter sample-by-sample?

2

Can’t you just call DSP.filt! on one sample (or, more efficiently, on a batch of a few samples) at a time, passing the si array to track the cumulative effect of the previous samples?

More simply, you can use a stateful filter wrapper object with something like:

butterF = DSP.Filters.DF2TFilter(butter)

and then call filt!(out, butterF, x) to perform the filtering on a chunk of inputs x to get a chunk of outputs out.

3 Likes
3

This works, but is very inefficient:

using ControlPlots, DSP

# Define the cut-off frequency in Hz
cut_off_freq = 2.0

# Define the sampling time in seconds
dt = 0.05
fs = 1/dt
sim_time = 4.0
N = Int(sim_time / dt)

# Design the filter
butter = Filters.digitalfilter(Filters.Lowpass(cut_off_freq; fs), Filters.Butterworth(2))

# Create an array of measurements
measurements = zeros(N)
for i in Int(N/2):N
    measurements[i] = 1.0
end

function butter_filter(butter, measurement, buffer, index)
    buffer[index] = measurement
    res = filt(butter, buffer[1:index])
    return res[index]
end

buffer = zeros(N)
results = zeros(N)
for i in 1:N
    results[i] = butter_filter(butter, measurements[i], buffer, i)
end

# Plot the results
plot((1:N)*dt, [measurements, results]; xlabel="Time (s)", ylabel="Amplitude", 
     fig="Forth order Butterworth Filter")

I was reading I should convert the ZPG notation to transfer function coefficients and use those to create a difference equation which would be short and efficient. No idea how to do that, though.

4

Of course it’s inefficient — you are re-doing the entire filter from the beginning for every new measurement.

Using the filter state (see my suggestion) allows you to avoid this — each batch of new measurements requires only the new filter terms to be applied.

5

NIce, but I do not know how to put your suggestion into practice. A big issue with the DSP package is that there are no examples provided.

6

What did you try? The documentation may lack examples written out for you, but it documents the relevant functions.

For example, this code applies your filter 1 element at a time:

butterF = DSP.Filters.DF2TFilter(butter)
results2 = zeros(N)
for i = 1:N
    @views filt!(results2[i:i], butterF, measurements[i:i])
end

(Bigger batches should be more efficient, but the complexity should still be linear in N.)

3 Likes
7

I am using now this function:

"""
    apply_filter(tfilter::DF2TFilter, measurement)

Apply the filter to the measurement.

# Arguments
- `tfilter`: The filter, created with `create_filter` and converted to a DF2TFilter
- `measurement`: The measurement value
"""
function apply_filter(tfilter::DF2TFilter, measurement)
    results = zeros(1)
    measurements = ones(1) * measurement
    @views filt!(results[1:1], tfilter, measurements)
    return results[1]
end

While it is fast enough for my purposes, I think it is ugly and a function that does not require to allocate one element arrays should be added to DSP.jl.

I created an issue: Improve online filter functions · Issue #605 · JuliaDSP/DSP.jl · GitHub