Skip to content

Commit

Permalink
new, faster version
Browse files Browse the repository at this point in the history
  • Loading branch information
Uwe Fechner committed Dec 4, 2024
1 parent e7a0ed8 commit 14e85f4
Show file tree
Hide file tree
Showing 2 changed files with 51 additions and 0 deletions.
1 change: 1 addition & 0 deletions examples/test_butter.jl
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ results = zeros(N)
for i in 1:N
results[i] = apply_filter(butter, measurements[i], buffer, i)
end
@time apply_filter(butter, measurements[N], buffer, N)

# Plot the step response
p = plot((1:N)*dt, [measurements, results]; xlabel="Time (s)", ylabel="Amplitude",
Expand Down
50 changes: 50 additions & 0 deletions examples/test_butter2.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,50 @@
using DiscreteFilters, ControlPlots, DSP, ControlSystemsBase, Printf, LaTeXStrings
include("plotting.jl")

function create_filter2(cut_off_freq; order=4, type=:Butter, dt)
if type == :Butter
return (Filters.digitalfilter(Filters.Lowpass(cut_off_freq; fs=1/dt), Filters.Butterworth(order)))
elseif type == :Cheby1
return (Filters.digitalfilter(Filters.Lowpass(cut_off_freq; fs=1/dt), Filters.Chebyshev1(order, 0.01)))
end
end
function apply_filter2(butterF, measurement, index)
results = zeros(1)
measurements = ones(1) * measurement
@views filt!(results[1:1], butterF, measurements)
return results[1]
end

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

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

# Design the filter
butter = create_filter2(cut_off_freq; order=4, dt)

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

# apply the filter
results = zeros(N)
tfilter = DSP.Filters.DF2TFilter(butter)
for i in 1:N
results[i] = apply_filter2(tfilter, measurements[i], i)
end
@time apply_filter2(tfilter, measurements[N], N)

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

# Plot the frequency response
bo = tf(butter, dt)
bode_plot(bo; hz=true, from=0.5, to=1.5, title="4th order Butterworth Filter")

0 comments on commit 14e85f4

Please sign in to comment.