|
| 1 | +""" pyplots.ai |
| 2 | +bode-basic: Bode Plot for Frequency Response |
| 3 | +Library: letsplot 4.9.0 | Python 3.14.3 |
| 4 | +Quality: 90/100 | Created: 2026-03-21 |
| 5 | +""" |
| 6 | + |
| 7 | +import numpy as np |
| 8 | +import pandas as pd |
| 9 | +from lets_plot import * # noqa: F403 |
| 10 | +from lets_plot.export import ggsave as export_ggsave |
| 11 | +from scipy import signal |
| 12 | + |
| 13 | + |
| 14 | +LetsPlot.setup_html() # noqa: F405 |
| 15 | + |
| 16 | +# Data - 3rd-order system: G(s) = 500 / ((s+20)(s^2+2s+25)) |
| 17 | +# Complex poles (wn=5, zeta=0.2) give a prominent resonance peak ~8 dB |
| 18 | +# Real pole at s=-20 preserves resonance shape |
| 19 | +# Stable system with GM≈5 dB and PM≈18° |
| 20 | +num = [500.0] |
| 21 | +den = np.polymul([1, 20], [1, 2.0, 25]) |
| 22 | +system = signal.TransferFunction(num, den) |
| 23 | + |
| 24 | +omega = np.logspace(-1, 2.5, 500) |
| 25 | +_, mag, phase_deg = signal.bode(system, omega) |
| 26 | +frequency_hz = omega / (2 * np.pi) |
| 27 | + |
| 28 | +# Build DataFrames |
| 29 | +df_mag = pd.DataFrame({"frequency_hz": frequency_hz, "magnitude_db": mag}) |
| 30 | +df_phase = pd.DataFrame({"frequency_hz": frequency_hz, "phase_deg": phase_deg}) |
| 31 | + |
| 32 | +# Gain crossover frequency (where magnitude crosses 0 dB) |
| 33 | +zero_crossings = np.where(np.diff(np.sign(mag)))[0] |
| 34 | +if len(zero_crossings) > 0: |
| 35 | + idx_gc = zero_crossings[0] |
| 36 | + t = abs(mag[idx_gc]) / (abs(mag[idx_gc]) + abs(mag[idx_gc + 1])) |
| 37 | + freq_gc = frequency_hz[idx_gc] + t * (frequency_hz[idx_gc + 1] - frequency_hz[idx_gc]) |
| 38 | + phase_at_gc = phase_deg[idx_gc] + t * (phase_deg[idx_gc + 1] - phase_deg[idx_gc]) |
| 39 | + phase_margin = 180 + phase_at_gc |
| 40 | +else: |
| 41 | + freq_gc = None |
| 42 | + phase_margin = None |
| 43 | + |
| 44 | +# Phase crossover frequency (where phase crosses -180 deg) |
| 45 | +phase_crossings = np.where(np.diff(np.sign(phase_deg + 180)))[0] |
| 46 | +if len(phase_crossings) > 0: |
| 47 | + idx_pc = phase_crossings[0] |
| 48 | + t_pc = abs(phase_deg[idx_pc] + 180) / (abs(phase_deg[idx_pc] + 180) + abs(phase_deg[idx_pc + 1] + 180)) |
| 49 | + freq_pc = frequency_hz[idx_pc] + t_pc * (frequency_hz[idx_pc + 1] - frequency_hz[idx_pc]) |
| 50 | + mag_at_pc = mag[idx_pc] + t_pc * (mag[idx_pc + 1] - mag[idx_pc]) |
| 51 | + gain_margin = -mag_at_pc |
| 52 | +else: |
| 53 | + freq_pc = None |
| 54 | + gain_margin = None |
| 55 | + |
| 56 | +# Colors - colorblind-safe palette (no red-green) |
| 57 | +COLOR_MAIN = "#1B4F72" |
| 58 | +COLOR_REF = "#AEB6BF" |
| 59 | +COLOR_GM = "#D35400" # Orange for gain margin |
| 60 | +COLOR_PM = "#2471A3" # Blue for phase margin |
| 61 | + |
| 62 | +# Common theme |
| 63 | +common_theme = flavor_high_contrast_light() + theme( # noqa: F405 |
| 64 | + axis_text=element_text(size=16), # noqa: F405 |
| 65 | + axis_title=element_text(size=20, face="bold"), # noqa: F405 |
| 66 | + panel_grid_major=element_line(color="#E0E0E0", size=0.3), # noqa: F405 |
| 67 | + panel_grid_minor=element_blank(), # noqa: F405 |
| 68 | + panel_border=element_blank(), # noqa: F405 |
| 69 | + axis_line=element_line(color="#AAAAAA", size=0.5), # noqa: F405 |
| 70 | + plot_margin=[30, 30, 10, 20], |
| 71 | +) |
| 72 | + |
| 73 | +# Magnitude plot |
| 74 | +mag_plot = ( |
| 75 | + ggplot(df_mag, aes(x="frequency_hz", y="magnitude_db")) # noqa: F405 |
| 76 | + + geom_line( # noqa: F405 |
| 77 | + color=COLOR_MAIN, |
| 78 | + size=2.5, |
| 79 | + tooltips=layer_tooltips() # noqa: F405 |
| 80 | + .format("frequency_hz", ".3f") |
| 81 | + .format("magnitude_db", ".1f") |
| 82 | + .line("Freq: @frequency_hz Hz") |
| 83 | + .line("Mag: @magnitude_db dB"), |
| 84 | + ) |
| 85 | + + geom_hline(yintercept=0, color=COLOR_REF, size=0.8, linetype="dashed") # noqa: F405 |
| 86 | + + scale_x_log10() # noqa: F405 |
| 87 | + + scale_y_continuous(limits=[-35, max(mag) + 3]) # noqa: F405 |
| 88 | + + labs( # noqa: F405 |
| 89 | + x="", y="Magnitude (dB)", title="bode-basic \u00b7 letsplot \u00b7 pyplots.ai" |
| 90 | + ) |
| 91 | + + theme(plot_title=element_text(size=26, hjust=0.5, face="bold")) # noqa: F405 |
| 92 | + + common_theme |
| 93 | + + ggsize(1600, 450) # noqa: F405 |
| 94 | +) |
| 95 | + |
| 96 | +# Add gain margin annotation if phase crossover exists |
| 97 | +if freq_pc is not None: |
| 98 | + gm_seg = pd.DataFrame({"x": [freq_pc], "y": [mag_at_pc], "xend": [freq_pc], "yend": [0.0]}) |
| 99 | + gm_pts = pd.DataFrame({"x": [freq_pc, freq_pc], "y": [0.0, mag_at_pc]}) |
| 100 | + gm_label = pd.DataFrame({"x": [freq_pc * 3.0], "y": [mag_at_pc - 5], "label": [f"GM = {gain_margin:.1f} dB"]}) |
| 101 | + mag_plot = ( |
| 102 | + mag_plot |
| 103 | + + geom_vline(xintercept=freq_pc, color=COLOR_GM, size=0.5, linetype="dotted", alpha=0.5) # noqa: F405 |
| 104 | + + geom_segment( # noqa: F405 |
| 105 | + aes(x="x", y="y", xend="xend", yend="yend"), # noqa: F405 |
| 106 | + data=gm_seg, |
| 107 | + color=COLOR_GM, |
| 108 | + size=3.0, |
| 109 | + arrow=arrow(type="closed", length=8, ends="both"), # noqa: F405 |
| 110 | + ) |
| 111 | + + geom_point( # noqa: F405 |
| 112 | + aes(x="x", y="y"), # noqa: F405 |
| 113 | + data=gm_pts, |
| 114 | + color=COLOR_GM, |
| 115 | + size=10, |
| 116 | + shape=18, |
| 117 | + ) |
| 118 | + + geom_label( # noqa: F405 |
| 119 | + aes(x="x", y="y", label="label"), # noqa: F405 |
| 120 | + data=gm_label, |
| 121 | + size=18, |
| 122 | + color=COLOR_GM, |
| 123 | + fill="#FDEBD0", |
| 124 | + label_padding=0.3, |
| 125 | + label_r=0.15, |
| 126 | + fontface="bold", |
| 127 | + ) |
| 128 | + ) |
| 129 | + |
| 130 | +# Add gain crossover point on magnitude plot |
| 131 | +if freq_gc is not None: |
| 132 | + gc_mag_pt = pd.DataFrame({"x": [freq_gc], "y": [0.0]}) |
| 133 | + mag_plot = mag_plot + geom_point( # noqa: F405 |
| 134 | + aes(x="x", y="y"), # noqa: F405 |
| 135 | + data=gc_mag_pt, |
| 136 | + color=COLOR_PM, |
| 137 | + size=10, |
| 138 | + shape=18, |
| 139 | + ) |
| 140 | + |
| 141 | +# Phase plot |
| 142 | +phase_plot = ( |
| 143 | + ggplot(df_phase, aes(x="frequency_hz", y="phase_deg")) # noqa: F405 |
| 144 | + + geom_line( # noqa: F405 |
| 145 | + color=COLOR_MAIN, |
| 146 | + size=2.5, |
| 147 | + tooltips=layer_tooltips() # noqa: F405 |
| 148 | + .format("frequency_hz", ".3f") |
| 149 | + .format("phase_deg", ".1f") |
| 150 | + .line("Freq: @frequency_hz Hz") |
| 151 | + .line("Phase: @phase_deg\u00b0"), |
| 152 | + ) |
| 153 | + + geom_hline(yintercept=-180, color=COLOR_REF, size=0.8, linetype="dashed") # noqa: F405 |
| 154 | + + scale_x_log10() # noqa: F405 |
| 155 | + + labs(x="Frequency (Hz)", y="Phase (\u00b0)") # noqa: F405 |
| 156 | + + common_theme |
| 157 | + + ggsize(1600, 450) # noqa: F405 |
| 158 | +) |
| 159 | + |
| 160 | +# Add phase margin annotation if gain crossover exists |
| 161 | +if freq_gc is not None: |
| 162 | + pm_seg = pd.DataFrame({"x": [freq_gc], "y": [phase_at_gc], "xend": [freq_gc], "yend": [-180.0]}) |
| 163 | + pm_pts = pd.DataFrame({"x": [freq_gc, freq_gc], "y": [-180.0, phase_at_gc]}) |
| 164 | + pm_label = pd.DataFrame( |
| 165 | + {"x": [freq_gc * 2.5], "y": [phase_at_gc + 10], "label": [f"PM = {phase_margin:.1f}\u00b0"]} |
| 166 | + ) |
| 167 | + phase_plot = ( |
| 168 | + phase_plot |
| 169 | + + geom_vline(xintercept=freq_gc, color=COLOR_PM, size=0.5, linetype="dotted", alpha=0.5) # noqa: F405 |
| 170 | + + geom_segment( # noqa: F405 |
| 171 | + aes(x="x", y="y", xend="xend", yend="yend"), # noqa: F405 |
| 172 | + data=pm_seg, |
| 173 | + color=COLOR_PM, |
| 174 | + size=3.0, |
| 175 | + arrow=arrow(type="closed", length=8, ends="both"), # noqa: F405 |
| 176 | + ) |
| 177 | + + geom_point( # noqa: F405 |
| 178 | + aes(x="x", y="y"), # noqa: F405 |
| 179 | + data=pm_pts, |
| 180 | + color=COLOR_PM, |
| 181 | + size=10, |
| 182 | + shape=18, |
| 183 | + ) |
| 184 | + + geom_label( # noqa: F405 |
| 185 | + aes(x="x", y="y", label="label"), # noqa: F405 |
| 186 | + data=pm_label, |
| 187 | + size=18, |
| 188 | + color=COLOR_PM, |
| 189 | + fill="#D4E6F1", |
| 190 | + label_padding=0.3, |
| 191 | + label_r=0.15, |
| 192 | + fontface="bold", |
| 193 | + ) |
| 194 | + ) |
| 195 | + |
| 196 | +# Add phase crossover point on phase plot |
| 197 | +if freq_pc is not None: |
| 198 | + pc_phase_pt = pd.DataFrame({"x": [freq_pc], "y": [-180.0]}) |
| 199 | + phase_plot = phase_plot + geom_point( # noqa: F405 |
| 200 | + aes(x="x", y="y"), # noqa: F405 |
| 201 | + data=pc_phase_pt, |
| 202 | + color=COLOR_GM, |
| 203 | + size=10, |
| 204 | + shape=18, |
| 205 | + ) |
| 206 | + |
| 207 | +# Combine vertically |
| 208 | +combined = ggbunch( # noqa: F405 |
| 209 | + plots=[mag_plot, phase_plot], regions=[(0, 0, 1, 0.53, 0, 0), (0, 0.47, 1, 0.55, 0, 0)] |
| 210 | +) |
| 211 | + |
| 212 | +# Save |
| 213 | +export_ggsave(combined, filename="plot.png", path=".", scale=3) |
| 214 | +export_ggsave(combined, filename="plot.html", path=".") |
0 commit comments