ControlSystems / Frequency in Hz for bodeplot

Specific Domains Modelling & Simulations
,
1

Hello,

I would like to know how to change the xaxis unit (rad/s) to frequency for bodeplot.

Thanks in advance.
Patrick.

2

I’m no expert, but a clarifying question…:

  • is it purely the xaxis text you want to change, or is it both the text and the scaling (so that the abcissa actually is the frequecy in Hz)?

One of my frustrations with MATLAB’s Control TB is that it defaults to β€œs” (step, impulse) and β€œrad/s” for Bodeplots… I often consider systems where β€œmin”, β€œh”, β€œday”, β€œweek” – or β€œms” for that matter - is the time unit. I’m sure it is possible to override this in MATLAB, but it is not intuitive how to do this – at least not the way my brain works.

Hopefully, this is simpler in ControlSystems of Julia. If it is only to change the text, my first attempt would be just specify xlabel = "Hz". If you need to scale the axis, a cumbersome way is to scale the time constants, etc with 2*pi.

3

Hi BLI,

I want also to change the scale. I am an electronic designer and we are dealing with frequency in Hz, not the rad/s.

Patrick.

4

You can always do

mag, phase, w = bode(sys)
plot(w./(2pi), mag)

and if on master branch, I think you can do

mag, phase, w = bode(sys)
plot(w./(2pi), mag, seriestype=:bodemag)

(not entierly sure about the keyword seriestype though)

5

Also you could check out the SampledSignals package, which is designed to represent regularly-sampled data and is samplerate-aware. the SampleBuf type holds time-series data, and when you take the fft of it you get a SpectrumBuf. the domain function will give you the x-axis values, either in seconds or Hz, respectively. So plotting the spectrum is just

using Plots
using SampledSignals

julia> buf = SampleBuf(rand(1024) .* sin.(linspace(0,2Ο€,1024)), 48000)
1024-frame, 1-channel SampleBuf{Float64, 1}
0.021333333333333333s sampled at 48000.0Hz
β–„β–…β–†β–†β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–†β–†β–†β–…β–„β–…β–†β–†β–†β–†β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–‡β–†β–†β–†β–…β–„


julia> spec = fft(buf)
1024-frame, 1-channel SpectrumBuf{Complex{Float64}, 1}
48000.0Hz sampled at 0.021333333333333333s
β–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆβ–ˆ


julia> plot(domain(spec), abs.(spec))

There are still some rough edges, but there are a couple of us working to smooth them out. I’ve also been pretty slow to get 0.7 compatibility working, but that should be coming in the next couple weeks at the latest.

6

This looks like an ok solution. I assume there is a similar :bodephase seriestype?

Is it ok to just set the xlabel? (I’m away from laptop, snd cannot run Julia on my smartphone… :slight_smile: ).

7

Yes, all plot commands works just as in Plots.jl

8

You have the code and can inspect the available options here. It’s all based on PlotRecipes.

9

Testing plotting of ControlSystems systems… I have used julia> using ControlSystems. Then

julia> sys = tfg("2*(1-3*s)/(1+0.1*s)/(1+10*s)");
#
julia> bodeplot(sys,xlabel=["" "Frequency (Hz)"],linecolor = [:red :blue])

The result is as follows:
image

… which looks fine, except that the specified linecolor is not adhered to. [The reason for specifying individual xlabel for magnitude and phase, is that I don’t want to waste space by inserting xlabel for the uppermost plot…]

Next, I tried to use the bode command in combination with standard plot.

mag,phase,w = bode(sys);

Turns out that mag and phase are not standard vectors, but 3D arrays:

julia> mag
400Γ—1Γ—1 Array{Float64,3}:
[:, :, 1] =
 1.99097  
 1.99055  
...

Thus, to plot the magnitude…

julia> plot(w,20*log10.(mag[:][:][:]),xscale=:log10)

leads to:
image

So, two questions:

  1. Why doesn’t it work to specify different colors for magnitude and phase?
  2. Why are mag and phase returned as 3D arrays? [This slightly complicates plotting…]

PS: with a little bit of fine tuning, it is possible to get decent Bode diagrams:

julia> plot(w,[20*log10.(mag[:][:][:]),phase[:][:][:]],xscale=:log10,linecolor=[:red :blue],layout=(2,1))
julia> plot!(title=["Bode plot" ""],xlabel=["" "Frequency [rad/s]"])
julia> plot!(ylabel=["Magnitude (dB)" "Phase (deg)"])

leads to:
image

10

We have been rewriting a lot of the functionality in ControlSystems.jl, some of the changes are on the β€œmaster” branch of the package and not in a tagged version. It seems to me that the color works as expected on the β€œmaster” branch. You could try

Pkg.checkout("ControlSystems", "master")

and see if it works. There are other major updates on its way on the β€œautodiff” branch, but they are not completely stable yet.

For the question regarding the 400Γ—1Γ—1 Array:
This is the result of the system being a 1-input 1-output system, so to access the response from input j to output i on a general system you do mag[:,i,j] and get a simple vector.

There have been discussions on automatically reduce this to a vector in some SISO cases, but there is no consensus yet, and it has some implications on the type system, input is welcome.

I would also advice against using the tfg method, this will probably be removed soon. It was intended only for cases where the system is not a simple finite order system, as in sys = tfg("e^(-s)/s") and offers limited functionality. The easiest way to create a TransferFunction is to first create the system β€œs” using the special constructor tf("s"):

s = tf("s") # Equivalent to tf([1, 0], [1])
sys = 2*(1-3*s)/(1+0.1*s)/(1+10*s)

or if you know the coefficients

sys = tf([-6.0, 2], [1, 10.1, 1]) # (-6s+2)*(1s^2+10.1s+1)

or as a zero-pole-gain system

sys = zpk([1/3], [-10, -0.1], -6) # -6*(s-1/3)/((s+10)*(s+0.1))

A last note: mag[:] will β€œvectorize” the array, so there is no need to do mag[:][:][:], they are equivalent.

Edit: this command should generated the plot you want on the β€œmaster” branch

bodeplot(sys, xlabel=["" "Frequency (Hz)"], linecolor=[:red :blue], title=["Bode plot" ""])
1 Like