From 860e34a068084d1865857b2ee3e48c2a871cd717 Mon Sep 17 00:00:00 2001 From: David Banas Date: Sun, 25 Feb 2024 08:03:51 -0500 Subject: [PATCH] Added windowing and explicit fmax/fstep. - Now applying a raised cosine filter to frequency domain channel descriptions, before `irfft()`ing them. - Made maximum frequency and frequency step explicit configuration parameters in the GUI. - Made it necessary to explicitly request automatic port renumbering for s4p files. - Removed some now unused functions from utility.py. --- src/pybert/gui/view.py | 23 +++-- src/pybert/models/bert.py | 33 ++++++-- src/pybert/pybert.py | 78 +++++++++++------ src/pybert/utility.py | 171 +++++++++++--------------------------- 4 files changed, 142 insertions(+), 163 deletions(-) diff --git a/src/pybert/gui/view.py b/src/pybert/gui/view.py index eb99bc6..4c2c4d6 100644 --- a/src/pybert/gui/view.py +++ b/src/pybert/gui/view.py @@ -110,6 +110,22 @@ label="Impulse Response Length (ns)", tooltip="Manual impulse response length override", ), + HGroup( + Item( + name="f_max", + label="Maximum Frequency", + tooltip="Max. frequency to use in generating H(f).", + ), + Item(label="GHz"), + ), + HGroup( + Item( + name="f_step", + label="Frequency Step", + tooltip="Frequency step to use in generating H(f).", + ), + Item(label="MHz"), + ), Item( name="thresh", label="Pj Threshold (sigma)", @@ -186,13 +202,6 @@ label="Use File", enabled_when="ch_file_valid == True", ), - Item( - name="f_step", - label="f_step", - tooltip="Frequency step to use in generating H(f).", - enabled_when="use_ch_file == True", - ), - Item(label="MHz"), ), HGroup( Item( diff --git a/src/pybert/models/bert.py b/src/pybert/models/bert.py index 1459b6d..2a1b109 100644 --- a/src/pybert/models/bert.py +++ b/src/pybert/models/bert.py @@ -27,6 +27,7 @@ mean, ones, pad, + pi, repeat, resize, std, @@ -37,6 +38,7 @@ from numpy.fft import fft, irfft from numpy.random import normal from scipy.signal import iirfilter, lfilter +from scipy.interpolate import interp1d from pybert.models.dfe import DFE from pybert.utility import ( @@ -47,6 +49,7 @@ import_channel, make_bathtub, make_ctle, + raised_cosine, safe_log10, trim_impulse, ) @@ -138,6 +141,7 @@ def _check_sim_status(): # Pull class variables into local storage, performing unit conversion where necessary. t = self.t + t_model = self.t_model w = self.w bits = self.bits symbols = self.symbols @@ -171,12 +175,17 @@ def _check_sim_status(): bandwidth = self.sum_bw * 1.0e9 rel_thresh = self.thresh mod_type = self.mod_type[0] + impulse_length = self.impulse_length try: # Calculate misc. values. fs = bit_rate * nspb Ts = t[1] ts = Ts + min_len = 30 * nspui + max_len = 100 * nspui + if impulse_length: + min_len = max_len = impulse_length / ts # Generate the ideal over-sampled signal. # @@ -291,7 +300,7 @@ def _check_sim_status(): self.status = "Tx GetWave() Error." self.log("ERROR: Never saw a rising step come out of Tx GetWave()!", alert=True) return - tx_h, _ = trim_impulse(diff(tx_s)) + tx_h, _ = trim_impulse(diff(tx_s), min_len=min_len, max_len=max_len) tx_out, _ = tx_model.getWave(self.x) else: # Init()-only. tx_out = convolve(tx_h, self.x) @@ -435,10 +444,20 @@ def _check_sim_status(): ctle_H = fft(ctle_h) ctle_H *= sum(ctle_h) / ctle_H[0] else: - _, ctle_H = make_ctle(rx_bw, peak_freq, peak_mag, w, ctle_mode, ctle_offset) - ctle_h = irfft(ctle_H) - ctle_h.resize(len(chnl_h), refcheck=False) - ctle_out = convolve(rx_in, ctle_h) + ctle_w, ctle_H = make_ctle(rx_bw, peak_freq, peak_mag, w, ctle_mode, ctle_offset) + ctle_h = irfft(raised_cosine(ctle_H)) + spl = interp1d(t_model, ctle_h, bounds_error=False, fill_value=0) + ctle_h = spl(t) + ctle_h *= abs(ctle_H[0]) / sum(ctle_h) + spl = interp1d(ctle_w/(2*pi), ctle_H, bounds_error=False, fill_value='extrapolate') + ctle_H = spl(self.f_plot) + ctle_h, _ = trim_impulse(ctle_h, front_porch=False, min_len=min_len, max_len=max_len) + try: + ctle_out = convolve(rx_in, ctle_h)[:len(rx_in)] + except: + print(f"rx_in: {rx_in}") + print(f"ctle_h: {ctle_h}") + raise ctle_out -= mean(ctle_out) # Force zero mean. if self.ctle_mode == "AGC": # Automatic gain control engaged? ctle_out *= 2.0 * decision_scaler / ctle_out.ptp() @@ -789,7 +808,7 @@ def update_results(self): eye_uis = self.eye_uis num_ui = self.nui clock_times = self.clock_times - f = self.f + f = self.f_plot t = self.t t_ns = self.t_ns t_ns_chnl = self.t_ns_chnl @@ -800,7 +819,7 @@ def update_results(self): ignore_samps = (num_ui - eye_uis) * samps_per_ui # Misc. - f_GHz = f[: len(f) // 2] / 1.0e9 + f_GHz = f / 1.0e9 len_f_GHz = len(f_GHz) len_t = len(t_ns) self.plotdata.set_data("f_GHz", f_GHz[1:]) diff --git a/src/pybert/pybert.py b/src/pybert/pybert.py index 4a826df..abd7d61 100755 --- a/src/pybert/pybert.py +++ b/src/pybert/pybert.py @@ -26,7 +26,7 @@ import skrf as rf from chaco.api import ArrayPlotData, GridPlotContainer from numpy import ( - append, array, convolve, cos, cumsum, diff, exp, linspace, logspace, + append, arange, array, convolve, cos, cumsum, diff, exp, linspace, logspace, log10, ones, pad, pi, roll, sinc, where, zeros) from numpy.fft import fft, irfft from numpy.random import randint @@ -65,10 +65,8 @@ calc_com, calc_gamma, calc_s2p_resp, - calc_sbr, com_opt, import_channel, - interp_s2p, interp_time, lfsr_bits, make_ctle, @@ -135,7 +133,8 @@ class PyBERT(HasTraits): do_xtalk = Bool(False) #: Include crosstalk, for s32p file? (Default = False) ch_is_s32p = Bool(False) #: Is channel Touchstone file 32-port? (Default = False) victim_chnl_ix = Range(low=1, high=8, value=1) #: Victim channel # when s32p given. (Default = 1) - f_step = Float(10) #: Frequency step to use when constructing H(f). (Default = 10 MHz) + f_step = Float(10) #: Frequency step to use when constructing H(f) (MHz). (Default = 10 MHz) + f_max = Float(40) #: Maximum frequency to use when constructing H(f) (GHz). (Default = 40 GHz) impulse_length = Float(0.0) #: Impulse response length. (Determined automatically, when 0.) Rdc = Float(0.1876) #: Channel d.c. resistance (Ohms/m). w0 = Float(10e6) #: Channel transition frequency (rads./s). @@ -402,6 +401,9 @@ class PyBERT(HasTraits): default_value="IEEE-802.3ck", ) com = Float(0) + com_Av = Float(0.4) + com_Afe = Float(0.4) + com_Ane = Float(0.6) com_Asig = Float(0) com_loc = Int(0) com_ser = Float(1e-4) @@ -490,8 +492,10 @@ class PyBERT(HasTraits): rel_opt = Property(Float, depends_on=["cost"]) t = Property(Array, depends_on=["ui", "nspb", "nbits"]) t_ns = Property(Array, depends_on=["t"]) - f = Property(Array, depends_on=["t"]) - w = Property(Array, depends_on=["f"]) + f_plot = Property(Array, depends_on=["t"]) + f_model = Property(Array, depends_on=["f_step", "f_max"]) + w = Property(Array, depends_on=["f_model"]) + t_model = Property(Array, depends_on=["f_model"]) bits = Property(Array, depends_on=["pattern", "nbits", "mod_type", "run_count"]) symbols = Property(Array, depends_on=["bits", "mod_type", "vod"]) ffe = Property(Array, depends_on=["tx_taps.value", "tx_taps.enabled"]) @@ -740,24 +744,41 @@ def _get_t_ns(self): return self.t * 1.0e9 @cached_property - def _get_f(self): - """Calculate the frequency vector appropriate for indexing non-shifted - FFT output, in Hz. - + def _get_f_plot(self): + """ + Calculate the frequency vector appropriate for indexing non-shiftedFFT output, in Hz. # (i.e. - [0, f0, 2 * f0, ... , fN] + [-(fN - f0), -(fN - 2 * f0), ... , -f0] - Note: Changed to positive freqs. only, in conjunction w/ irfft() usage. + Note: Changed to positive freqs. only, in conjunction w/ `irfft()` usage. """ t = self.t npts = len(t) f0 = 1.0 / (t[1] * npts) - half_npts = npts // 2 - return array([i * f0 for i in range(half_npts + 1)]) + return array([i * f0 for i in range(npts//2 + 1)]) + + @cached_property + def _get_f_model(self): + """ + Calculate the frequency vector for channel model construction. + """ + fstep = self.f_step * 1e6 + fmax = self.f_max * 1e9 + return arange(0, fmax+fstep, fstep) # "+fstep", so fmax gets included @cached_property def _get_w(self): """System frequency vector, in rads./sec.""" - return 2 * pi * self.f + return 2 * pi * self.f_model + + @cached_property + def _get_t_model(self): + """ + Calculate the equivalent time vector to `f_model`. + """ + f_model = self.f_model + tmax = 1 / f_model[1] + tstep = 0.5 / f_model[-1] + return arange(0, tmax, tstep) @cached_property def _get_bits(self): @@ -1563,8 +1584,7 @@ def calc_chnl_h(self): """ t = self.t - f = self.f - w = self.w + f = self.f_model nspui = self.nspui impulse_length = self.impulse_length * 1.0e-9 Rs = self.rs @@ -1581,7 +1601,8 @@ def calc_chnl_h(self): if self.use_ch_file: ch_s2p_pre, self.aggressors = import_channel( self.ch_file, ts, f, self.com_Av, self.com_Afe, self.com_Ane, - vic_chnl_ix=self.victim_chnl_ix) + vic_chnl_ix=self.victim_chnl_ix, + renumber=True) else: # Construct PyBERT default channel model (i.e. - Howard Johnson's UTP model). # - Grab model parameters from PyBERT instance. @@ -1612,7 +1633,9 @@ def add_ondie_s(s2p, ts4f, isRx=False): """ ts4N = rf.Network(ts4f) # Grab the 4-port single-ended on-die network. ntwk = sdd_21(ts4N) # Convert it to a differential, 2-port network. - ntwk2 = interp_s2p(ntwk, s2p.f) # Interpolate to system freqs. + # Interpolate to system freqs. + ntwk2 = ntwk.extrapolate_to_dc().windowed(normalize=False).interpolate( + s2p.f, coords='polar', bounds_error=False, fill_value='extrapolate') if isRx: res = s2p**ntwk2 else: # Tx @@ -1654,11 +1677,11 @@ def apply_Zt(s2p, _Zt, _min_len, _max_len, fstep=10e6, front_porch=True): spl = interp1d(chnl_t, chnl_h, bounds_error=False, fill_value=0) chnl_h = spl(t) chnl_dly = where(chnl_h == max(chnl_h))[0][0] * ts - chnl_h, start_ix = trim_impulse( + chnl_h_trim, start_ix = trim_impulse( chnl_h, min_len=_min_len, max_len=_max_len, front_porch=front_porch) - chnl_s = cumsum(chnl_h) - chnl_h *= abs(chnl_H[0]) / chnl_s[-1] - return chnl_h, chnl_f, chnl_H, start_ix, chnl_dly + chnl_s_trim = cumsum(chnl_h_trim) + chnl_h_trim *= abs(chnl_H[0]) / chnl_s_trim[-1] + return chnl_h_trim, chnl_f, chnl_H, start_ix, chnl_dly, chnl_h # Do it for PyBERT. min_len = 30 * nspui @@ -1671,13 +1694,16 @@ def Zt(f): """ return RL / (1 + 1j * 2*pi*f * RL * Cp) - chnl_h, chnl_f, chnl_H, start_ix, chnl_dly = apply_Zt( + chnl_h, chnl_f, chnl_H, start_ix, chnl_dly, chnl_h_orig = apply_Zt( ch_s2p, Zt, min_len, max_len, fstep=self.f_step*1e6) chnl_s = chnl_h.cumsum() chnl_p = chnl_s - pad(chnl_s[:-nspui], (nspui, 0), "constant", constant_values=(0, 0)) temp = chnl_h.copy() temp.resize(len(t), refcheck=False) - chnl_trimmed_H = fft(temp) + chnl_trimmed_H = fft(temp) # Has a different fundamental freq.! + chnl_trimmed_H *= abs(chnl_H[0]) / abs(chnl_trimmed_H[0]) + spl = interp1d(chnl_f, chnl_H, bounds_error=False, fill_value=0) + chnl_H = spl(self.f_plot) # Do it for COM, which takes care of modeling the Rx AFE itself. min_len = max_len = len(chnl_h) @@ -1687,7 +1713,7 @@ def Zt(f): """ return RL - com_chnl_h, com_chnl_f, com_chnl_H, _, _ = apply_Zt( + com_chnl_h, com_chnl_f, com_chnl_H, _, _, _ = apply_Zt( ch_s2p, Zt, min_len, max_len, fstep=self.f_step*1e6) com_chnl_s = com_chnl_h.cumsum() com_chnl_p = com_chnl_s - pad(com_chnl_s[:-nspui], (nspui, 0), "constant", constant_values=(0, 0)) @@ -1701,7 +1727,7 @@ def init_edge(x): return where(x_diff_abs > 0.2*x_diff_abs_max)[0][0] com_chnl_p_init = init_edge(com_chnl_p) - self.agg_hs, self.agg_fs, self.agg_Hs, _, _ = list( + self.agg_hs, self.agg_fs, self.agg_Hs, _, _, _ = list( zip(*[ apply_Zt(agg, Zt, min_len, max_len, fstep=self.f_step*1e6) for agg in self.aggressors ]) ) diff --git a/src/pybert/utility.py b/src/pybert/utility.py index bfed362..f381dee 100644 --- a/src/pybert/utility.py +++ b/src/pybert/utility.py @@ -343,6 +343,8 @@ def calc_jitter( t_jitter.append(ideal_xing) jitter = array(jitter) + assert len(jitter) > 0 and len(t_jitter) > 0, "No crossings found!" + if debug: print("mean(jitter):", mean(jitter)) print("len(jitter):", len(jitter)) @@ -386,7 +388,12 @@ def calc_jitter( # -- Make TIE vector uniformly sampled in time, via interpolation, for use as input to `fft()`. # spl = UnivariateSpline(t_jitter, jitter) # Way of the future, but does funny things. :( - spl = interp1d(t_jitter, jitter, bounds_error=False, fill_value="extrapolate") + try: + spl = interp1d(t_jitter, jitter, bounds_error=False, fill_value="extrapolate") + except: + print(f"t_jitter: {t_jitter}") + print(f"jitter: {jitter}") + raise tie_interp = spl(t) y = fft(tie_interp) jitter_spectrum = abs(y[:half_len]) @@ -949,7 +956,7 @@ def H_2_s2p(H, Zc, fs, Zref=50): return rf.Network(s=tmp, f=fs / 1e9, z0=[Zref, Zref]) # `f` is presumed to have units: GHz. -def import_channel(filename, sample_per, fs, Av, Afe, Ane, zref=100, vic_chnl_ix=1): +def import_channel(filename, sample_per, fs, Av, Afe, Ane, zref=100, vic_chnl_ix=1, renumber=False): """Read in a channel description file. Args: @@ -958,8 +965,12 @@ def import_channel(filename, sample_per, fs, Av, Afe, Ane, zref=100, vic_chnl_ix fs([real]): (Positive only) frequency values being used by caller (Hz). KeywordArgs: - zref(real): Reference impedance (Ohms), for time domain files. (Default = 100) - vic_chnl_ix(int): Victim channel number (from 1). (Default = 1) + zref(real): Reference impedance (Ohms), for time domain files. + Default = 100 + vic_chnl_ix(int): Victim channel number (from 1). + Default = 1 + renumber(bool): Automatically detect/fix "1=>3/2=>4" port numbering, when True. + Default = False Returns: (skrf.Network, [skrf.Network]): Pair consisting of: @@ -984,23 +995,23 @@ def import_channel(filename, sample_per, fs, Av, Afe, Ane, zref=100, vic_chnl_ix port 2 of all returned networks corresponds to the same physical node, the victim Rx node. """ - assert vic_chnl_ix > 0 and vic_chnl_ix <= 8, f"Victim index ({vic_chnl_ix}) out of range!" aggs = [] extension = os.path.splitext(filename)[1][1:] tstone_ext = re.match(r"^s(\d+)p$", extension, re.ASCII | re.IGNORECASE) if tstone_ext: # Touchstone file? n_ports = int(tstone_ext.group(1)) if n_ports == 32: + assert vic_chnl_ix > 0 and vic_chnl_ix <= 8, f"Victim index ({vic_chnl_ix}) out of range!" ntwks = import_s32p(filename, Av, Afe, Ane, vic_chnl_ix) assert ntwks[0].f[-1] < 1e12, f"Maximum frequency > 1 THz!\n\tfs = {fs}\n\tntwks[0] = {ntwks[0]}" - # chnls = list(map(lambda ntwk: interp_s2p(ntwk, fs), ntwks)) - # ts2N = chnls[0] ts2N = ntwks[0] assert ts2N.f[-1] < 1e12, f"Maximum frequency > 1 THz!\n\tfs = {fs}\n\tts2N = {ts2N}" - # aggs = chnls[1:] aggs = ntwks[1:] else: - ts2N = interp_s2p(import_freq(filename), fs) + ts2N = import_freq( + filename, renumber=renumber + ).extrapolate_to_dc().windowed(normalize=False).interpolate( + fs, coords='polar', bounds_error=False, fill_value='extrapolate') else: # simple 2-column time domain description (impulse or step). h = import_time(filename, sample_per) # Fixme: an a.c. coupled channel breaks this naive approach! @@ -1122,7 +1133,6 @@ def se2mm(ntwk, norm=0.5, renumber=False): assert rs == 4, "Touchstone file must have 4 ports!" # Detect/correct "1 => 3" port numbering if requested. - # ix = ntwk.s.shape[0] // 5 # So as not to be fooled by d.c. blocking. if renumber: ix = 1 if abs(ntwk.s21.s[ix, 0, 0]) < abs(ntwk.s31.s[ix, 0, 0]): # 1 ==> 3 port numbering? @@ -1158,13 +1168,17 @@ def se2mm(ntwk, norm=0.5, renumber=False): return rf.Network(frequency=f, s=s, z0=z) -def import_freq(filename): +def import_freq(filename, renumber=False): """Read in a 1, 2, or 4-port Touchstone file, and return an equivalent 2-port network. Args: filename(str): Name of Touchstone file to read in. + KeywordArgs: + renumber(bool): Automatically detect/fix "1=>3/2=>4" port numbering, when True. + Default = False + Returns: skrf.Network: 2-port network. @@ -1183,7 +1197,7 @@ def import_freq(filename): # Convert to a 2-port network. if rs == 4: # 4-port Touchstone files are assumed single-ended! - return sdd_21(ntwk) + return sdd_21(ntwk, renumber=renumber) if rs == 2: return ntwk return rf.network.one_port_2_two_port(ntwk) @@ -1383,97 +1397,6 @@ def mon_mag(zs): return zs_flat.reshape(zs.shape) -def interp_s2p(ntwk, f): - """Safely interpolate a 2-port network, by applying certain constraints to - any necessary extrapolation. - - Args: - ntwk(skrf.Network): The 2-port network to be interpolated. - f([real]): The list of new frequency sampling points (Hz). - - Returns: - skrf.Network: The interpolated/extrapolated 2-port network. - - Raises: - ValueError: If `ntwk` is _not_ a 2-port network. - """ - (fs, rs, cs) = ntwk.s.shape - assert rs == cs, "Non-square Touchstone file S-matrix!" - assert rs == 2, "Touchstone file must have 2 ports!" - - extrap = ntwk.interpolate( - # f, fill_value="extrapolate", coords="polar", assume_sorted=True) - f, coords="polar", assume_sorted=True, kind="rational") - assert extrap.f[-1] < 1e12, f"Maximum frequency > 1 THz!\n\tf: {f}\n\textrap: {extrap}" - s11 = cap_mag(extrap.s[:, 0, 0]) - s22 = cap_mag(extrap.s[:, 1, 1]) - # s12 = ntwk.s12.interpolate(f, fill_value=0, bounds_error=False, coords="polar", assume_sorted=True).s.flatten() - # s21 = ntwk.s21.interpolate(f, fill_value=0, bounds_error=False, coords="polar", assume_sorted=True).s.flatten() - s12 = cap_mag(extrap.s[:, 0, 1]) - s21 = cap_mag(extrap.s[:, 1, 0]) - s = np.array(list(zip(zip(s11, s12), zip(s21, s22)))) - if ntwk.name is None: - ntwk.name = "s2p" - return rf.Network(f=f, s=s, z0=extrap.z0, name=(ntwk.name + "_interp"), f_unit="Hz") - - -# ToDo: Are there any uses of this function remaining? Can we eliminate them? -def renorm_s2p(ntwk, zs): - """Renormalize a simple 2-port network to a new set of port impedances. - - This function was originally written as a check on the - `skrf.Network.renormalize()` function, which I was attempting to use - to model the Rx termination when calculating the channel impulse - response. (See lines 1640-1650'ish of `pybert.py`.) - - In my original specific case, I was attempting to model an open - circuit termination. And when I did the magnitude of my resultant - S21 dropped from 0 to -44 dB! - I didn't think that could possibly be correct. - So, I wrote this function as a check on that. - - Args: - ntwk(skrf.Network): A 2-port network, which must use the same - (singular) impedance at both ports. - - zs(complex array-like): The set of new port impedances to be - used. This set of frequencies may be unique for each port and at - each frequency. - - Returns: - skrf.Network: The renormalized 2-port network. - """ - (Nf, Nr, Nc) = ntwk.s.shape - assert Nr == 2 and Nc == 2, "May only be used to renormalize a 2-port network!" - assert all(ntwk.z0[:, 0] == ntwk.z0[0, 0]) and all( - ntwk.z0[:, 0] == ntwk.z0[:, 1] - ), f"May only be used to renormalize a network with equal (singular) reference impedances! z0: {ntwk.z0}" - assert zs.shape == (2,) or zs.shape == ( - len(ntwk.f), - 2, - ), "The list of new impedances must have shape (2,) or (len(ntwk.f), 2)!" - - if zs.shape == (2,): - zt = zs.repeat(len(Nf)) - else: - zt = np.array(zs) - z0 = ntwk.z0[0, 0] - S = ntwk.s - I = np.identity(2) - Z = [] - for s in S: - Z.append(inv(I - s).dot(I + s)) # Resultant values are normalized to z0. - Z = np.array(Z) - Zn = [] - for z, zn in zip(Z, zt): # Iterration is over frequency and yields: (2x2 array, 2-element vector). - Zn.append(z.dot(z0 / zn)) - Zn = np.array(Zn) - Sn = [] - for z in Zn: - Sn.append(inv(z + I).dot(z - I)) - return rf.Network(s=Sn, f=ntwk.f / 1e9, z0=zs) - - def getwave_step_resp(ami_model): """Use a model's GetWave() function to extract its step response. @@ -1912,8 +1835,8 @@ def calc_sbr(s2p, ui, Zt=None, fstep=10e6): ui(real): Unit interval (s). KeywordArgs: - Zt(complex or [complex]): Termination impedance, either scalar or as function of frequency. - (Default: None) + Zt(real -> complex): Function from frequency (Hz) to response. + Default: None fstep(real): Frequency step to use. Returns: @@ -1947,31 +1870,25 @@ def calc_s2p_resp(s2p, Zt=None, fstep=10e6): h: impulse response, and H: frequency response. """ - s2p_term = s2p.copy() + # Ensure uniform frequency sampling down to d.c. + fmax = s2p.f[-1] + f = arange(0, fmax+fstep, fstep) # "+fstep" so fmax will actually be included. + s2p_term = s2p.extrapolate_to_dc().windowed(normalize=False).interpolate( + f, bounds_error=False, fill_value='extrapolate') + if Zt is not None: - s2p_term_z0 = s2p.z0.copy() - s2p_term_z0[:, 1] = array(list(map(Zt, s2p.f))) + new_z0 = s2p_term.z0.copy() + new_z0[:, 1] = array(list(map(Zt, f))) try: - s2p_term.renormalize(s2p_term_z0) + s2p_term.renormalize(new_z0) except: - print(s2p_term_z0) + print(new_z0) + print(f) raise - # Ensure uniform frequency sampling down to d.c. - fmax = s2p_term.f[-1] - f = arange(fstep, fmax+fstep, fstep) # "+fstep" so fmax will actually be included. - s2p_resamp = s2p_term.interpolate(f).extrapolate_to_dc() - f = s2p_resamp.f - - # ToDo: Use SciKit-RF `impulse_response()` and `step_response()`? - # We take the transfer function, H, to be a ratio of voltages. # So, we must normalize our (now generalized) S-parameters. - H = s2p_resamp.s21.s.flatten() * np.sqrt(s2p_resamp.z0[:, 1] / s2p_resamp.z0[:, 0]) - - # Apply windowing, to eliminate artifacts. - win = (cos(pi * f / fmax) + 1)/2 - Hwin = H * win + Hwin = s2p_term.s21.s.flatten() * np.sqrt(s2p_term.z0[:, 1] / s2p_term.z0[:, 0]) # Calculate impulse response and associated time vector. h = irfft(Hwin) @@ -1981,3 +1898,11 @@ def calc_s2p_resp(s2p, Zt=None, fstep=10e6): t = array([n*ts for n in range(len(h))]) return (t, h, f, Hwin) + + +def raised_cosine(x): + """Apply raised cosine window to input. + """ + len_x = len(x) + w = (array([cos(pi * n / len_x) for n in range(len_x)]) + 1) / 2 + return w * x