diff --git a/src/WRPLL_freqsim_example.py b/src/WRPLL_freqsim_example.py
new file mode 100644
index 0000000..29528f7
--- /dev/null
+++ b/src/WRPLL_freqsim_example.py
@@ -0,0 +1,55 @@
+# %%
+import control as ct
+import numpy as np
+
+# %%
+# Gain term
+N = 0x8000
+F_s = 125e6 / N  # Hz
+T = 1 / F_s
+
+K_dxco = 0.0001164 * (2**np.pi) * N / (1e6)
+K_pd = N / (2**np.pi)
+
+# PI controller
+K_P, K_I = 6, 2
+
+sims = [
+    ["KP=6,KI=2", (K_P + K_I), -K_P, 0, 1, -1, 0],
+    [
+        # from https://www.earlevel.com/main/2021/09/02/biquad-calculator-v3/
+        "Biquad Low pass",
+        0.07209205036273991,
+        0.14418410072547982,
+        0.07209205036273991,
+        1,
+        -0.6114078511562919,
+        -0.10022394739274834,
+    ],
+]
+
+open_loops = []
+close_loops = []
+
+for [name, b0, b1, b2, a0, a1, a2] in sims:
+    ph_gain = ct.tf([K_pd], [1], dt=T)
+    loop_filter = ct.tf([b0, b1, b2], [a0, a1, a2], dt=T)
+    dcxo = ct.tf([K_dxco], [1, -1], dt=T)
+
+    open_tf = ct.series(ph_gain, loop_filter, dcxo)
+    close_tf = ct.feedback(open_tf, 1)
+    gm, pm, _, _ = ct.margin(open_tf)
+    print(f"{name} | gain margin:{gm} phase margin:{pm}")
+
+    # HACK labeling graph
+    open_tf.name = name
+    close_tf.name = name
+    open_loops.append(open_tf)
+    close_loops.append(close_tf)
+
+ct.bode_plot(open_loops, dB=True, Hz=True, title="Open loop bode")
+ct.bode_plot(
+    close_loops, dB=True, Hz=True, plot_phase=False, title="Close loop response"
+)
+
+# %%