How can I implement a filtered, discrete time derivative following the equations:
u_1 = u − y_1
dy_1/dt = (1 / T_D) * u_1
y = (K_D / T_D) * u_1
filtertau = 0.01
as described on page 128 of the following document: https://home.engineering.iastate.edu/~jdm/ee554/IEEEstd421.5-2016RecPracExSysModsPwrSysStabStudies.pdf
?
Can you clarify for us non-EEs what you mean by implement?
Like, physically with a soldering iron? Or for a simulation in code? But what do you want there?
To me it looks like you have an input current u(t), which linearly drives an internal state \partial_ty_1(t) = -(1/TD)* (y_1(t) - u(t)). From that you obtain an output x(t) = F(u(t), y_1(t)) with a nonlinear function F.
You don’t need a differential equation solver for that, it’s mere cubature via variation of constants, y_1(t_2) = e^{-(t_2-t_1)/TD}y_1(t_1) + \frac 1 {TD} \int_{t_1}^{t_2}e^{-(t_2-t)/TD} u(t) dt.
Now everything depends on how your input current is given in your setting.
Step function? Solve by hand and use the explicit solution.
Numerical function where you have sample points? Use numerical integrator.
Fourier modes? Something else?
(edit: all formulas were wrong on typing then the first time…)
My input is a time series of wind turbine speed measurements, which could be real values or results of a simulation. A time series with 500000 values. I have the speed and want to calculate the acceleration.
I try to port a Simulink model that uses a filtered derivative to Julia to verify my results, but this means I should use the same algorithms…
This is a simple linear filter system with a saturation on the output, DSP.filter or ControlSystems.lsim does this, you just need to add the saturation using clamp on the filter output.
Do you have code available? Or some blackbox implementation? Or at least a sample of input-output values to validate your replication?
Or just a paper with plots?
If its just a paper with plots, then I fear that there is lots of freedom in the details. If you want your replication to be faithful then this will be pain (unless there is an unstated convention for all these details – ask an engineer with domain knowledge or an author of the thing you’re replicating?).
Already “a time series of measurements” – yeah, the filtered derivative thingy you pointed at is a continuous model, so they don’t fit together.
Either somebody must make a modeling decision how to embed “time series of measurements” into a time-continuous function space, i.e. a modeling decision about interpolation.
Or somebody must make a discretization decision for the filtered derivative thingy (instead of deriving it from an interpolation procedure).
But the point of a replication is that you adamantly refuse to make any decisions at all.
PS. Or ask @baggepinnen , maybe I’m overthinking it, coming from the maths side…
No you are correct, the filter shown is indeed in continuous time and some additional assumptions like you describe are required. I didn’t read the whole chapter there, maybe the OP could add all the required details to this discourse post?
Here is a script implementing filtered derivative:
using GLMakie
mutable struct FilteredDerivative{T <: AbstractFloat}
# Filter state
current_state::T
# Filter paramters
Ts::T # Sample time
Kd::T # Filter Gain
Td::T # Filter time constant
A::T # Filter lower bound
B::T # Filter upper bound
function FilteredDerivative(Ts::T, Kd::T, Td::T, A::T, B::T) where T <: AbstractFloat
fd = new{T}()
fd.current_state = 0.0 # Assumes zero initial condition
fd.Ts = Ts
fd.Kd = Kd
fd.Td = Td
fd.A = A
fd.B = B
return fd
end
end
function step(u::T, FD::FilteredDerivative) where T <: AbstractFloat
# Compute the output
y1 = FD.current_state
u1 = (u-y1)
y = (FD.Kd / FD.Td)*u1
if y > B
y = B
elseif y < A
y = A
end
# Update the internal state
FD.current_state = FD.current_state + FD.Ts*(1/FD.Td)*u1
return y
end
Ts = 0.01
Td = 6.0
Kd = 1.0
A = -1.5
B = 1.5
Fd = FilteredDerivative(Ts, Kd, Td, A, B)
signal(t) = t
t = range(start = 0.0, stop = 10.0, step = Ts)
u_in = signal.(t) + 0.2*sin.(100*t) # Input signal with noise
n = length(t)
u_f = similar(u_in) # Filtered output
u_d = (signal.(t[2:end]) - signal.(t[1:end-1]))/Ts
u_nd = (u_in[2:end] - u_in[1:end-1])/Ts
for i = 1:n
global Fd
u_f[i] = step(u_in[i], Fd)
end
f, a, l = lines(t, u_in, label ="Input signal")
lines!(a, t, u_f, label = "FilteredDerivative")
lines!(a, t[1:end-1], u_d, label = "Derivative of signal whithout noise")
axislegend()
f
f1, a1, l1 = lines(t[1:end-1], u_nd, label = "Noisy Derivative")
axislegend()
f1
Well, the Simulink block is both for discrete time and continues time, and they link in the help to this paper… They write, in discrete time the block would be equivalent to the following transfer function:
(K/T)(z-1)/(z+Ts/T-1)
with T the time constant and Ts the sampling time and K the gain.
@A_C Does you example match this discrete time notation?
The formula you quote corresponds a first-order IIR filter. Is there a reason why you want to use this particular filter, vs. a finite difference filter (= FIR filter, as in your other recent thread) in conjunction with smoothing (low-pass filtering) to reduce noise, or a higher-order IIR filter, or …?
In general, filter design is very context dependent. How noisy is your data? Over what bandwidth do you need to be accurate? How will it be implemented, and how efficient does it need to be?
If you have 50000 points that you need to approximately differentiate in Julia, that is a very small-scale problem where efficiency doesn’t sound like a primary concern (as opposed to something that needs to operate continuously in an embedded system). So if I were you, I would implement something straightforward that you can understand completely, e.g. smoothing followed by a finite-difference approximation, instead of trying to copy formulas whose properties you don’t really understand.
But if you want to learn more, I would pick up a textbook on discrete-time signal processing and start by understanding IIR and FIR filters. They are extremely simple to implement once you understand them (and the DSP.jl package, for example, can do this for you) — the hard part is designing the filter weights for a particular application.
Yes…but do compare with MATLAB results and check. This implementation is suitable if you want to implement the D part of the PID controller (It is from the paper that you have linked in the first post). Since the implementation is causal if you have data coming in real time and want to compute the derivative then it can be used. But do note that it will reach the value of derivative asymptotically (as can be seen from the graphs). If you have the complete data you can just filter it appropriately and compute the derivative numerically.
Well, I have the complete data in the simulation, but I am designing a filter that shall be used online…
The point is, I am trying to improve an existing implementation, and therefore I need a correct Julia implementation of what exists first. Yes, I think an FIR filter might be better, but first I need to implement the existing IIR filter correctly as reference…
Furthermore, to have 1:1 translations of Simulink components to Julia is valuable in itself…
Well, your code matches the results from Matlab/Simulink very closely… But it does not look like a discrete time implementation…
How can I implement the discrete time transfer function:
with T being the time constant and T_s the sampling time in Julia directly?
Look at the Wikipedia article on IIR filters, for example.
If you have a transfer function
then it corresponds in (discrete) time domain to the formula:
which you can implement by straightforward loops:
y[1:Q] = # initial values, typically by assuming x[k] = y[k] = 0 for k ≤ 0
for n = Q+1:length(y)
@views y[n] = dot(b, x[n:-1:n-P]) - dot(a, y[n-1:-1:n-Q])
end
In your case, you have
corresponding P = Q = 1 and (in Julia) to b = [1, -1] / T and a = [Tₛ/T - 1]. The loops above simplify to:
y[1] = # initial value, typically b[1] * x[1] - 0
for n = 2:length(y)
y[n] = b[1] * x[n] + b[2] * x[n-1] - a[1] * y[n-1]
end
Looks like magic to me… But I am learning…
And how do you create this nicely formatted formulas?
Enclosing the formula in $$ did the trick… Thank you!
It is an discrete time implementation. I have discritized the differential equations that you have written in original post using euler method. Discretization can be done in many ways. Here, are some ways you can convert continous transfer function to discrete transfer function. So I dont think you can get an unique implementation.
One more question: Is this just a first order high pass filter, or is there anything special about it?