Bode Plot Example

Setup

Mix.install([
  {:nx, "~> 0.13"},
  {:nx_lapack, path: System.fetch_env!("HOME") <> "/Documents/forks/nx_lapack"},
  {:nx_control, path: System.fetch_env!("HOME") <> "/Documents/forks/nx_control"},
  {:vega_lite, "~> 0.1"},
  {:kino_vega_lite, "~> 0.1"}
])

alias VegaLite, as: Vl

System Definition

Consider an open-loop transfer function:

[ G(s) = \frac{1}{s(s+1)(s+10)} ]

tf = NxControl.TransferFunction.new([1], [1, 11, 10, 0])

Frequency Response Data

{freqs, mag, phase} = NxControl.FrequencyResp.bode_data(tf)

Magnitude Plot

data = %{
  log_freq: Nx.to_flat_list(freqs) |> Enum.map(&:math.log10/1),
  magnitude: Nx.to_flat_list(mag)
}

Vl.new(width: 600, height: 300, title: "Bode Plot — Magnitude")
|> Vl.data_from_values(data)
|> Vl.mark(:line)
|> Vl.encode_field(:x, "log_freq", type: :quantitative, title: "log₁₀(ω)")
|> Vl.encode_field(:y, "magnitude", type: :quantitative, title: "Magnitude (dB)")

Phase Plot

data2 = %{
  log_freq: Nx.to_flat_list(freqs) |> Enum.map(&:math.log10/1),
  phase: Nx.to_flat_list(phase)
}

Vl.new(width: 600, height: 300, title: "Bode Plot — Phase")
|> Vl.data_from_values(data2)
|> Vl.mark(:line)
|> Vl.encode_field(:x, "log_freq", type: :quantitative, title: "log₁₀(ω)")
|> Vl.encode_field(:y, "phase", type: :quantitative, title: "Phase (deg)")

Stability Margins

{gm_mag, gm_phase, pm_freq, pm_phase} = NxControl.FrequencyResp.margins(tf)

MATLAB Time-Constant Form Example

Transfer function defined in MATLAB’s (Ts + 1) format:

f0 = [0.82 1];
f1 = [0.337 1];
f2 = [1.64 1];
f3 = [0.07 1];
den1 = conv(f0, f1)
den2 = conv(f2, f3)
sys = tf(den1, den2)
bode(sys)

Equivalent in nx_control:

# Polynomial multiplication (conv)
f0 = [0.82, 1.0]; f1 = [0.337, 1.0]
f2 = [1.64, 1.0]; f3 = [0.07, 1.0]

poly_mul = fn a, b ->
  len = length(a) + length(b) - 1
  for i <- 0..(len - 1) do
    Enum.reduce(0..i, 0.0, fn j, acc ->
      (Enum.at(a, i - j) || 0.0) * (Enum.at(b, j) || 0.0) + acc
    end)
  end
end

num = poly_mul.(f0, f1)  # 0.27634s² + 1.157s + 1.0
den = poly_mul.(f2, f3)  # 0.1148s² + 1.71s + 1.0

tf_m = NxControl.TransferFunction.new(num, den)
{freqs_m, mag_m, phase_m} = NxControl.FrequencyResp.bode_data(tf_m)

Bode Plot

data_m = %{
  log_freq: Nx.to_flat_list(freqs_m) |> Enum.map(&:math.log10/1),
  magnitude: Nx.to_flat_list(mag_m)
}

Vl.new(width: 600, height: 300, title: "MATLAB Transfer Function — Magnitude")
|> Vl.data_from_values(data_m)
|> Vl.mark(:line, color: "green")
|> Vl.encode_field(:x, "log_freq", type: :quantitative, title: "log₁₀(ω)")
|> Vl.encode_field(:y, "magnitude", type: :quantitative, title: "Magnitude (dB)")

data_m2 = %{
  log_freq: Nx.to_flat_list(freqs_m) |> Enum.map(&:math.log10/1),
  phase: Nx.to_flat_list(phase_m)
}

Vl.new(width: 600, height: 300, title: "MATLAB Transfer Function — Phase")
|> Vl.data_from_values(data_m2)
|> Vl.mark(:line, color: "green")
|> Vl.encode_field(:x, "log_freq", type: :quantitative, title: "log₁₀(ω)")
|> Vl.encode_field(:y, "phase", type: :quantitative, title: "Phase (deg)")