Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Equivalent to MATLAB's freqz in Julia

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?

like image 237
t.o. Avatar asked Sep 17 '25 13:09

t.o.


2 Answers

@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

enter image description here

like image 133
t.o. Avatar answered Sep 19 '25 06:09

t.o.


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)

enter image description here

The bodeplot plot recipe handles phase unrolling etc. for you.

like image 21
Fredrik Bagge Avatar answered Sep 19 '25 06:09

Fredrik Bagge