MATLAB has a nice library of digital signals processing functions that give you a lot of control over linear-time invariant systems. My goal is to view the frequency response of a system governed by two vectors a=[1, -1.1, 0.46]
and b=[1, 0.04, 0.76]
.
To do this in MATLAB, I'd just use the freqz(b,a)
and get back the frequency response h
and frequencies evaluated w
. I am working in Julia and am using the DSP
package, so I'm attempting to use the DSP.Filters.freqresp
function. To set this up, I define the vectors and import the relevant package:
a=[1, -1.1, 0.46]
b=[1, 0.04, 0.76]
using DSP.Filters : freqresp
h, w = freqresp([b,a])
When I run this I get the error
MethodError: no method matching freqresp(::Vector{Int64})
How should I be implementing this function to obtain a valid result?
@Shayan pointed me in the right direction by mentioning PolynomialRatio
, as plugging in the a
and b
vectors give a representation that freqresp
will accept. As a result, I get the same answer as in MATLAB
a=[1, -1.1, 0.46]
b=[1, 0.04, 0.76]
using DSP
z = PolynomialRatio(b,a)
h, w = freqresp(z)
mag = abs.(h)
phase = atan.(imag.(h)./real.(h))*180/pi
p1 = plot(w, mag)
xaxis!("")
yaxis!("Magnitude")
p2 = plot(w, phase)
xaxis!("Normalized Frequency")
yaxis!("Wrapped Phase [deg]")
plot(p1, p2, layout=(2,1))
Which gives
You can also make use of the ControlSystems.jl package for this:
julia> using ControlSystemsBase, Plots
julia> a=[1, -1.1, 0.46]
3-element Vector{Float64}:
1.0
-1.1
0.46
julia> b=[1, 0.04, 0.76]
3-element Vector{Float64}:
1.0
0.04
0.76
julia> Z = tf(b,a,1)
TransferFunction{Discrete{Int64}, ControlSystemsBase.SisoRational{Float64}}
1.0z^2 + 0.04z + 0.76
---------------------
1.0z^2 - 1.1z + 0.46
Sample Time: 1 (seconds)
Discrete-time transfer function model
julia> bodeplot(Z)
The bodeplot
plot recipe handles phase unrolling etc. for you.
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With