import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
# ============================================================
# SETTINGS
# ============================================================
filename = Path(
r"C:\Users\roger\Documents\cst_prj\diff_line_new_diff_two_4_lines_1.s4p"
)
VIN_STEP = 1.0
DF = 10e6 # 10 MHz frequency spacing
T_MIN_PS = -500
T_MAX_PS = 3000
SETTLING_TOL = 0.02 # ±2%
# ============================================================
# READ S4P
#
# Touchstone ordering:
#
# S11 S21 S31 S41
# S12 S22 S32 S42
# S13 S23 S33 S43
# S14 S24 S34 S44
#
# ============================================================
def read_s4p(filename):
values = []
freq_unit = "ghz"
data_format = "ma"
unit_scale = {
"hz": 1.0,
"khz": 1e3,
"mhz": 1e6,
"ghz": 1e9
}
with open(filename, "r") as file:
for line in file:
# Remove comments
line = line.split("!")[0].strip()
if not line:
continue
# Touchstone option line
if line.startswith("#"):
parts = line.lower().split()
freq_unit = parts[1]
data_format = parts[3]
continue
# Ignore Touchstone 2.0 keyword lines
if line.startswith("["):
continue
values.extend(
float(x) for x in line.split()
)
values = np.asarray(values)
# Each 4-port frequency point has:
#
# 1 frequency
# +
# 16 complex S-parameters × 2 numbers
#
# = 33 numbers
if len(values) % 33 != 0:
raise ValueError(
"Could not interpret the file as a standard S4P file."
)
data = values.reshape(-1, 33)
f = data[:, 0] * unit_scale[freq_unit]
def convert(a, b):
if data_format == "ma":
return (
a
* np.exp(
1j * np.deg2rad(b)
)
)
elif data_format == "db":
mag = 10 ** (a / 20)
return (
mag
* np.exp(
1j * np.deg2rad(b)
)
)
elif data_format == "ri":
return a + 1j*b
else:
raise ValueError(
f"Unknown Touchstone format: {data_format}"
)
# First column of S-matrix
#
# S11 -> columns 1,2
# S21 -> columns 3,4
# S31 -> columns 5,6
# S41 -> columns 7,8
s11 = convert(
data[:, 1],
data[:, 2]
)
s21 = convert(
data[:, 3],
data[:, 4]
)
s31 = convert(
data[:, 5],
data[:, 6]
)
s41 = convert(
data[:, 7],
data[:, 8]
)
return f, s11, s21, s31, s41
# ============================================================
# LOAD FILE
# ============================================================
f, s11, s21, s31, s41 = read_s4p(filename)
print("\nFILE INFORMATION")
print("----------------")
print(f"Points = {len(f)}")
print(f"Start frequency = {f[0]/1e9:.3f} GHz")
print(f"Stop frequency = {f[-1]/1e9:.3f} GHz")
print(f"Frequency step = {(f[1]-f[0])/1e6:.3f} MHz")
# ============================================================
# CHECK FREQUENCY SPACING
# ============================================================
measured_df = np.diff(f)
if not np.allclose(
measured_df,
DF,
rtol=1e-6,
atol=1.0
):
raise ValueError(
"Frequency points are not uniformly spaced by 10 MHz."
)
# ============================================================
# ESTIMATE DC
#
# Linear complex extrapolation from first two measured points.
#
# For a real time-domain response, the DC bin must be real.
# ============================================================
def estimate_dc(f, s):
slope = (
s[1] - s[0]
) / (
f[1] - f[0]
)
dc_complex = (
s[0]
- slope * f[0]
)
return complex(
np.real(dc_complex),
0.0
)
s11_dc = estimate_dc(f, s11)
s21_dc = estimate_dc(f, s21)
s31_dc = estimate_dc(f, s31)
s41_dc = estimate_dc(f, s41)
print("\nESTIMATED DC")
print("------------")
print("S11(0) =", s11_dc)
print("S21(0) =", s21_dc)
print("S31(0) =", s31_dc)
print("S41(0) =", s41_dc)
# ============================================================
# ADD DC POINT
# ============================================================
f = np.insert(
f,
0,
0.0
)
s11 = np.insert(
s11,
0,
s11_dc
)
s21 = np.insert(
s21,
0,
s21_dc
)
s31 = np.insert(
s31,
0,
s31_dc
)
s41 = np.insert(
s41,
0,
s41_dc
)
# ============================================================
# CREATE TWO-SIDED HERMITIAN SPECTRUM
#
# Positive:
#
# 0, 10 MHz, ... 67 GHz
#
# Negative:
#
# -67 GHz, ... -10 MHz
#
# This produces an odd-length FFT record.
# ============================================================
def make_full_spectrum(s):
return np.concatenate(
[
s,
np.conj(
s[-1:0:-1]
)
]
)
S11_full = make_full_spectrum(s11)
S21_full = make_full_spectrum(s21)
S31_full = make_full_spectrum(s31)
S41_full = make_full_spectrum(s41)
N = len(S21_full)
# ============================================================
# TIME AXIS
# ============================================================
dt = 1.0 / (N * DF)
# After fftshift:
#
# -50 ns ........ 0 ........ +50 ns
#
t = (
np.arange(N)
- N // 2
) * dt
t_ps = t * 1e12
print("\nTIME DOMAIN")
print("-----------")
print(f"FFT points = {N}")
print(f"Time spacing = {dt*1e12:.4f} ps")
print(f"Full time period = {1/DF*1e9:.2f} ns")
print(
f"Centered axis = "
f"{t[0]*1e9:.2f} ns "
f"to {t[-1]*1e9:.2f} ns"
)
# ============================================================
# FFT FREQUENCY AXIS
# ============================================================
freq_full = np.fft.fftfreq(
N,
d=dt
)
FMAX = f[-1]
# ============================================================
# HAMMING FREQUENCY WINDOW
#
# w(0) = 1
#
# w(±FMAX) ≈ 0.08
#
# ============================================================
x = (
np.abs(freq_full)
/ FMAX
)
w = (
0.54
+
0.46
* np.cos(
np.pi * x
)
)
w[x > 1.0] = 0.0
# ============================================================
# STEP RESPONSE
#
# IMPORTANT:
#
# Raw IFFT ordering:
#
# 0 -> positive time -> negative time
#
# fftshift changes this to:
#
# negative time -> 0 -> positive time
#
# THEN cumsum performs the step integral in chronological order.
# ============================================================
def step_response(S):
# Frequency-domain window
S_windowed = (
S * w
)
# Raw IFFT
impulse_raw = np.fft.ifft(
S_windowed
)
# Put time in chronological order
impulse_shifted = np.fft.fftshift(
impulse_raw
)
# Hermitian spectrum should give real time response
impulse_shifted = np.real(
impulse_shifted
)
# Step response = integral of impulse response
step = np.cumsum(
impulse_shifted
)
return VIN_STEP * step
# ============================================================
# CALCULATE RESPONSES
# ============================================================
V11 = step_response(
S11_full
)
V21 = step_response(
S21_full
)
V31 = step_response(
S31_full
)
V41 = step_response(
S41_full
)
# ============================================================
# FINAL VALUE CHECK
#
# For a 1 V step:
#
# final voltage = Re[Sij(0)] × 1 V
#
# ============================================================
print("\nFINAL VALUE CHECK")
print("-----------------")
print(
f"S21 final = {V21[-1]:.8f} V"
)
print(
f"S21 expected = {np.real(s21_dc)*VIN_STEP:.8f} V"
)
print()
print(
f"S31 final = {V31[-1]:.8f} V"
)
print(
f"S31 expected = {np.real(s31_dc)*VIN_STEP:.8f} V"
)
print()
print(
f"S41 final = {V41[-1]:.8f} V"
)
print(
f"S41 expected = {np.real(s41_dc)*VIN_STEP:.8f} V"
)
print()
print(
f"S11 final = {V11[-1]:.8f} V"
)
print(
f"S11 expected = {np.real(s11_dc)*VIN_STEP:.8f} V"
)
# ============================================================
# 2% SETTLING TIME
#
# Settling time is still CALCULATED,
# but it will NOT be drawn automatically on the graphs.
# ============================================================
def settling_time(
response,
final_value,
tolerance=0.02
):
final_value = (
np.real(final_value)
* VIN_STEP
)
band = (
tolerance
* abs(final_value)
)
lower = (
final_value
- band
)
upper = (
final_value
+ band
)
# Only examine t >= 0
positive_indices = np.where(
t >= 0
)[0]
response_positive = response[
positive_indices
]
inside = (
(response_positive >= lower)
&
(response_positive <= upper)
)
# Work backward:
#
# True means:
# every sample after this point
# remains inside the band.
stable = np.logical_and.accumulate(
inside[::-1]
)[::-1]
valid = np.where(
stable
)[0]
if len(valid) == 0:
return (
None,
final_value,
lower,
upper
)
index = positive_indices[
valid[0]
]
return (
t_ps[index],
final_value,
lower,
upper
)
# ============================================================
# CALCULATE SETTLING VALUES
# ============================================================
ts21, ss21, low21, high21 = settling_time(
V21,
s21_dc,
SETTLING_TOL
)
ts31, ss31, low31, high31 = settling_time(
V31,
s31_dc,
SETTLING_TOL
)
ts41, ss41, low41, high41 = settling_time(
V41,
s41_dc,
SETTLING_TOL
)
ts11, ss11, low11, high11 = settling_time(
V11,
s11_dc,
SETTLING_TOL
)
print("\n2% SETTLING TIMES")
print("-----------------")
print(
f"S21 = {ts21} ps"
)
print(
f"S31 = {ts31} ps"
)
print(
f"S41 = {ts41} ps"
)
print(
f"S11 = {ts11} ps"
)
# ============================================================
# DATA FOR PLOTS
# ============================================================
responses = [
V21,
V31,
V41,
V11
]
names = [
"S21",
"S31",
"S41",
"S11"
]
ylabels = [
"Port 2 voltage [V]",
"Port 3 voltage [V]",
"Port 4 voltage [V]",
"Reflected voltage [V]"
]
steady_values = [
ss21,
ss31,
ss41,
ss11
]
lower_values = [
low21,
low31,
low41,
low11
]
upper_values = [
high21,
high31,
high41,
high11
]
# ============================================================
# CREATE FOUR STACKED PLOTS
# ============================================================
fig, ax = plt.subplots(
4,
1,
figsize=(13, 13),
sharex=True
)
for i in range(4):
# --------------------------------------------------------
# Step response
# --------------------------------------------------------
ax[i].plot(
t_ps,
responses[i],
linewidth=1.6
)
# --------------------------------------------------------
# t = 0 reference
# --------------------------------------------------------
ax[i].axvline(
0,
linestyle=":",
linewidth=1.0,
alpha=0.5
)
# --------------------------------------------------------
# Final DC value
# --------------------------------------------------------
ax[i].axhline(
steady_values[i],
linestyle="--",
linewidth=1.2,
label=(
f"DC final = "
f"{steady_values[i]:.6f} V"
)
)
# --------------------------------------------------------
# ±2% band
# --------------------------------------------------------
ax[i].axhspan(
lower_values[i],
upper_values[i],
alpha=0.12,
label="±2% band"
)
# --------------------------------------------------------
# IMPORTANT:
#
# NO settling-time vertical line
# NO settling-time annotation
#
# Nothing is marked automatically at startup.
# --------------------------------------------------------
ax[i].set_ylabel(
ylabels[i]
)
ax[i].set_title(
f"{names[i]} Step Response"
)
ax[i].grid(
True,
alpha=0.3
)
ax[i].minorticks_on()
ax[i].legend(
loc="best"
)
# ============================================================
# X AXIS
# ============================================================
ax[-1].set_xlabel(
"Time [ps]"
)
ax[-1].set_xlim(
T_MIN_PS,
T_MAX_PS
)
# ============================================================
# MAIN TITLE
# ============================================================
fig.suptitle(
"1 V Step at Port 1 — Hamming Window\n"
"Shifted IFFT Time Data",
fontsize=14
)
# ============================================================
# PERMANENT SYNCHRONIZED CLICK MARKERS
#
# LEFT CLICK ON ANY GRAPH:
#
# same time point is marked on all four graphs.
#
# There are NO markers before the first click.
# ============================================================
marker_number = 0
def onclick(event):
global marker_number
# Left mouse only
if event.button != 1:
return
# Click must be inside one of the plots
if event.inaxes not in ax:
return
if event.xdata is None:
return
# --------------------------------------------------------
# Find nearest actual IFFT time sample
# --------------------------------------------------------
index = np.argmin(
np.abs(
t_ps
- event.xdata
)
)
selected_time = (
t_ps[index]
)
values = [
V21[index],
V31[index],
V41[index],
V11[index]
]
marker_number += 1
# --------------------------------------------------------
# Put same-time marker on all four plots
# --------------------------------------------------------
for i in range(4):
ax[i].plot(
selected_time,
values[i],
marker="o",
markersize=7
)
ax[i].axvline(
selected_time,
linestyle=":",
linewidth=0.8,
alpha=0.5
)
ax[i].annotate(
(
f"M{marker_number}\n"
f"t = {selected_time:.2f} ps\n"
f"V = {values[i]:.6f} V"
),
xy=(
selected_time,
values[i]
),
xytext=(
12,
15
),
textcoords="offset points",
bbox=dict(
boxstyle="round",
alpha=0.8
),
arrowprops=dict(
arrowstyle="->"
)
)
# --------------------------------------------------------
# Terminal output
# --------------------------------------------------------
print(
f"\nMarker M{marker_number}"
)
print(
f"Time = {selected_time:.3f} ps"
)
print(
f"S21 = {values[0]:.8f} V"
)
print(
f"S31 = {values[1]:.8f} V"
)
print(
f"S41 = {values[2]:.8f} V"
)
print(
f"S11 = {values[3]:.8f} V"
)
fig.canvas.draw_idle()
# ============================================================
# ACTIVATE MOUSE CALLBACK
# ============================================================
fig.canvas.mpl_connect(
"button_press_event",
onclick
)
# ============================================================
# DISPLAY
# ============================================================
plt.tight_layout()
plt.show()