KNOW-HOW

Summary of Equalizations

32 to 56 Gbps Serial Link Analysis and Optimization Methods for Pathological Channels [https://docs.keysight.com/eesofapps/files/678068240/678068273/1/1629077956000/tutorial-32-to-56-gbps-serial-link-analysis-optimization-methods-pathological-channels.pdf]

Qasim Chaudhari. A Classification of Equalization Techniques [https://wirelesspi.com/a-classification-of-equalization-techniques/]

CTLE vs. FFE

Keysight Signal Integrity Educational Posts [Post 5: Root Cause of Eye Closure], [Post 6: Eye-opening Experience with CTLE], [Post 7: Eye-opening Experience with FFE]

in the time domain

• CTLE provide only limited improvement in the pre-cursor ISI, because of the continuous-time, analog nature of CTLE
• FFE to reduce ISI in both the pre-cursor and post-cursor, because of operating digitally in discrete-time

in the frequency domain

• CTLE focuses on boosting frequency content at the Nyquist frequency
• FFE algorithm is selecting taps that effectively amplify the odd harmonics of the Nyquist frequency

In the case of FFE, because of the nature of finite impulse response filter, we would expect amplification and attenuation of different harmonics of Nyquist Frequency

Until 6.5 dB of CTLE DC attenuation, the spread of the single pulse is positive and reaches almost zero at 6.5 dB. As the DC attenuation increases to more than 6.5 dB, the single pulse spectrum is restored too much, resulting in a negative dip at the end of the pulse

the maximum eye opening does not happen at maximum DC attenuation at 9 dB

DFE

Keysight Signal Integrity Educational Posts [Post 8: Eye-opening Experience with DFE]

There are kinks in the eye diagram, the signature of an opened DFE eye is different than other equalizations

DFE algorithm is reducing ISI based on the detected data (symbol)

Since DFE assumes that past symbol decisions are correct. Incorrect decisions from the symbol detector corrupt the filtering of the feedback loop. As a result, the inclusion of the feedforward filter on the front end is crucial in minimizing the probability of error

• because symbol detection is nonlinear, decision feedback equalization is also nonlinear
• because of the nonlinearity of the DFE response, it must be modeled in the time domain

FFE vs. DFE

FFE: convolution with input waveform/symbols. Linear.

DFE: convolution with past detected symbols. Recursive and nonlinear because of the slicer.

That is why DFE is not a simple LTI convolution system from input to output.

1
2
3
4
5

received sample r[n] ──┬── subtract ──> slicer ──> detected symbol a_hat[n]
                       ▲                    │
                       │                    │
                       └── DFE filter <─────┘
                          past decisions

FIR Coefficient Selection

Jose E. Schutt-Aine, Spring 2024 ECE 546 Lecture - 27 Equalization [http://emlab.uiuc.edu/ece546/Lect_27.pdf]

Sam Palermo. Lecture 7 - Equalization Intro & TX FIR EQ [https://people.engr.tamu.edu/spalermo/ecen689/lecture7_ee720_eq_intro_txeq.pdf]

Kevin Zheng, Boris Murmann, Hongtao Zhang, and Geoff Zhang. Feedforward Equalizer Location Study for High-Speed Serial Systems [https://www.signalintegrityjournal.com/articles/1228-feedforward-equalizer-location-study-for-high-speed-serial-systems]

—, "System-Driven Circuit Design for ADC-Based Wireline Data Links", Ph.D. Dissertation, Stanford University, 2018 [https://purl.stanford.edu/hw458fp0168]

Hanumolu, P. K., Wei, G. Y., & Moon, Y. K. (2005). Equalizers for high-speed serial links. International Journal of High Speed Electronics and Systems [https://people.engr.tamu.edu/spalermo/ecen689/hslink_eq_overview_hanumolu_jhses05.pdf]

with MMSE

Lecture 7: Equalization Introduction & TX FIR Eq [https://people.engr.tamu.edu/spalermo/ecen689/lecture7_ee720_eq_intro_txeq.pdf]

Toeplitz matrix: transforms discrete convolution into \(y=Ax\), where \(x\) is a flattened input vector

Lone-Pulse Equalization

\[\begin{align} E^TE &=(W^T H^T - Y_{des}^T)(HW-Y_{des})=W^TH^THW+Y_{des}^TY_{des}-W^TH^TY_{des}-Y_{des}^THW \\ &=W^TH^THW+Y_{des}^TY_{des}-2Y_{des}^THW \end{align}\]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

h=[0.004, 0.0010, 0.0023, 0.0052, 0.0812, 0.3437, 0.1775, 0.0917, 0.0526,...
    0.0360, 0.0224, 0.0162, 0.0152, 0.0097, 0.0090, 0.0067];
k = length(h);
n = 3;
l = 1;
m2 = 5;
m1 = 1;

H = zeros([k+n+l-2, n+l-1]);
H(1:end-2,1) = h;
H(2:end-1,2) = h;
H(3:end,3) = h;

Ydes = zeros([k+n+l-2, 1]);
Ydes(m1+m2+1,1) = 1;

HT = transpose(H);

Wls = inv(HT*H)*HT*Ydes;

% Wls =
% 
%    -0.8177
%     3.7239
%    -1.7181

Wlsnorm = Wls/sum(norm(Wls,1));

% Wlsnorm =
% 
%    -0.1306
%     0.5949
%    -0.2745

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

fcsvf = readtable("hsample_pre10post20.csv");
h= fcsvf.hsample_Design_Point_1_Y;
k = length(h);
n = 8;
l = 1;
m2 = 10;  % channel pre-cursor sample#
m1 = 1;

H = zeros([k+n+l-2, n+l-1]);
for i =1:n
    H(i:i+k-1,i) = h;
end

Ydes = zeros([k+n+l-2, 1]);
Ydes(m1+m2+1,1) = 1;

HT = transpose(H);
Wls = inv(HT*H)*HT*Ydes;


Wlsnorm = Wls/sum(norm(Wls,1));
% Wlsnorm =
% 
%    -0.0926
%     0.6383
%    -0.2691

with ZFS

Zero Forcing Solution (ZFS)

The number of channel samples may exceed the number of equalizer taps to accurately compute the optimal tap coefficients

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

ht = [0.3, 1.0, -0.2, 0.1, 0.0, 0.0];
ytarget = [0;1;0];

x1 = [[1.0 0.3 0.0];
    [-0.2 1.0 0.3];
    [0.1 -0.2 1.0]];  % better

x2 = [[1.0 0.3 0.0];
    [-0.2 1.0 0.3];
    [0.0 -0.2 1.0]];

p1 = inv(x1)*ytarget; % -0.2657    0.8857    0.2037
p2 = inv(x2)*ytarget; % -0.2679    0.8929    0.1786

ht_p1 = conv(ht, p1); % p1 better -> x1
ht_p2 = conv(ht, p2);

% ht_p1 =
%
%    -0.0797         0    1.0000         0    0.0478    0.0204         0         0
%
% ht_p2 =
%
%    -0.0804    0.0000    1.0000   -0.0268    0.0536    0.0179         0         0

subplot(3,1,1)
stem(ht, 'LineWidth', 2); grid on; xlim([0,10])
subplot(3,1,2)
stem(p1, 'LineWidth', 2); hold on; stem(p2, 'LineWidth', 2);
grid on; legend(["p1" "p2"]); xlim([0,10])
subplot(3,1,3)
stem(ht_p1, 'LineWidth', 2); hold on; stem(ht_p2, 'LineWidth', 2)
grid on; legend(["h\_p1" "h\_p2"]); xlim([0,10])

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

h = [0.01 -0.02 0.05 -0.1 0.2 1 0.15 -0.15 0.05 -0.02 0.005];
[val, idx] = max(h);

htc = h(idx-2:idx+2)';

H1 = [1 0.2 -0.1 0.05 -0.02;
    0.15 1 0.2 -0.1 0.05;
    -0.15 0.15 1 0.2 -0.1;
    0.05 -0.15 0.15 1 0.2;
    -0.02 0.05 -0.15 0.15 1];  % better

H2 = [1 0.2 -0.1 0 0;
    0.15 1 0.2 -0.1 0;
    -0.15 0.15 1 0.2 -0.1;
    0 -0.15 0.15 1 0.2;
    0 0 -0.15 0.15 1];


ytgt = zeros(5,1);
ytgt(3) = 1;

heq1 = inv(H1)*ytgt;
heq2 = inv(H2)*ytgt;

ZFS vs MMSE

minimum mean squared error (MMSE)

There are three major MMSE-based algorithms:

• least mean square (LMS),
• normalized least mean square (NLMS)
• recursive least square (RLS)

• ZFS eliminates the ISI only at the sampling points that correspond to the equalizer taps. The equalized pulse shows ISI in the intervals between the sample points and at sample points outside the equalizer
• The Minimum Mean-Square Error Linear Equalizer (MMSE-LE) balances ISI reduction and noise enhancement. The MMSE-LE always performs as well as, or better than, the ZFE

LMS (Least-Mean-Square)

Qasim Chaudhari. Least Mean Square (LMS) Equalizer – A Tutorial [https://wirelesspi.com/least-mean-square-lms-equalizer-a-tutorial/]

CC Chen, Why Background EQ Adaptation? [https://youtu.be/l46OesuNfp4]

TX with SS-LMS

V. Stojanovic et al., "Autonomous dual-mode (PAM2/4) serial link transceiver with adaptive equalization and data recovery," IEEE Journal of Solid-State Circuits, vol. 40, no. 4, pp. 1012–1026, Apr. 2005 [https://sci-hub.ru/10.1109/JSSC.2004.842863]

—, "Channel-Limited High-Speed Links: Modeling, Analysis and Design," PhD. Thesis, Stanford University, Sep. 2004. [pdf]

—, US7423454B2. High speed signaling system with adaptive transmit pre-emphasis [pdf]

\[ dLev_{n+1} = dLev_n - \frac{\Delta_{dLev}}{2}\left(\frac{\partial e_n^2}{\partial dLev_n}\right) = dLev_n - \Delta _{dLev} e_n\left(\frac{\partial (dLev_n-y_n)}{\partial dLev_n}\right) = \color{red} dLev_n - \Delta _{dLev} e_n \] note \(e_n = dLev_n-y_n\)

J. T. Stonick, Gu-Yeon Wei, J. L. Sonntag and D. K. Weinlader, "An adaptive PAM-4 5-Gb/s backplane transceiver in 0.25-μm CMOS," in IEEE Journal of Solid-State Circuits, vol. 38, no. 3, pp. 436-443, March 2003, [https://sci-hub.st/10.1109/JSSC.2002.808282]

RX with SS-LMS

E. -H. Chen et al., "Near-Optimal Equalizer and Timing Adaptation for I/O Links Using a BER-Based Metric," in IEEE Journal of Solid-State Circuits, vol. 43, no. 9, pp. 2144-2156, Sept. 2008 [https://sci-hub.ru/10.1109/JSSC.2008.2001871]

Sam Palermo. ECEN720: High-Speed Links Circuits and Systems [Lecture 7 - Equalization Intro & TX FIR EQ], [Lecture 8 - RX FIR, CTLE, DFE, & Adaptive Eq.]

Jinhyung Lee, Design of High-Speed Receiver for Video Interface with Adaptive Equalization; Phd thesis, August 2019. [thesis link]

Kwangho Lee, Design of Receiver with Offset Cancellation of Adaptive Equalizer and Multi-Level Baud-Rate Phase Detector; Phd thesis, August 2021.[pdf]

\(e[n] = d[n] - Dlev_n\cdot tx[n]\)

Bang-Bang CDR

Alexander PD or !!PD

By definition the edge sample will be zero at a zero crossing, given \(a_na_{n+1}=-1\)

By proper equalization choice, the pulse response may approximate even symmetry

\[ f(t) = g(t+T/2) - g(t-T/2) \]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

# https://share.google/aimode/l0gnYPTxlyUa7WUed

import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(-1.5, 1.5, 1000)

# Pulse Response and its continuous derivative
p = np.exp(-4 * t ** 2)
p_prime = -8 * t * np.exp(-4 * t ** 2)

# Discrete Approximation (Finite Difference over T=1 UI)
T = 1
discrete_approx = (np.exp(-4 * (t + T / 2) ** 2) - np.exp(-4 * (t - T / 2) ** 2)) / T

# Alexander Bang-Bang Output
g_tau = np.sign(discrete_approx)

plt.figure(figsize=(10, 6))
plt.plot(t, p, label="Pulse Response", color='k', alpha=0.8)
plt.plot(t, p_prime, label="Continuous Derivative $P'(t)$", color='red', alpha=0.6)
plt.plot(t, discrete_approx, label="Discrete Finite Difference Approx", color='blue', ls='--')
plt.step(t, g_tau, where='mid', color='green', lw=1, label="Alexander PD Timing Function", alpha=0.6)
plt.axhline(0, color='black', lw=1)
plt.axvline(0, color='gray', ls=':', lw=1, label='Lock Point')
plt.title('Discrete Approximation of Pulse Derivative (Page 42)')
plt.xlabel('Phase Error (UI)')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

Alexander (Bang-Bang) PD does not typically lock at the maximum pulse value when the pulse is asymmetric.

For an asymmetric pulse (like one with a slow trailing edge caused by ISI), this lock point shifts toward the slower-decaying side of the pulse.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

# https://share.google/aimode/l0gnYPTxlyUa7WUed

import numpy as np
import matplotlib.pyplot as plt


t = np.linspace(-1, 1.5, 1000)
T = 1  # Sampling width

# Asymmetric Pulse: Fast rise, slow tail
p = 0.8 * np.exp(-15 * (t + 0.1) ** 2) + 0.4 * np.exp(-1.5 * (t - 0.4) ** 2)

# Calculate Slope Approximation (Alexander PD Function)
p_plus = np.interp(t + T / 2, t, p)
p_minus = np.interp(t - T / 2, t, p)
slope = (p_plus - p_minus) / T

# Find Critical Points
t_peak = t[np.argmax(p)]
t_lock = t[np.argmin(np.abs(slope))]

plt.figure(figsize=(10, 5))
plt.plot(t, p, label="Asymmetric Pulse", lw=2)
plt.plot(t, slope, label="Discrete Finite Difference Approx", color='blue', ls='--')
plt.axvline(t_peak, color='red', ls='--', label=f'True Peak: {t_peak:.2f} UI')
plt.axvline(t_lock, color='green', ls='-', label=f'PD Lock Point: {t_lock:.2f} UI')
plt.title("Sampling Error in Alexander PD due to Pulse Asymmetry")
plt.legend(); plt.grid(True, alpha=0.3)
plt.show()

Kwangho Lee, "Design of Receiver with Offset Cancellation of Adaptive Equalizer and Multi-Level Baud-Rate Phase Detector" [https://s-space.snu.ac.kr/bitstream/10371/177584/1/000000167211.pdf]

Shahramian, Shayan, "Adaptive Decision Feedback Equalization With Continuous-time Infinite Impulse Response Filters" [https://tspace.library.utoronto.ca/bitstream/1807/77861/3/Shahramian_Shayan_201606_PhD_thesis.pdf]

MENIN, DAVIDE, "Modelling and Design of High-Speed Wireline Transceivers with Fully-Adaptive Equalization" [https://air.uniud.it/retrieve/e27ce0ca-15f7-055e-e053-6605fe0a7873/Modelling%20and%20Design%20of%20High-Speed%20Wireline%20Transceivers%20with%20Fully-Adaptive%20Equalization.pdf]

Oh, Kyung Suk, and Xingchao Yuan. High Speed Signaling Jitter Modeling, Analysis, and Budgeting. Pearson Education, 2012

Mueller-Muller CDR

Faisal A. Musa. "HIGH-SPEED BAUD-RATE CLOCK RECOVERY" [https://www.eecg.utoronto.ca/~tcc/thesis-musa-final.pdf]

—."CLOCK RECOVERY IN HIGH-SPEED MULTILEVEL SERIAL LINKS" [https://www.eecg.utoronto.ca/~tcc/faisal_iscas03.pdf]

K. Yadav, P. -H. Hsieh and A. C. Carusone, "Loop Dynamics Analysis of PAM-4 Mueller–Muller Clock and Data Recovery System," in IEEE Open Journal of Circuits and Systems, vol. 3, pp. 216-227, 2022 [https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=9910561]

Jaeduk Han, "Design and Automatic Generation of 60Gb/s Wireline Transceivers" [https://www2.eecs.berkeley.edu/Pubs/TechRpts/2019/EECS-2019-143.pdf]

S. Kim, K. K. Tokgoz and G. Kim, "Modeling and Simulation of Mueller-Muller Clock Data Recovery System for PAM-4 Wireline Transceivers," 2025 IEEE/IEIE International Conference on Consumer Electronics-Asia (ICCE-Asia), Busan, Korea, Republic of, 2025, pp. 1-3, doi: 10.1109/ICCE-Asia67487.2025.11263607

Mueller-Muller type A timing function

Mueller-Muller type B timing function

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35

# https://share.google/aimode/ajIRVJNOatjPnY2zp

import numpy as np
import matplotlib.pyplot as plt

# 1. Asymmetric Pulse Modeling
t_cont = np.linspace(-3, 8, 2000)
h_cont = np.where(t_cont >= -0.8, (t_cont + 0.8) * np.exp(-0.6 * (t_cont + 0.8)), 0)
h_cont /= np.max(h_cont)

def get_h(phase):
    return np.interp(phase, t_cont, h_cont)

# 2. Timing Error Calculations
phases = np.linspace(-0.5, 3.5, 1000)
tef_a = np.array([get_h(ts + 1) - get_h(ts - 1) for ts in phases])
tef_b = np.array([get_h(ts - 1) for ts in phases])

# 3. Robust Lock Logic
t_lock_a = phases[np.argmin(np.abs(tef_a))]

# Type B: Finding the LAST minimal error before h(-1) rises
indices_b = np.where(np.abs(tef_b) == np.min(np.abs(tef_b)))[0]
t_lock_b = phases[indices_b][-1]

# 4. Plotting Results
plt.figure(figsize=(10, 6))
plt.plot(t_cont, h_cont, 'k-', lw=2, label='Pulse Response $h(t)$')
plt.axvline(t_lock_a, color='red', linestyle='--', label=f'Type A Lock: {t_lock_a:.2f}T')
plt.axvline(t_lock_b, color='green', linestyle='--', label=f'Type B Lock: {t_lock_b:.2f}T')
plt.plot(t_lock_b - 1, get_h(t_lock_b - 1), 'go', label='Type B: $h_{-1}=0$')
plt.title('Baud Rate Lock: Type A (Symmetry) vs. Type B (Last-Zero Precursor)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()

Suppose 1-precursor, 1-postcursor — \(y_k = d_{k-1}h_1 + d_k h_0 + d_{k+1}h_{-1}\) \[ \color{red}E[y_k\cdot d_{k-1}] - E[y_k\cdot d_{k+1}] = E[|d_{k-1}|^2h_{1}] - E[|d_{k+1}|^2h_{-1}] =h_1-h_{-1} \] MMPD infers the channel response from baud-rate samples of the received data, the adaptation aligns the sampling clock such that pre-cursor is equal to the post-cursor in the pulse response

note \(E[y_k\cdot d_{k+1}] = E[y_{k-1}\cdot d_{k}] = h_{-1}\)

SS-MMPD

F. Spagna et al., "A 78mW 11.8Gb/s serial link transceiver with adaptive RX equalization and baud-rate CDR in 32nm CMOS," 2010 IEEE International Solid-State Circuits Conference - (ISSCC), San Francisco, CA, USA, 2010, [https://sci-hub.ru/10.1109/ISSCC.2010.5433823]

[https://people.engr.tamu.edu/spalermo/ecen689/lecture12_ee720_cdrs.pdf]

Suppose \(x_k = d_{k-1}h_1 + d_k h_0 + d_{k+1}h_{-1}\) and \(x_{k-1} = d_{k-2}h_1 + d_{k-1} h_0 + d_{k}h_{-1}\) \[ \color{red}E\{z_k\} = \frac{1}{2} E\{|d_{k-1}|^2h_1\} - \frac{1}{2} E\{|d_{k}|^2h_{-1}\} = \frac{1}{2}(h_1 - h_{-1}) \]

Reference Level Choice

Chen, J., Gu, Y., Feng, X., Chi, R., Wu, J., & Chen, Y. (2024). Analysis of Mueller–Muller Clock and Data Recovery Circuits with a Linearized Model. Electronics, 13(21), 4218 [https://www.mdpi.com/2079-9292/13/21/4218]

Liu, Tao & Li, Tiejun & Lv, Fangxu & Liang, Bin & Zheng, Xuqiang & Wang, Heming & Wu, Miaomiao & Lu, Dechao & Zhao, Feng. (2021). Analysis and Modeling of Mueller-Muller Clock and Data Recovery Circuits. Electronics. [10. 1888. 10.3390/electronics10161888.]

Gu, Youzhi & Feng, Xinjie & Chi, Runze & Chen, Yongzhen & Wu, Jiangfeng. (2022). Analysis of Mueller-Muller Clock and Data Recovery Circuits with a Linearized Model. [10.21203/rs.3.rs-1817774/v1]

Adjust Locking Point

Avago Technologies, US8649476B2 Adjusting sampling phase in a baud-rate CDR using timing skew [pdf]

Y. Jung, H. -J. Shin, J. Kim, S. Lee, J. -S. Park and K. Park, "A 28-Gb/s Receiver with Baud-Rate CDR Employing Integrated Pattern-Based Phase Detector Achieving ISI Invariant Phase Locking," 2025 IEEE Asian Solid-State Circuits Conference (A-SSCC), Daejeon, Korea, Republic of, 2025,

R. Dokania et al., "10.5 A 5.9pJ/b 10Gb/s serial link with unequalized MM-CDR in 14nm tri-gate CMOS," 2015 IEEE International Solid-State Circuits Conference - (ISSCC) Digest of Technical Papers, San Francisco, CA, USA, 2015 [https://sci-hub.jp/10.1109/ISSCC.2015.7062987]

H. Zhang, B. Jiao, Y. Liao, and G. Zhang, "A Tutorial on PAM4 Signaling for 56G Serial Link," [DesignCon 2016], [DesignCon 2017]

TODO 📅

Kwangho Lee, "Design of Receiver with Offset Cancellation of Adaptive Equalizer and Multi-Level Baud-Rate Phase Detector" [https://s-space.snu.ac.kr/bitstream/10371/177584/1/000000167211.pdf]

\(h_1\) is necessary

• without DFE

  SS-MMPD locks at the point (\(h_1=h_{-1}\)​)

• With a 1-tap DFE

  1-tap adaptive DFE that forces the \(h_1\) to be zero, the SS-MMPD locks wherever the \(h_{-1}\)​ is zero and drifts eventually.

  Consequently, it suffers from a severe multiple-locking problem with an adaptive DFE

Pattern filter

During adapting, we make

• \(s_{011}\) & \(s_{110}\) are approaching to each other
• \(s_{100}\) & \(s_{001}\) are approaching to each other

Then, \(h_{-1}\) and \(h_1\) are same, which is desired

CDR Loop Latency

A. Amirkhany, "Tutorial: Basics of Clock-and-Data Recovery Circuits," 2019 IEEE International Solid-State Circuits Conference - (ISSCC), San Francisco, CA, USA, 2019, pp. 1-105, doi: 10.1109/ISSCC19946.2019.11005770.

Nhat Nguyen and Masum Hossain, ISSCC 2021 Forum 6.7: 112Gb/s-and-Beyond Long-Reach and Short-Reach Electrical Interfaces

For high-BW CDR \(\omega_{ug}\), latency \(D\) constrain \(\phi_M\)

loop latency is represented as \(\color{red}e^{-sD}\) in linear model

Sensitivity to Loop Latency

Loop Latency model

CC Chen. Why A Low Loop Latency in A CDR Design? [https://youtu.be/io9WZbhlahU]

—. Why Understanding and Optimizing Loop Latency for A CDR Design? [https://youtu.be/Jyy18865jv8]

Walker, Richard. (2003). Designing Bang-Bang PLLs for Clock and Data Recovery in Serial Data Transmission Systems. [paper,slides]

\[ q[n]\rightarrow \left( K_P+\frac{K_I}{1-z^{-1}} \right) \rightarrow \frac{1}{1-z^{-1}} \rightarrow\phi_{CK} \]

The last term is already the VCO integration from frequency to phase.

So the two paths have different effects:

\[ K_P q[n] \quad\stackrel{\mathrm{VCO}}{\longrightarrow}\quad \text{one integration} \]

whereas

\[ K_I\sum q[n] \quad\stackrel{\mathrm{VCO}}{\longrightarrow}\quad \text{two integrations} \]

Using the delayed BBPD output

\[ q_D[n]=q[n-T_D],\qquad q_D[n]\in\{-1,+1\} \]

the loop filter is

\[ \boxed{ f_I[n+1]=f_I[n]+K_I q_D[n] } \]

and

\[ \boxed{ \Delta f_{\rm DCO}[n] = K_P q_D[n]+f_I[n] } \]

The DCO phase then evolves as

\[ \boxed{ \phi_{\rm DCO}[n+1] = \phi_{\rm DCO}[n] + \Delta\phi_{\rm nom} + \Delta f_{\rm DCO}[n] } \]

If we remove the nominal \(2\pi\) rotation and only track phase error, this becomes

\[ \boxed{ \Delta\phi_{\rm DCO}[n+1] = \Delta\phi_{\rm DCO}[n] + K_Pq_D[n]+f_I[n] } \]

So the structure is

\[ q_D \rightarrow \boxed{K_P+\frac{K_I}{1-z^{-1}}} \rightarrow \boxed{\frac{1}{1-z^{-1}}} \rightarrow \phi_{\rm DCO} \]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82

"""
Behavioral model of the experiment on the slides
"Max dphi w/ SJ & Sweeping loop Latency" (pages 5 and 6):

    8 Gb/s, PRBS7, no input RJ
    SJ = 10 %UI p-p @ 100 MHz
    T_D (loop latency) = 2 / 5 / 10 UI  ->  dphi = 12 / 19 / 23 %UI p-p

Units: time step = 1 UI (125 ps @ 8 Gb/s), all phases in UI.

Loop (digital type-II, latency T_D inserted in the feedback path):

    e[n]       = sign(dphi[n]) on data transitions, 0 otherwise   (BBPD)
    acc[n+1]   = acc[n]   + Ki * e[n-TD]                          (integral / freq)
    ph_ck[n+1] = ph_ck[n] + acc[n+1] + Kp * e[n-TD]               (VCO phase)
    dphi[n]    = ph_din[n] - ph_ck[n]

Fitted gains that reproduce the three measured points:

    Kp = 0.60 %UI per decision,  Ki = Kp/1024
    -> dphi(p-p) = 12.8 / 19.3 / 23.4 %UI   (measured 12 / 19 / 23)

Key mechanism (confirmed by the page-6 zoom): the loop is SLEW-RATE LIMITED,
not bandwidth limited.  The CDR can move clock phase at most

    SR    = Kp * eta                    ~ 3.0 mUI/UI  (eta = PRBS7 edge density)

while the SJ demands

    SJ_SR = pi * A_pp * f_sj / f_b      ~ 3.9 mUI/UI

Since SR < SJ_SR, ph_ck ramps at constant slope, falls behind, and then keeps
ramping past the input after the sine turns around; T_D extends that turnaround.
The result is a coherent triangular oscillation at f_sj whose amplitude EXCEEDS
the input jitter.  The regime boundary is

    f_SL  = Kp * eta * f_b / (pi * A_pp) ~ 77 MHz
"""

def bbcdr(TD_UI=2, Kp=KP_FIT, zeta_div=ZD_FIT, ki_en=None,
          sj_pp=0.10, f_sj=100e6, f_b=F_B, ppm=0.0,
          n_ui=20000, n_settle=4000, seed=0x7F):
    """
    One CDR run.  Returns the waveforms plus p-p / peak / rms phase error
    measured after n_settle UI.

    dphi = ph_din - ph_ck   (>0 : data late -> clock must be delayed)

    ki_en : True -> type-II (Ki = Kp/zeta_div), False -> type-I (Ki = 0),
            None -> follow the module-level KI_EN.
    """
    use_ki = KI_EN if ki_en is None else ki_en
    Ki = Kp / zeta_div if use_ki else 0.0
    TD = int(round(TD_UI))

    n = np.arange(n_ui)
    ph_din = 0.5 * sj_pp * np.sin(2 * np.pi * (f_sj / f_b) * n) + ppm * 1e-6 * n

    bits = prbs7(n_ui + 1, seed)
    trans = (bits[1:] != bits[:-1]).astype(np.int8)   # BBPD updates on edges only

    ph_ck = np.empty(n_ui)
    dphi = np.empty(n_ui)
    e_buf = np.zeros(TD + 1)                          # latency pipeline
    acc = 0.0
    ck = 0.0

    for i in range(n_ui):
        d = ph_din[i] - ck
        dphi[i] = d
        ph_ck[i] = ck

        e_buf[1:] = e_buf[:-1]
        e_buf[0] = np.sign(d) * trans[i]              # BBPD decision
        e_del = e_buf[TD]                             # delayed by T_D UI

        acc += Ki * e_del
        ck += acc + Kp * e_del

    d = dphi[n_settle:]
    return dict(dphi=dphi, ph_ck=ph_ck, ph_din=ph_din, t=n,
                pp=d.max() - d.min(), peak=np.abs(d).max(), rms=d.std())

\(K_P\) is the normalized DCO frequency deviation caused by one BBPD decision, and because that deviation lasts for one UI, it produces \(K_P\) UI of excess phase over that interval:

\[ \Delta\phi_{\rm excess} = \frac{\Delta f_{\rm DCO}}{f_{\rm data}} = K_P \qquad \boxed{\text{unit}: \space\mathrm{UI/UI}} \]

another model

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72

def simulate(td_ui: int, kp: float, ki: float):
    """
    BB-CDR recursion.

    e[n]    = phi_data[n] - phi_clk[n]

    q[n]    = sign(e[n]) when a PRBS data transition exists.
              During no-transition intervals, the binary detector keeps
              the previous Early/Late state.

    qd[n]   = q[n-TD]

    fi[n+1] = fi[n] + Ki*qd[n]

    df[n]   = Kp*qd[n] + fi[n+1]

    phi_clk[n+1] = phi_clk[n] + df[n]

    The last line is the VCO frequency-to-phase integration.
    """
    n = np.arange(N_UI)

    # 100 MHz SJ at 8 Gb/s => period = 80 UI
    phi_data = SJ_PK * np.sin(2*np.pi*(F_SJ/RB)*n)

    phi_clk = np.zeros(N_UI)
    fi = np.zeros(N_UI)        # integral path; frequency correction
    df = np.zeros(N_UI)        # total normalized frequency correction
    bb = np.zeros(N_UI)

    last_bb = 0.0

    for k in range(N_UI - 1):
        phase_err_now = wrap_ui(phi_data[k] - phi_clk[k])

        # Transition-aware binary phase detector
        if TRANSITION[k]:
            if phase_err_now > 0:
                last_bb = +1.0
            elif phase_err_now < 0:
                last_bb = -1.0

        # Binary PD: no HOLD state
        bb[k] = last_bb

        # Explicit loop latency
        kd = k - td_ui
        cmd = bb[kd] if kd >= 0 else 0.0

        # PI loop filter
        fi[k + 1] = fi[k] + ki * cmd
        df[k + 1] = kp * cmd + fi[k + 1]

        # VCO: frequency -> phase
        phi_clk[k + 1] = phi_clk[k] + df[k + 1]

    phase_err = wrap_ui(phi_data - phi_clk)

    ss = slice(N_BURN, None)
    dphi_pp = np.ptp(phase_err[ss])

    return {
        "n": n,
        "phi_data": phi_data,
        "phi_clk": phi_clk,
        "phase_err": phase_err,
        "bb": bb,
        "fi": fi,
        "df": df,
        "dphi_pp": dphi_pp,
    }

reference

Hall, Stephen H., and Howard L. Heck. Advanced Signal Integrity for High-speed Digital Designs. Wiley : IEEE, 2009 [pdf]

Oh, Kyung, and Xing Yuan. High-Speed Signaling: Jitter Modeling, Analysis, and Budgeting. 1st edition. Prentice Hall, 2011. [pdf]

John M. Cioffi, [Chapter 3 - Equalization], [Chapter 6 - Fundamentals of Synchronization]

David Johns. ECE1392H - Integrated Circuits for Digital Communications - Fall 2001: [Equalization], [Timing Recovery]

B. Kim, "Tutorial: Basics of Equalization Techniques: Channels, Equalization, and Circuits," 2022 IEEE International Solid-State Circuits Conference (ISSCC), San Francisco, CA, USA, 2022

Masum Hossain, ISSCC2023 T11: "Digital Equalization and Timing Recovery Techniques for ADC-DSP-based Highspeed Links" [https://www.nishanchettri.com/isscc-slides/2023%20ISSCC/TUTORIALS/T11.pdf]

—, "LOW POWER DIGITAL EQUALIZATION FOR HIGH SPEED SERDES" [https://www.ieeetoronto.ca/wp-content/uploads/2020/06/SSCS_invited_talk.pdf]

Vivek Telang, 2012, Equalization for High-Speed Serdes: System-level Comparison of Analog and Digital Techniques [https://ewh.ieee.org/r5/denver/sscs/Presentations/2012_08_Telang.pdf]

Gain Kim, 2023. Equalization, Architecture, and Circuit Design for High-Speed Serial Link Receiver [pdf]

—, CICC2022 ES4: Equalization, Architecture, and Circuit Design for High-Speed Serial Link Receiver

S. Laxman, "Equalization algorithms in Millimeter wave communication systems," 2017 IEEE Custom Integrated Circuits Conference (CICC), Austin, TX, USA, 2017 [pdf]

A. Amirkhany, "Basics of Clock and Data Recovery Circuits: Exploring High-Speed Serial Links," in IEEE Solid-State Circuits Magazine, vol. 12, no. 1, pp. 25-38, Winter 2020 [https://sci-hub.jp/10.1109/MSSC.2019.2939342]

—, ISSCC2019 T6: "Basics of Clock and Data Recovery Circuits"

Fulvio Spagna, CICC2018 Clock and Data Recovery Systems [pdf]

Wei-Zen Chen, ISSCC2026. T9: Clocking and CDR Techniques for High-Performance Wireline Transceiver

B. Razavi, "The Design of a Clock and Data Recovery Circuit [The Analog Mind]," in IEEE Solid-State Circuits Magazine, vol. 18, no. 3, pp. 11-116, Summer 2026, doi: 10.1109/MSSC.2026.3706674.

A. A. Bazargani, H. Shakiba and D. A. Johns, "MMSE Equalizer Design Optimization for Wireline SerDes Applications," in IEEE Transactions on Circuits and Systems I: Regular Papers [https://www.eecg.utoronto.ca/~johns/nobots/papers/pdf/2024_bazaragani.pdf]

A. Sharif-Bakhtiar, A. Chan Carusone, "A Methodology for Accurate DFE Characterization," IEEE RFIC Symposium, Philadelphia, Pennsylvania, June 2018. [PDF] [Slides – PDF]

Tony Chan Carusone. High Speed Communications Part 11 – SerDes DSP Interactions [https://youtu.be/YIAwLskuVPc]

—, 2022 Optimization Tools for Future Wireline Transceivers [https://www.ieeetoronto.ca/wp-content/uploads/2022/12/UofT-Future-of-Wireline-Workshop-2022.pdf]

Alphawave IP CEO. How DSP is Killing the Analog in SerDes [https://youtu.be/OY2Dn4EDPiA]

S. Kiran, S. Cai, Y. Zhu, S. Hoyos and S. Palermo, "Digital Equalization With ADC-Based Receivers: Two Important Roles Played by Digital Signal Processingin Designing Analog-to-Digital-Converter-Based Wireline Communication Receivers," in IEEE Microwave Magazine, vol. 20, no. 5, pp. 62-79, May 2019 [https://sci-hub.se/10.1109/MMM.2019.2898025]

K. K. Parhi, "Design of multigigabit multiplexer-loop-based decision feedback equalizers," in IEEE Transactions on Very Large Scale Integration (VLSI) Systems, vol. 13, no. 4, pp. 489-493, April 2005 [http://sci-hub.se/10.1109/TVLSI.2004.842935]

T. Toifl et al., "A 3.5pJ/bit 8-tap-feed-forward 8-tap-decision feedback digital equalizer for 16Gb/s I/Os," ESSCIRC 2014 - 40th European Solid State Circuits Conference (ESSCIRC), Venice Lido, Italy, 2014 [https://sci-hub.se/10.1109/ESSCIRC.2014.6942120]

Daniel Friedman, 2018 Considerations and Implementations for High data Rate Serial Link Design [https://www.ieeetoronto.ca/wp-content/uploads/2020/06/DL-Toronto-Nov-2018.pdf]

Hongtao Zhang, DesignCon 2016. PAM4 Signaling for 56G Serial Link Applications − A Tutorial [https://www.xilinx.com/publications/events/designcon/2016/slides-pam4signalingfor56gserial-zhang-designcon.pdf]

Tony Chan Carusone Integrated Systems Laboratory, University of Toronto [https://isl.utoronto.ca/publications/]

Tony Chan Carusone 2022. Optimization Tools for Future Wireline Transceivers [https://www.ieeetoronto.ca/wp-content/uploads/2022/12/UofT-Future-of-Wireline-Workshop-2022.pdf]

Aleksey Tyshchenko, SeriaLink Systems Clinton Walker, Alphawave IP. DesignCon 2022. IBIS-AMI Modeling and Correlation Methodology for ADC-Based SerDes Beyond 100 Gb/s [https://static1.squarespace.com/static/5fb343ad64be791dab79a44f/t/63d807441bcd266de258b975/1675102025481/SLIDES_Track02_IBIS_AMI_Modeling_and_Correlation_Tyshchenko.pdf]

[https://ibis.org/summits/apr22/tyshchenko.pdf]

[https://www.mathworks.com/content/dam/mathworks/conference-or-academic-paper/ibis-ami-modeling-and-correlation.pdf]

Ali Sheikholeslami Electronics Group, University of Toronto [https://www.eecg.utoronto.ca/~ali/]

J. Liang, A. Sheikholeslami, H. Tamura, Y. Ogata and H. Yamaguchi, "Loop Gain Adaptation for Optimum Jitter Tolerance in Digital CDRs," in IEEE Journal of Solid-State Circuits, vol. 53, no. 9, pp. 2696-2708, Sept. 2018 [https://sci-hub.jp/10.1109/JSSC.2018.2839038]

BB PD

Youngdon Choi, Deog-Kyoon Jeong and W. Kim, "Jitter transfer analysis of tracked oversampling techniques for multigigabit clock and data recovery," in IEEE Transactions on Circuits and Systems II: Analog and Digital Signal Processing, vol. 50, no. 11, pp. 775-783, Nov. 2003 [https://sci-hub.st/10.1109/TCSII.2003.819070]

John T. Stonick, ISSCC 2011 TUTORIALS T5: DPLL-Based Clock and Data Recovery

Walker, Richard. (2003). Designing Bang-Bang PLLs for Clock and Data Recovery in Serial Data Transmission Systems. [pdf]

—, Clock and Data Recovery for Serial Data Communications, focusing on bang-bang CDR design methodology, ISSCC Short Course, February 2002. [slides]

It's ternary, because early, late and no transition

notice the transition density = 1 in digital PLL

Linearization

The effective PD gain is a function of the input jitter pdf, it enables one to anticipate the effects of input jitter on loop characteristics

BB Gain is the slope of average BB output \(\mu\), versus phase offset \(\phi\), i.e. \(\frac {\partial \mu}{\partial \phi}\),

BB only produces output for a transition and this de-rates the gain. Transition density = 0.5 for random data

\[ K_{BB} = \frac{1}{2}\frac {\partial \mu}{\partial \phi} \]

where \(\mu = (1)\times \mathrm{P}(\text{late}|\phi) + (-1)\times \mathrm{P}(\text{early}|\phi)\)

Both jitter and amplitude noise distribution are same, just scaled by slope

Self-Noise Term

One price we pay for BB PD versus linear PD is the self-noise term. For small phase errors BB output noise is the full magnitude of the sliced data

The PD output should be almost 0 for small phase errors. i.e. ideal PD output noise should be 0

\[ \sigma_{BB}^2 = 1^2 \cdot \mathrm{P}(\text{trans}) + 0^2\cdot (1-\mathrm{P}(\text{trans})) = 0.5 \]

Input referred jitter from BB PD is proportional to incoming jitter

gain simulation

L. Avallone, M. Mercandelli, A. Santiccioli, M. P. Kennedy, S. Levantino and C. Samori, "A Comprehensive Phase Noise Analysis of Bang-Bang Digital PLLs," in IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 68, no. 7, pp. 2775-2786, July 2021 [https://sci-hub.st/10.1109/TCSI.2021.3072344]

T. -K. Kuan and S. -I. Liu, "A Bang Bang Phase-Locked Loop Using Automatic Loop Gain Control and Loop Latency Reduction Techniques," in IEEE Journal of Solid-State Circuits, vol. 51, no. 4, pp. 821-831, April 2016 [https://sci-hub.st/10.1109/JSSC.2016.2519391]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

import matplotlib.pyplot as plt
import numpy as np

N = 2**10
sigma = 0.1
dt = np.random.normal(0, sigma, size=N)
et = np.sign(dt)

# Eq-(2)
coef_form = np.mean(np.abs(dt)) / np.mean(np.power(dt, 2))
print(f'coef_form: {coef_form}')

# Eq-(9)
coef_gauss = (2/np.pi)**0.5/sigma
print(f'coef_gauss: {coef_gauss}')

# polyfit
coef_fit = np.polyfit(dt, et, 1)
print(f'coef_fit: {coef_fit}')

x = np.linspace(-3.5, 3.5, 1000)
y = coef_fit[0]*x + coef_fit[1]

plt.figure(figsize=(12,6))
plt.plot(dt, et, 'o')
plt.plot(x, y, linewidth=2, linestyle='--')

# Calculate histogram counts and bin edges
counts, bin_edges = np.histogram(dt, bins=100)
# Find the maximum count
max_count = counts.max()
# Create weights to normalize the maximum height to 1
weights = np.ones_like(dt) / max_count
plt.hist(dt, bins=100, weights=weights)

plt.xlabel(r'$\Delta t$')
plt.grid(True)
plt.legend([r'$\Delta t \sim \varepsilon $', r'$x_{fit} \sim y_{fit}$', r'$ \text{hist}_{\Delta t}$'])
plt.show()


# coef_form: 7.82197790742685
# coef_gauss: 7.978845608028654
# coef_fit: [7.82251511 0.01010586]

Pavan, Shanthi, Richard Schreier, and Gabor Temes. (2016). Understanding Delta-Sigma Data Converters. 2nd ed. Wiley. - 2.2.1 Quantizer Modeling

\[ \frac{\mathrm{d}\sigma_e^2}{\mathrm{d}k} =0\space\space\Rightarrow\space\space k=\frac{\left\langle v,y\right\rangle}{\left\langle y,y \right\rangle} \]

DCO Quantization Noise

TODO 📅

TDC Quantization Noise

TODO 📅

Hunting Jitter

S. Jang, S. Kim, S. -H. Chu, G. -S. Jeong, Y. Kim and D. -K. Jeong, "An Optimum Loop Gain Tracking All-Digital PLL Using Autocorrelation of Bang–Bang Phase-Frequency Detection," in IEEE Transactions on Circuits and Systems II: Express Briefs, vol. 62, no. 9, pp. 836-840, Sept. 2015 [https://sci-hub.st/10.1109/TCSII.2015.2435691] [phd thesis]

Deog-Kyoon Jeong. Topics in IC (Wireline Transceiver Design). Lec 3 - All-Digital PLL [https://ocw.snu.ac.kr/sites/default/files/NOTE/Lec%203%20-%20ADPLL.pdf]

—. Topics in IC (Wireline Transceiver Design). Lec 6 - Clock and Data Recovery [https://ocw.snu.ac.kr/sites/default/files/NOTE/Lec%206%20-%20Clock%20and%20Data%20Recovery.pdf]

Lee Hae-Chang.: ‘An estimation approach to clock and data recovery’, PhD Thesis, Stanford University, November 2006 [https://www-vlsi.stanford.edu/people/alum/pdf/0611_HaechangLee_Phase_Estimation.pdf]

J. Kim, Design of CMOS Adaptive-Supply Serial Links, Ph.D. Thesis, Stanford University, December 2002. [https://vlsiweb.stanford.edu/people/alum/pdf/0212_Kim_______Design_Of_CMOS_AdaptiveSu.pdf)]

High-speed Serial Interface 2013. Lect. 16 – Clock and Data Recovery 3 [http://tera.yonsei.ac.kr/class/2013_1_2/lecture/Lect16_CDR-3.pdf]

CC Chen. Why Hunting Jitter Happens in CDR: The Role of Input Jitter and Latency? [https://youtu.be/hPDielPsFgY]

Hunting jitter is often referred to as dithering jitter, the periodic time error between data clock and input data, which exhibits a limit-cycle behavior

DT & CT Spectral Density

Sam Palermo, ECEN620: Network Theory Broadband Circuit Design Fall 2025 Lecture 9: Digital PLLs [https://people.engr.tamu.edu/spalermo/ecen620/lecture09_ee620_digital_PLLs.pdf]

Michael Perrott, August 14, 2008. Short Course On Phase-Locked Loops and Their Applications Day 4, AM Lecture Digital Frequency Synthesizers [https://www.cppsim.com/PLL_Lectures/day4_am.pdf]

—, "A modeling approach for Sigma Delta fractional-N frequency synthesizers allowing straightforward noise analysis," in IEEE Journal of Solid-State Circuits, vol. 37, no. 8, pp. 1028-1038, Aug. 2002 [https://www.cppsim.com/Publications/JNL/perrott_jssc02.pdf]

—. "Techniques for high data rate modulation and low power operation of fractional-N frequency synthesizers." 1997. [https://www.cppsim.com/Publications/Theses/perrott_phdthesis.pdf]

Hsu, Chun-Ming, Ph. D. Massachusetts Institute of Technology. "Techniques for high-performance digital frequency synthesis and phase control." 2008. [http://hdl.handle.net/1721.1/45870]

J. R. Barry, E. A. Lee, and D. G. Messerschmitt, Digital Communication, 3rd ed., Boston, MA: Kluwer Academic Publishers, 2003.

• \(\color{red}\frac{1}{T}\) of CT-DT originates from CT to sampled sequence to impulse train Fourier transform, which is explicit in frequency-domain model
• \(\color{red}\frac{1}{T}\) originates from DT spectrum to CT impulse spectrum is implicit in frequency-domain model
• \(\color{red}T\) of DT-CT originates from ZOH approximation following impulse train, which is explicit in frequency-domain model

Divider Sampling Operation & ITM

Impulse Train Modulator (ITM)

The double outline of the box in the figure is meant to serve as a reminder that a sampling operation is taking place

DT -> CT

\(\boxed{S_x(e^{j2\pi fT}) = S_d(f)\cdot \textcolor{blue}{\frac{1}{T}}=S_c(f)\cdot \textcolor{blue}{\frac{1}{T}}}\), For example, the quantization noise spectrum is given by \(S_x(e^{j2\pi fT}) = \frac{1}{12f_s}\cdot \frac{1}{T} = \frac{1}{12}\). To convert the sequence spectrum into a continuous impulse train spectrum, we multiply by \(\color{blue}\frac{1}{T^2}\) \[ S_y(f) = S_x(e^{j2\pi fT}) \cdot \textcolor{blue}{T\cdot \frac{1}{T^2}\cdot |H(f)|^2 } = S_x(e^{j2\pi fT}) \cdot \textcolor{blue}{\frac{1}{T}|H(f)|^2} \]

CT -> DT -> CT

Assume that the lower \(m\) bits of the digital filter output are discarded by truncation. The truncation error is therefore modeled as a uniformly distributed random variable,

\[ E_t \sim U[0,2^m\text{LSB}] \]

Because the effective output resolution is reduced, the new least significant bit becomes

\[ \text{LSB}_t = 2^m\text{LSB} \]

and the DAC quantization error is correspondingly modeled as

\[ Q_{DAC} \sim U[0,1] \space \text{in}\space \text{LSB}_t \]

Enhancing Resolution w/ DSM

J. Stonick. ISSCC 2011 tutorials, T5: "DPLL-Based Clock and Data Recovery"

Amir Amirkhany. ISSCC 2019 "Basics of Clock and Data Recovery Circuits"

\(M+1\) bits ensure on overflow or underflow in the signed adder

MRDT (Multi-rate Discrete-Time) Modeling

Y. Hu, T. Siriburanon and R. B. Staszewski, "Multirate Timestamp Modeling for Ultralow-Jitter Frequency Synthesis: A Tutorial," in IEEE Transactions on Circuits and Systems II: Express Briefs, vol. 69, no. 7, pp. 3030-3036, July 2022 [https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=9765581]

There are two key features associated with the behavior of DPLLs, namely, the multi-rate and discrete-time properties

reference & DCO model

timestamps with synchronous jitter for reference clock signal

periods with period jitter for free-running DCO

output jitter, PN & input jitter

Determine quantitatively the system jitters and PN from behavioral simulation of the MRDT DPLL

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

%% Wang, Xu and Michael Peter Kennedy. “Jitter and Spur Minimization in Fractional-N Digital Frequency Synthesizers - Modeling, Simulation, Analysis, and Design Methodologies.” *Analog Circuits and Signal Processing* (2026).

function [F_mean,t_n_rms,phi_n_rms,f_PSD,PSD_phi] = acquisition(t)
%% Ensure column vector
[~,ncol]=size(t);
if ncol ~= 1
t = t';
end
%% Periods
T          = diff(t);
%% Average period
T_mean     = mean(T);
F_mean     = 1/T_mean;
%% Period jitter
T_n        = T - T_mean;
%% Zero-mean accumulated time jitter
t_n        = cumsum(T_n);
t_n        = t_n - mean(t_n);
%% RMS time jitter
t_n_rms    = std(t_n);
%% Phase error
phi_n      = 2*pi/T_mean * t_n;
phi_n_rms  = std(phi_n);
%% PSD of phase error
npsd = 2^(nextpow2(length(phi_n)/8)-1);
Fpsd = 1/(npsd*T_mean);
f_PSD = 0:Fpsd:Fpsd*(floor(npsd/2)-1);
[PSD_phi,~] = pwelch(phi_n,window(@hann,npsd),npsd/2,npsd,1/T_mean, 'two-sided');
PSD_phi = PSD_phi(1:floor(npsd/2));

cross-domain scaling for PSD-consistent

N. Da Dalt, "Linearized Analysis of a Digital Bang-Bang PLL and Its Validity Limits Applied to Jitter Transfer and Jitter Generation," in IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 55, no. 11, pp. 3663-3675, Dec. 2008 [https://sci-hub.st/10.1109/TCSI.2008.925948]

—, “Theory and implementation of digital bang-bang frequency synthesizers for high speed serial data communications,” Ph.D. dissertation, RWTH Aachen University, Aachen, Germany, 2007. [https://publications.rwth-aachen.de/record/62439/files/DaDalt_Nicola.pdf]

H. Lu and P. P. Mercier, "Linear Periodically Time-Variant Digital PLL Phase Noise Modeling Using Conversion Matrices and Uncorrelated Upsampling," in IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 71, no. 12, pp. 6021-6033, Dec. 2024, doi: 10.1109/TCSI.2024.3415001

At low frequencies, the transfer function from \(t_{r}\) to \(T_{v}\) can be approximated as: \[ H_{t_r,T_v} \approx \frac{1-z^{-1}}{Nz^{-1}} \]

\[ \boxed{ \text{For deriving Eq. (10), the }\uparrow N\text{ block must be treated as unity-amplitude rate conversion} } \]

The statement

\[ \uparrow N_{\text{Da Dalt}}=N\,\uparrow N_{\text{standard}} \]

cannot be inserted as a scalar gain \(N\) when deriving Eq. (10). If we do that, the explicit \(1/N\) in Fig. 10 cancels:

\[ N\times \frac1N=1 \]

and Eq. (10) would become \(N\) times larger. Its DC phase gain would then be \(N^2\), instead of the correct \(N\).

Da Dalt only ever uses the second. Hence the printed \(1/N\).

Equation (10) is never used to trace a waveform, while exists to be squared and multiplied into a spectrum, via (9):

\[ S_{\phi_v}(f) = \left|H_{\phi_r,\phi_v}(f)\right|^2\cdot\left(S_{\phi_r}+S_{\phi_{\mathrm{bpd}}}(f)\right) + \left|H_{\phi_{\mathrm{dco}},\phi_v}(f)\right|^2 S_{\phi_{\mathrm{dco}}}(f) \]

and to be integrated for jitter variance in (13). Both uses require a PSD-consistent transfer function, and for a slow-in/fast-out path that is not the same as the transform ratio.

Write the cross-domain path as

\[ y[n] = \sum_k h[n-kN]\,x[k] \]

with \(x\) zero-mean, white, and of variance \(\sigma^2\) on the slow grid \(T_{r0}\), and with real-valued \(h\). Then

\[ E\{y^2[n]\} = \sigma^2\sum_k h^2[n-kN] \]

The variance is \(N\)-periodic, but it need not differ between output phases. The output is generally cyclostationary; even when its variance is constant, its autocorrelation can depend on \(n\bmod N\).

Averaging the variance over \(N\) fast samples: \[ \overline{\sigma_y^2} = \frac{\sigma^2}{N}\sum_j h^2[j] \]

Now demand the ordinary form \[ S_y = c\,|H|^2 S_x \]

With the paper's convention \(S_x(f) = T_{r0}\sigma^2\), integrating over the fast Nyquist band \(F = N/T_{r0}\) and applying Parseval \(\int_{-F/2}^{F/2}|H|^2\,df = F\sum_j h^2[j]\):

\[ \overline{\sigma_y^2} = c\,T_{r0}\sigma^2\cdot\frac{N}{T_{r0}}\sum_j h^2[j] \quad\Longrightarrow\quad \textcolor{red}{c = \frac{1}{N^2}} \]

Finally

\[ \boxed{\;S_y(f) = \frac{|H|^2}{N^2}\,S_x(f)\;} \]

So the transfer function you may legitimately plug into \(S_y = |H|^2 S_x\) is \(\textcolor{red}{H/N}\), not \(H\). That is equation (10)

For the hold stage alone, \(H_{\mathrm{ZOH}}(z)=(1-z^{-N})/(1-z^{-1})\), whose impulse response is rectangular. For the complete cross-domain path, \(h\) denotes the full effective impulse response.

The root cause in one line: \(S_x\) is normalized on \(T_{r0}\) while \(S_y\) is normalized on \(T_{v0} = T_{r0}/N\). The \(1/N\) reconciles the two normalizations.

Assume real-valued \(h\) and zero-mean white input:

\[ E\{x[k]\}=0,\qquad E\{x[k]x[\ell]\}=\sigma^2\delta_{k\ell}. \]

Here \(k\) indexes slow samples, while \(n\) indexes fast samples; one slow interval contains \(N\) fast samples.

Starting from

\[ y[n]=\sum_k h[n-kN]x[k], \]

we have \(E\{y[n]\}=0\), so its variance equals its second moment. Expanding the square:

\[ \begin{aligned} \operatorname{Var}(y[n]) &=E\!\left\{ \left(\sum_k h[n-kN]x[k]\right) \left(\sum_\ell h[n-\ell N]x[\ell]\right) \right\}\\ &=\sum_k\sum_\ell h[n-kN]h[n-\ell N]\, E\{x[k]x[\ell]\}\\ &=\boxed{\sigma^2\sum_k h^2[n-kN]}. \end{aligned} \]

Every term with \(k\ne\ell\) vanishes because distinct input samples are uncorrelated. Independence is unnecessary.

Write \(n=qN+r\), where \(r\in\{0,\ldots,N-1\}\). Then

\[ \operatorname{Var}(y[qN+r]) =\sigma^2\sum_k h^2[r-(k-q)N] =\sigma^2\sum_m h^2[r-mN]. \]

The result depends only on the phase \(r\), not the period number \(q\).

Thus, at each output phase, the variance uses only the filter coefficients whose indices have that particular remainder modulo \(N\).

Take the arithmetic average of these \(N\) phase variances:

\[ \overline{\sigma_y^2} =\frac1N\sum_{r=0}^{N-1}\operatorname{Var}(y[r]) =\frac{\sigma^2}{N} \sum_{r=0}^{N-1}\sum_k h^2[r-kN]. \]

Every integer \(j\) has exactly one representation

\[ j=r-kN,\qquad 0\le r<N. \]

Consequently, the double sum includes every \(h^2[j]\) exactly once:

\[ \boxed{\overline{\sigma_y^2} =\frac{\sigma^2}{N}\sum_j h^2[j]}. \]

This is the average of the variances, not the variance of an averaged output signal.

\[ \boxed{\text{Physically: zero stuff }T_v\rightarrow\text{ZOH}\rightarrow\text{accumulate; no }1/N} \]

Because \(f_{\mathrm{mod}} = 1\,\mathrm{kHz}\) is near DC relative to the loop bandwidth, you should get approximately \[ j_v(t) \approx j_r(t) \]

and therefore

\[ \frac{A_{\mathrm{out}}}{A_{\mathrm{in}}} \approx 1 \]

or

\[ 20\log_{10}\left|\frac{J_v}{J_r}\right| \approx 0\,\mathrm{dB}. \]

Here \(J_v\) and \(J_r\) are the Fourier components of the output and reference timing jitter at \(f_{\mathrm{mod}}\).

At the same time,

\[ \Delta t = j_r - j_v \approx 0. \]

[Github Gist — dpll_in_out.py]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

# ============================================================
# DCO OUTPUT TIMESTAMPS
#
# The DLF output is held for N DCO cycles.
# ============================================================

# DCO period error due to loop control
dTctrl = (KT * u)

# Total DCO period error for every fast-clock cycle
dT = np.repeat(dTctrl, N)

# Limit to exactly Nsim*N DCO cycles
dT = dT[:Nsim * N]

# Accumulate timing error
jv_fast = np.concatenate(
    ([0.0], np.cumsum(dT))
)

# ------------------------------------------------------------
# Actual DCO timestamps
# ------------------------------------------------------------
hdco = np.arange(Nsim * N + 1)

tv = hdco * Tv0 + jv_fast

# ============================================================
# OUTPUT JITTER FROM tv
#
# jv[h] = tv[h] - h*Tv0
# ============================================================
jv_from_tv = tv - hdco * Tv0

So, specifically for Python model:

\[ \boxed{ \texttt{np.repeat(x,N)} = \uparrow N+ \frac{1-z^{-N}}{1-z^{-1}} } \]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

============================================================
DPLL LOW-FREQUENCY INPUT JITTER TEST
============================================================
Reference frequency        : 100.000 MHz
DCO frequency              : 2.400 GHz
Frequency division         : 24

Input jitter frequency     : 1000.000 Hz
Input jitter amplitude     : 1.004880 ps

Output jitter amplitude    : 1.004533 ps
Residual error amplitude   : 0.015916 ps

|Jv/Jr|                    : 0.99965447
Jv/Jr gain                 : -0.003002 dB

Input RMS jitter           : 708.703382 fs
Output RMS jitter          : 729.245449 fs
Residual RMS error         : 172.582105 fs

|Error/Jr|                 : 1.583866e-02
Error transfer gain        : -36.006 dB
============================================================

!!DPLL time-domain model

L. Avallone, M. Mercandelli, A. Santiccioli, M. P. Kennedy, S. Levantino and C. Samori, "A Comprehensive Phase Noise Analysis of Bang-Bang Digital PLLs," in IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 68, no. 7, pp. 2775-2786, July 2021 [https://sci-hub.st/10.1109/TCSI.2021.3072344]

—, “Contributions to the Theory and Development of Low-Jitter Bang-Bang Integrated Frequency Synthesizers.” University College Dublin. School of Electrical and Electronic Engineering, 2022. [http://hdl.handle.net/10197/13372]

The paper defines the BPD input as

\[ \boxed{\Delta t[k]=t_r[k]-t_d[k]} \]

where \(t_r[k]\) is the reference-edge timestamp and \(t_d[k]\) is the divider-output timestamp.

The signal chain is essentially \[ \epsilon[k] \rightarrow \underbrace{\beta\epsilon[k]+\psi[k]}_{u[k]} \rightarrow \underbrace{K_T u[k]}_{\text{DCO period change}} \rightarrow T_v. \]

A useful distinction is that \(K_T\) is not the usual DCO frequency gain \(K_{\mathrm{DCO}}\) in Hz/code. This paper models the DCO in the period domain, so its gain is period/code. Around the nominal operating point, \(f=\frac{1}{T}\)

hence for a small period change,

\[ \Delta f \approx -\frac{\Delta T}{T_0^2}. \]

Therefore the corresponding frequency gain would be approximately

\[ \boxed{ K_{\mathrm{DCO}} \approx -\frac{K_T}{T_0^2} = -K_T f_v^2 } \]

in Hz/code. The minus sign means increasing the period lowers the frequency.

with \(t_A = t_B - d_t[k]=j_r[k]-d_t[k]\) \[ t_C = t_A + NT_{v0}+NK_T u[k]+W_v[k] = j_r[k]-d_t[k] + NT_{v0}+NK_T u[k]+W_v[k] \] with \(t_D = N T_{v0} + j_r[k+1]\) \[ \textcolor{red}{d_t[k+1]} = t_D - t_C = \textcolor{red}{\boxed{d_t[k] + (j_r[k+1] - j_r[k]) - NK_T u[k]-W_v[k]}} \]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166

%% DPLL parameters from the paper
fr    = 100e6;
N     = 24;
fv    = N*fr;
Tv0   = 1/fv;             % nominal free-running DCO period [s]

beta  = 70;
alpha = beta/2^8;

KT = 0.145e-15;          % [s/bit]

df = 1e6;                % DCO PN specified at 1 MHz

%% Table-I cases
%           Ref PN       DCO PN @ 1 MHz
cases = [   -155,        -Inf;
            -Inf,        -108;
            -155,        -108;
            -150,        -108;
            -155,        -103];

Nsim  = 1e6;
Nburn = 2e4;


rng(1);

sigma_sim = zeros(size(cases,1),1);

for icase = 1:size(cases,1)

    Lref = cases(icase,1);
    Ldco = cases(icase,2);

    %% ------------------------------------------------------------
    % Reference white phase noise -> absolute timestamp jitter
    %% ------------------------------------------------------------
    if isinf(Lref)
        sigma_ref = 0;
    else
        Sphi_ref = 10^(Lref/10);

        sigma_ref = sqrt( ...
            Sphi_ref/(4*pi^2*fr) );
    end

    %% ------------------------------------------------------------
    % DCO 1/f^2 phase noise -> DCO cycle jitter
    % Eq. (15)
    %% ------------------------------------------------------------
    if isinf(Ldco)
        sigma_Tv = 0;
    else
        Sphi_dco = 10^(Ldco/10);

        sigma_Tv = sqrt( ...
            Sphi_dco * df^2 / fv^3 );
    end

    fprintf('\nCase %d\n',icase);
    fprintf('sigma_ref = %.3f fs\n',sigma_ref/1e-15);
    fprintf('sigma_Tv  = %.3f fs\n',sigma_Tv/1e-15);

    %% Reference absolute jitter
    jr = sigma_ref * randn(Nsim+1,1);

    %% DCO cycle-period errors at the fast clock rate
    %
    % Tv(:,k) holds the N zero-mean period-noise samples in reference
    % interval k.  Their sum is the Wv[k] term in the reference-rate
    % recursion.  Keeping the individual samples makes it possible to
    % construct the path-consistent DCO timestamps tv[h] from Eq. (4).
    %
    if sigma_Tv == 0
        Tv = [];
        Wv = zeros(Nsim,1);
    else
        Tv = sigma_Tv*randn(N,Nsim);
        Wv = sum(Tv,1).';
    end

    %% ------------------------------------------------------------
    % Bang-bang DPLL
    %% ------------------------------------------------------------

    dt  = zeros(Nsim+1,1);
    u   = zeros(Nsim,1);
    psi = 0;

    % t_r[0] - t_d[0]
    dt(1) = jr(1);  % or w/ Zero-initialization

    for k = 1:Nsim

        %% Binary phase detector
        if dt(k) >= 0
            epsilon = +1;
        else
            epsilon = -1;
        end

        %% Integral path
        psi = psi + alpha*epsilon;

        %% PI loop-filter output
        u(k) = beta*epsilon + psi;

        %% Time-error recursion
        dt(k+1) = dt(k) ...
                + (jr(k+1)-jr(k)) ...	% absolute jitter to period jitter
                - N*KT*u(k) ...
                - Wv(k);

    end

    %% ------------------------------------------------------------
    % DCO-output timestamps, Eq. (4)
    %% ------------------------------------------------------------
    % The DLF output is zero-order held for N DCO cycles.  Add that
    % control contribution to the SAME fast-rate noise samples whose
    % block sums drove the recursion above.  Accumulate the period error
    % first so nominal-period roundoff cannot mask femtosecond jitter.
    dTctrl = (KT*u).';
    if isempty(Tv)
        Tv = repmat(dTctrl,N,1);
    else
        for m = 1:N
            Tv(m,:) = Tv(m,:) + dTctrl;
        end
    end
    % Reuse Tv for the accumulated timing error to limit the peak to
    % roughly two full DCO-rate traces during the leading-zero prepend.
    clear Wv u dTctrl
    Tv = cumsum(Tv(:));             % t_v[h] - h*Tv0, h = 1,...,N*Nsim
    tv = [0; Tv];                   % include t_v[0]
    clear Tv

    % Form literal edge timestamps without cumulatively adding Tv0.
    % Chunking avoids another full-size temporary vector.
    Ndco = N*Nsim;
    edge_chunk = 1e6;
    for h0 = 0:edge_chunk:Ndco
        h1 = min(h0+edge_chunk-1,Ndco);
        h = (h0:h1).';
        idx = h+1;
        tv(idx) = tv(idx) + h*Tv0;
    end


    %% Remove startup transient
    dt_ss = dt(Nburn+1:end);

    sigma_sim(icase) = std(dt_ss);

end

%% Results
fprintf('\n------------------------------------------\n');
fprintf('Case       sigma_Dt [fs]\n');
fprintf('------------------------------------------\n');

for k = 1:length(sigma_sim)
    fprintf('(%c)        %8.2f\n', ...
        'a'+k-1, sigma_sim(k)/1e-15);
end

There are two different random sequences: \[ \boxed{j_r[k] = \text{absolute reference edge jitter}} \]

versus

\[ \boxed{\delta T_r[k]=j_r[k+1]-j_r[k] =\text{reference period jitter}}. \]

Our simulation generates

1

jr = sigma_ref * randn(...);

because the paper assumes white absolute reference jitter. The paper explicitly calls \(\sigma_{t_r}^2\) the absolute jitter variance of the reference.

Then the model naturally converts that absolute jitter into reference-period variation using

1

jr(k+1) - jr(k)

There is also a subtle but important consequence. If \(j_r[k]\sim\mathcal N(0,\sigma_{t_r}^2)\) is i.i.d., then \(\operatorname{Var}\{j_r[k+1]-j_r[k]\} = 2\sigma_{t_r}^2.\)

But consecutive period errors are correlated:

\[ \operatorname{Cov} \left( j_r[k+1]-j_r[k], j_r[k+2]-j_r[k+1] \right) = -\sigma_{t_r}^2. \]

So you should not replace the code with independent samples such as

1

djr = sqrt(2)*sigma_ref*randn(...);

because that gets the variance right but loses the required correlation.

In short:

\[ \boxed{ \texttt{jr[k]}=\text{edge jitter} \quad\Rightarrow\quad \texttt{jr[k+1]-jr[k]}=\text{period jitter} } \]

and the DPLL recursion evolves from one edge interval to the next, which is why the difference appears.

\(\boxed{\sigma_{\Delta t}}\) at the BPD input, not directly the RMS DCO-output jitter. To obtain actual DCO-output jitter, we should also simulate/store \(t_v[h]\), rather than only the reference-rate recursion for \(\Delta t[k]\).

causality z^-1 in linear model

\(\textcolor{red}{z^{-1}}\) makes the loop causal, and it represents the one-reference-cycle latency that is physically unavoidable in a digital PLL

With the indexing used in this slide,

\[ \boxed{\Delta t[k]=t_r[k]-t_v[k-1]} \]

so in the \(z\)-domain,

\[ \Delta T(z)=T_r(z)-z^{-1}T_v(z) \]

That is exactly the \(z^{-1}\) shown in the feedback path.

The timing sequence is essentially

\[ t_v[k-1] \;\longrightarrow\; \Delta t[k] \;\longrightarrow\; \text{TDC/filter} \;\longrightarrow\; u[k] \;\longrightarrow\; t_v[k] \]

Meanwhile the DCO timing recursion is

\[ t_v[k]=t_v[k-1]+N K_T u[k] \]

which gives

\[ \frac{T_v(z)}{U(z)} = \frac{N K_T}{1-z^{-1}} \]

So the two appearances of \(z^{-1}\) have related but different meanings:

• Inside \(\displaystyle \frac{NK_T}{1-z^{-1}}\): the DCO accumulates period/time increments.
• In the feedback path: the TDC at iteration \(k\) sees the previous DCO/divider timing state \(t_v[k-1]\).

A useful way to read the lower figure is therefore:

\[ \boxed{ t_r[k] - \underbrace{t_v[k-1]}_{\text{available feedback edge}} \rightarrow \text{TDC} \rightarrow H(z) \rightarrow \text{DCO} \rightarrow t_v[k] } \]

Model at Reference Rate

S. Levantino, "Digital phase-locked loops," 2018 IEEE Custom Integrated Circuits Conference (CICC), San Diego, CA, USA, 2018

—, "Advanced digital phase-locked loops," Proceedings of the IEEE 2013 Custom Integrated Circuits Conference, San Jose, CA, USA, 2013, pp. 1-95 [https://sci-hub.jp/10.1109/CICC.2013.6658505]

\[ \boxed{\frac{T_V}{W}(z) = \frac{NK_T}{z-1} = \frac{NK_T}{\textcolor{red}{1-z^{-1}}}\cdot\textcolor{red}{z^{-1}}} \]

\[\begin{align} y[n] &= y[n-1] + x[n-1] \quad &\Rightarrow\quad\quad \bbox[yellow]{\frac{Y}{X}(z) = \frac{z^{-1}}{1-z^{-1}}} \\ y[n] &= y[n-1] + x[n] \quad &\Rightarrow\quad\quad \frac{Y}{X}(z) = \frac{1}{1-z^{-1}} \end{align}\]

N. Da Dalt, "Linearized Analysis of a Digital Bang-Bang PLL and Its Validity Limits Applied to Jitter Transfer and Jitter Generation," in IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 55, no. 11, pp. 3663-3675, Dec. 2008 [https://sci-hub.st/10.1109/TCSI.2008.925948]

Model at DCO Rate

L. Avallone, M. Mercandelli, A. Santiccioli, M. P. Kennedy, S. Levantino and C. Samori, "A Comprehensive Phase Noise Analysis of Bang-Bang Digital PLLs," in IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 68, no. 7, pp. 2775-2786, July 2021 [https://sci-hub.st/10.1109/TCSI.2021.3072344]

t_r[h] in Fig. 5 is up-sampled and zero-padded from t_r[k] in Fig. 2. — zero-padded transform of a signal in the slow domain \[ \boxed{\Delta t[h] = \left\{ \begin{array}{cl} t_r[h]-t_d[h] & \text{when } h=0,\pm N,\pm 2N,... \\ 0 & \text{otherwise} \end{array} \right.} \]

Taking the \(\mathcal{Z}\)-transform of the discrete-time accumulator

\[ t_v^\prime[k] = t_v^\prime[k-1] + w[k-1] \]

results in the system function

\[ \boxed{\frac{T_v^\prime}{W}(z) = \frac{z^{-1}}{1-z^{-1}}} \]

reference

Wang, Xu and Michael Peter Kennedy. “Jitter and Spur Minimization in Fractional-N Digital Frequency Synthesizers - Modeling, Simulation, Analysis, and Design Methodologies.” Analog Circuits and Signal Processing (2026).

Brandonisio, F., & Kennedy, M. P. (2014). Noise-Shaping All-Digital Phase-Locked Loops: Modeling, Simulation, Analysis and Design. Springer.

Staszewski, Robert Bogdan and Poras T. Balsara. “All-digital frequency synthesizer in deep-submicron CMOS.” (2006).

Topics in IC (Wireline Transceiver Design) [https://ocw.snu.ac.kr/sites/default/files/NOTE/Lec%203%20-%20ADPLL.pdf]

Michael H. Perrott, ISSCC 2008 Tutorial on Digital Phase-Locked Loops

—, CICC 2009 Tutorial on Digital Phase-Locked Loops [https://www.cppsim.com/PLL_Lectures/digital_pll_cicc_tutorial_perrott.pdf]

Robert Bogdan Staszewski, CICC 2020: Beyond All-Digital PLL for RF and Millimeter-Wave Frequency Synthesis [link]

Akihide Sai, ISSCC 2023 T5: All-digital PLLs From Fundamental Concepts to Future Trends

Mike Shuo-Wei Chen, CICC 2020 ES2-3: Low-Spur PLL Architectures and Techniques [https://youtu.be/sgPDchYhN-4]

Saurabh Saxena, IIT Madras. Phase-Locked Loops: Noise Analysis in Digital PLL [https://youtu.be/mddtxcqfiKU]

Neil Robertson. Digital PLL's -- Part 1 [https://www.dsprelated.com/showarticle/967.php]

—. Digital PLL's -- Part 2 [https://www.dsprelated.com/showarticle/973.php]

—. Digital PLL's -- Part 3 [https://www.dsprelated.com/showarticle/1177.php]

Daniel Boschen. GRCon24 - Quick Start on Control Loops with Python Workshop [video, slides]

M. Zanuso, D. Tasca, S. Levantino, A. Donadel, C. Samori and A. L. Lacaita, "Noise Analysis and Minimization in Bang-Bang Digital PLLs," in IEEE Transactions on Circuits and Systems II: Express Briefs, vol. 56, no. 11, pp. 835-839, Nov. 2009 [https://sci-hub.st/10.1109/TCSII.2009.2032470]

N. Da Dalt, "Markov Chains-Based Derivation of the Phase Detector Gain in Bang-Bang PLLs," in IEEE Transactions on Circuits and Systems II: Express Briefs, vol. 53, no. 11, pp. 1195-1199, Nov. 2006 [https://sci-hub.st/10.1109/TCSII.2006.883197]

—, "A design-oriented study of the nonlinear dynamics of digital bang-bang PLLs," in IEEE Transactions on Circuits and Systems I: Regular Papers, vol. 52, no. 1, pp. 21-31, Jan. 2005 [https://sci-hub.se/10.1109/TCSI.2004.840089]

—, "Theory and Implementation of Digital Bang-Bang Frequency Synthesizers for High Speed Serial Data Communications", PhD Dissertation, RWTH Aachen University, Aachen, North Rhine-Westphalia, Germany, 2007 [pdf]

TDR Impedance

Amit Bahl, How TDR Impedance Measurements Work [https://www.protoexpress.com/blog/tdr-impedance-measurements/]

Minh Quach. Signal Integrity Consideration and Analysis 4/30/2004 Frequency & Time Domain Measurements/Analysis [https://ewh.ieee.org/r5/denver/sscs/Presentations/2004_04_Quach.pdf]

江上渔樵, 在ADS中查看TDR的3种方法 [https://zhuanlan.zhihu.com/p/420350734]

There are two primary approaches for obtaining a TDR response:

1. Measurement in the time domain using a sampling scope
2. Simulation from measured frequency-domain S-parameters (Not every SI lab is equipped with a TDR instrument, but most have a VNA)

Why TDR Impedance?

Measured time domain TDR

Abhargava, TDR Analysis using Agilent ADS [https://abhargava.wordpress.com/wp-content/uploads/2014/01/performing-tdr-analysis-using-agilent-ads.pdf]

Mike Steinberger, TDR: Reading the Tea Leaves [https://siguys.com/wp-content/uploads/2016/01/TDR_TeaLeaves.pdf]

\[ \color{red}Z_T(t) = Z_0\cdot \frac{1+\Gamma(t)}{1-\Gamma(t)} \]

S-parameter-converted TDR

Vladimir Dmitriev-Zdorov, Mentor Graphics, DesignCon 2014, Computation of Time Domain Impedance Profile from S-Parameters: Challenges and Methods [link]

Samtec, High Speed Characterization Report PCIEC-064-1000-EC-EM-P-85 [https://suddendocs.samtec.com/testreports/hsc-report_pciec-85_web.pdf]

w/ IFFT

比尔盖子, Frequency domain S11 conversion to time domain TDR [https://electronics.stackexchange.com/a/626063/233816]

HFSS™ 3D Layout Window Functions and Time Domain Plotting [https://ansyshelp.ansys.com/public/Views/Secured/Electronics/v252/en/Subsystems/HFSS3DLayout/Content/ReportsandPostProc/WindowFunctionsandTimeDomainPlotting.htm]

Time Domain Measurements using Vector Network Analyzer ZVR [https://scdn.rohde-schwarz.com/ur/pws/dl_downloads/dl_application/application_notes/1ez44/1ez44_0e.pdf]

scikit-rf plot_z_time_step

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

S₁₁(f)
   │  (extrapolate_to_dc → uniform grid starting at 0 Hz)
   ▼
W(f) · S₁₁(f)                 ← windowed(): half-window, 1 at DC, 0 at f_max,
                                normalize=False (preserves S(0))
   │  (np.fft.irfft + fftshift)
   ▼
h(t)  – impulse response       ← impulse_response()
   │  (cumulative_trapezoid)
   ▼
Γ_step(t)  ∈ [−1, 1]           ← step_response()
   │  Z(t) = z0 · (1+Γ)/(1−Γ)   (clamp Γ=1)
   ▼
Z(t)  – plotted vs t (in ns)   ← plot_attribute() with attribute='z',
                                conversion='time_step'

Window function with \(w[0] = 1\) and \(w[f_{max}]\)=0 ensure \(\Gamma(+\infty)\) and \(\Gamma(0)\) are correct

Jim Nadolny, Samtec. Technical Note Transformation of Samtec Connector Test Data For 85 ohm Differential Impedance Applications, [https://suddendocs.samtec.com/notesandwhitepapers/technical-note_85ohm-reference-z-xform_web.pdf]

[ADS: 1-10] TDR Impedance (Part 2) TDRインピーダンス解析 [https://youtu.be/ACINktqpM50]

ADS tdr_sp_imped

Peter Goossens, Transformation of time domain TDR to its frequency domain S11 (Return Loss) using FFT [https://www.gquipment.com/blog/transformation-of-time-domain-tdr-to-its-frequency-domain-s11-return-loss-using-fft]

\[ \color{green} \Gamma(t) \to S_{ii} \]

1
2
3
4
5
6
7
8
9
10
11
12
13
14

# Pseudo code, laying out the essential steps only
# Trace data is in x_values
sampling_interval = np.mean(np.diff(x_values))
#  First derivative
diff_gamma_values = np.gradient(gamma_values)
# FFT
fourier_data = fft(diff_gama_values)
# Get the frequencies corresponding to the FFT result
frequencies = np.fft.fftfreq(len(diff_gama_values), d=sampling_interval)
# Calculate the magnitude of the complex Fourier transform data
magnitude = np.abs(fourier_data)
# Return loss
magnitude = 20 * np.log10(magnitude)
# Plot (frequency, magnitude)

w/ RFE

rational fraction expansion (RFE)

TODO 📅

S-Parameter Properties Shape Simulated TDR

Tim Wang Lee, Ph.D., Keysight Technologies, When Frequency Shapes Time: How S-Parameter Properties Shape Simulated TDR Behavior [https://www.signalintegrityjournal.com/articles/4287-when-frequency-shapes-time-how-s-parameter-properties-shape-simulated-tdr-behavior]

In principle, both methods should yield identical impedance profiles. In practice, differences in S-parameter quality, bandwidth, and simulation setup can lead to noticeable discrepancies. The poor quality of the S-parameter data and improper simulation parameters caused the two traces to differ

Non-causal data — Early arrival, pre-ringing

\[ \boxed{S(t) = \int_{-\infty}^{t} h(\tau) d\tau, \text{ with } S(\infty) = H(0)} \] pure phase manipulation (linear phase or nonlinear phase) preserves the total area, and the TDR settles to the same final impedance on the flat sections

Time-Domain Transmission (TDT)

TODO 📅

Reading S-parameters

teledynelecroy. Reading S-parameters [https://blog.teledynelecroy.com/2020/05/]

keysight. How to Interpret Ripple in an S Parameters Measurement [https://docs.keysight.com/kkbopen/how-to-interpret-ripple-in-an-s-parameters-measurement-849642201.html]

You Measured What? Four Must-Know Checks Before Trusting Your Trace S-Parameters [https://www.signalintegrityjournal.com/articles/4083-you-measured-what-four-must-know-checks-before-trusting-your-trace-s-parameters]

TODO 📅

Ripple in an S Parameters

high-reflection channel

Christiaan Bil, Intel. Challenges in Setting Appropriate NUI Value in Measured Step Response-based 32 GT/s Rx Calibration

Using S Parameters to Estimate Q

Jeff Walling. ECE 5984 Using S Parameters to Estimate Q [https://youtu.be/PXgM6pGIRvk]

TODO 📅

oscilloscope bandwidth

Realtime oscilloscope bandwidth considerations for 25 Gbps PAM4 patterns, [https://www.ieee802.org/3/cy/public/adhoc/chang_3cy_01_12_20_22.pdf]

Bessel Thomson Filter Bandwidth

reference

Bogatin, Eric. 2020. Bogatin’s Practical Guide to Transmission Line Design and Characterization for Signal Integrity Applications / .Eric Bogatin. Artech House.

keysight, Signal Integrity Characterization Techniques [pdf]

Tim Wang-Lee, DesignCon 2026 KEF: Mastering TDR and De-embedding Through Simulation and Measurement [link]

Csaba SOOS, Signal and Power Integrity Design Practices [https://indico.cern.ch/event/358837/attachments/714663/1930957/Signal_and_Power_Integrity_Practices.pdf]

Sam Palermo, ECEN689: Special Topics in High-Speed Links Circuits and Systems Spring 2012 Lecture 3: Time-Domain Reflectometry & S-Parameter Channel Models [https://people.engr.tamu.edu/spalermo/ecen689/lecture3_ee689_tdr_spar.pdf]

Martin Stumpf, Preparing for PCIe® Electrical Measurements Beyond 64 GT/s

Electric field coupling (also called capacitive coupling) occurs when energy is coupled from one circuit to another through an electric field

Magnetic field coupling (also called inductive coupling) occurs when energy is coupled from one circuit to another through a magnetic field

For instance

• magnetic coupling between multiple inductors
• capacitive coupling between multiple transmission lines

Transmission Line

param. extraction from ABCD matrix

Chapter 4.5. High Frequency Passive Devices [https://www.cambridge.org/il/files/7713/6698/2369/HFIC_chapter_4_passives.pdf]

for lossless T-line, \(\gamma = j\beta\)

Capacitive Coupling

Faraday cage

Magnetic Coupling

Darabi H. Radio Frequency Integrated Circuits and Systems. 2nd ed. Cambridge University Press; 2020.

NEXT & FEXT

Backward (near-end) crosstalk & Forward (far-end) crosstalk

Mohammad Abu Khater, ISCAS2019 tutorial: High-Performance Printed Circuit Boards (PCBs)

Consider a small section at distance (x) from the input:

1. The aggressor’s edge reaches it at time \(x/v\), generating a small noise pulse.
2. That pulse travels forward on the victim through the remaining distance \(\ell-x\)

Its arrival time at the far end is therefore

\[ t_{\text{arrival}} =\underbrace{\frac{x}{v}}_{\text{aggressor reaches section}} +\underbrace{\frac{\ell-x}{v}}_{\text{noise reaches far end}} =\frac{\ell}{v} \]

Noise generated earlier travels farther; noise generated later travels less. The pulses overlap, so adding more coupled sections increases their summed amplitude

Each short section contributes in proportion to its length and the edge slope: \[ dV_F\propto dx\,\frac{\Delta V}{t_r} \quad\Longrightarrow\quad \boxed{V_{F,\text{peak}}\propto \ell\,\frac{\Delta V}{t_r}} \]

Here, \(\ell\) means the length over which the traces run alongside each other.

For comparison, backward noise arrives at \(x/v+x/v=2x/v\), so contributions from different locations spread out in time. That explains why extending a sufficiently long coupled line mainly increases the backward pulse’s duration

relative dielectric constant vs permittivity

The relative dielectric constant characterizes some of the electrical properties of an insulator

Return Path

ISSCC2002. Special Topic Evening Discussion Sessions SE1: Inductance: Implications and Solutions for High-Speed Digital Circuits [vSE1_Blaauw], [vSE1_Gauthier], [vSE1_Morton, [vSE1_Restle]]

Current return paths are frequency dependent \(Z = R +j\omega L\)

• Low frequency
  • \(R\) dominates - current use as many returns as possible to have parallel resistances
• High frequency
  • \(j\omega L\)​ dominates - current use the closest possible return path to form the smallest possible loop inductance
• Very high frequency
  • The current would be confined to the nearest possible return only at ultra-high frequencies (skin effect)

skin effect & Dielectric loss

EMX simulation

setup:

frequency sweep:

Cadence October 2020, Analysis of a Figure-Eight Inductor with EMX RAK

Tline Approximation

[https://web.stanford.edu/class/archive/ee/ee371/ee371.1066/handouts/markChapt.pdf]

N-section LC Model

Eric Bogatin. Pop Quiz: When is an Interconnect Not a Transmission Line? [https://www.signalintegrityjournal.com/blogs/4-eric-bogatin-signal-integrity-journal-technical-editor/post/265-pop-quiz-when-is-an-interconnect-not-a-transmission-line]

RLGC by Open/Short Circuit

RLGC can be extracted from measurements of a transmission line's input impedance under open-circuit and short-circuit terminations at a specific frequency

Dr. Muehlhaus Consulting & Software GmbH, lumpedmodel [https://github.com/VolkerMuehlhaus/lumpedmodel]

Transmission line from S2P data into RLGC lumped model

\[ \boxed{R= \text{Re}(\gamma Z_c)} \qquad \boxed{L= \frac{\text{Im}(\gamma Z_c)}{\omega}} \qquad \boxed{G= \text{Re}\left(\frac{\gamma}{Z_c}\right)} \qquad \boxed{C= \frac{\text{Im}\left(\frac{\gamma}{Z_c}\right)}{\omega}} \]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39

# https://github.com/VolkerMuehlhaus/lumpedmodel/blob/main/rlgc_from_s2p/rlgc_from_s2p.py


# input data, must be 2-port S2P data
sub = rf.Network(args.s2p, z0=args.z0_ohm)

# physical length must be supplied by user
length = args.l_um*1e-6

# target frequency for pi model extraction
f_target = args.f_ghz*1e9

assert f_target < freq.stop

# get index for exctraction
f = freq.f
ftarget_index = rf.find_nearest_index(freq.f, f_target)
omega = 2*np.pi*f[ftarget_index]

z11=sub.z[0::,0,0]          # z11, the open impedance
y11=sub.y[0::,0,0]          # 1/y11, the short impedance
Zline = np.sqrt(z11/y11)    # characteristic impedance of the line

gamma0 = 1/length * np.arctanh(1/(Zline*y11))   # propagation constant

# electrical phase  = β · length        = gamma0.imag * length   ✓
# attenuation       = α                 = gamma0.real            (no wrap)
# period=np.pi/2 instead of the default 2π reflects the branch period of arctanh
gamma_wideband =  gamma0.real + 1j*np.unwrap(gamma0.imag*length, period=np.pi/2)/length


gamma_ftarget = gamma_wideband[ftarget_index]
Zline_ftarget = Zline[ftarget_index]


R = (gamma_ftarget*Zline_ftarget).real
L = (gamma_ftarget*Zline_ftarget).imag / omega
G = (gamma_ftarget/Zline_ftarget).real
C = (gamma_ftarget/Zline_ftarget).imag / omega

Transmission Line [pdf]

Decoupling Capacitor

Grounding

Chapter 11 Layout and grounding [http://ieb-srv1.upc.es/gieb/tecniques/doc/EMC/pdfs/ScienceDirect_articles_27Jul2018_12-16-10.699/Chapter-11---Layout-and-grounding_2007_EMC-for-Product-Designers.pdf]

TODO

90^o Turns

Mohammad Abu Khater, ISCAS2019 tutorial: High-Performance Printed Circuit Boards (PCBs)

Hyperbolic Functions

reference

信号完整性揭秘:于博士SI设计手记

Bogatin, E. (2018). Signal and power integrity, simplified. Prentice Hall. [pdf]

High-speed Serial Interface Lect. 9 – Noise [http://tera.yonsei.ac.kr/class/2017_2_2/lecture/Lect%209%20Noise.pdf]

Yuriy Shlepnev. How Interconnects Work: Characteristic Impedance and Reflections [https://www.linkedin.com/pulse/how-interconnects-work-characteristic-impedance-yuriy-shlepnev/]

—. How Interconnects Work: Bandwidth for Modeling and Measurements [https://www.linkedin.com/pulse/how-interconnects-work-bandwidth-modeling-yuriy-shlepnev/?trackingId=874kpm3XuNyV9D0eP6IioA%3D%3D]

Eric Bogatin. Pop Quiz: When is an Interconnect Not a Transmission Line? [https://www.signalintegrityjournal.com/blogs/4-eric-bogatin-signal-integrity-journal-technical-editor/post/265-pop-quiz-when-is-an-interconnect-not-a-transmission-line]

TeledyneLeCroy/SignalIntegrity Python tools for signal integrity applications [SignalIntegrityApp]

A Look at Transmission-Line Losses [http://blog.teledynelecroy.com/2018/06/a-look-at-transmission-line-losses.html]

How Much Transmission-Line Loss is Too Much? [http://blog.teledynelecroy.com/2018/06/how-much-transmission-line-loss-is-too.html]

Raymond Y. Chen, Raymond Y. Chen. Fundamentals of S Fundamentals of S-Parameter Parameter Modeling for Power Distribution Modeling for Power Distribution System (PDS) and SSO Analysis System (PDS) and SSO Analysis [https://ibis.org/summits/jun05/chen.pdf]

Sam Palermo, ECEN720: High-Speed Links Circuits and Systems Spring 2025 Lecture 9: Noise Sources [https://people.engr.tamu.edu/spalermo/ecen689/lecture9_ee720_noise_sources.pdf]

Basis of Inductance

Eric Bogatin. What Really Is Inductance? [https://speedingedge.com/wp-content/uploads/BTS006_What_Is_Inductance-2.pdf]

Self-Inductance & Mutual Inductance

\[\begin{align} N_a &= L_a I_a + \color{red}M_{ab}I_b \\ N_b &= L_b I_b + \color{red}M_{ab}I_a \end{align}\]

Induced voltage

Partial Inductance

Loop Self & Mutual Inductance

Loop Mutual Inductance

Equivalent Inductance of Multiple Inductors

\[ L_\text{series} = L_1 + L_{12} + L_2 +L_{12} = L_1 + L_2 + 2L_{12} \]

\[ L_\text{parallel} = (L_1 + L_{12})\parallel (L_2 + L_{12}) = \frac{L_1L_2 +L_{12}(L_1+L_2)+L_{12}^2}{L_1+L_2+2L_{12}} \]

internal self-inductance

cross-sectional area

at low frequency

Partial Inductance

Loop Inductance is the sum of partial self-inductance and partial-mutual inductance

Magnetic Vector Potential (磁矢势)

Youjin Deng. 5-3 静磁场的基本规律 [http://staff.ustc.edu.cn/~yjdeng/EM2022/pdf/5-2(2022).pdf]

磁场不能用标量势描述

self partial inductance

mutual partial inductance

Ex. two-wire

[https://www.oldfriend.url.tw/Q3D/ansys_ch_Partial_Loop_Inductance.html]

Chapter 4.5. High Frequency Passive Devices [https://www.cambridge.org/il/files/7713/6698/2369/HFIC_chapter_4_passives.pdf]

On-Chip Spiral Inductors

Designer's Tips on RFIC Inductors [https://www.rfinsights.com/concepts/tips-on-rfic-inductors/]

substrate-induced eddy-current loss

8-Shaped Inductor

J. H. Mikkelsen, O. K. Jensen and T. Larsen, "Crosstalk coupling effects of CMOS co-planar spiral inductors," Proceedings of the IEEE 2004 Custom Integrated Circuits Conference [https://sci-hub.jp/10.1109/CICC.2004.1358825]

P. Guan et al., "8-Shaped Inductors: An Essential Addition to RFIC Designers' Toolbox," in IEEE Open Journal of the Solid-State Circuits Society, vol. 4, pp. 131-146, 2024 [pdf]

• Field picture
• Partial-inductance picture

Transformer

ideal transformer

For an ideal N:1 transformer, \[ \boxed{V_1 = N V_2} \qquad \boxed{I_1 = -\frac{I_2}{N}} \] The minus sign comes from power conservation \[ V_1 I_1 + V_2 I_2 = 0 \] so \[ N V_2 I_1 + V_2 I_2 = 0 \] therefore \[ I_1 = -\frac{I_2}{N} \]

So Fig. 2.3 has exactly the same two-port equations as the original coupled inductors

Electromotive Force (EMF) \(\mathcal{E}\)

back emf

mutual inductance

\[ M_{12}=M_{21}=M \qquad \frac{N_1\Phi_{12}}{i_2}=\frac{N_2\Phi_{21}}{i_1} = M \]

Important: this does not mean \[ \Phi_{12} = \Phi_{21} \]

\(k\) vs. \(M\)

Coupling coefficient \(k\) is based on flux linkage ratios, not directly magnetic flux ratios

任何封闭电路中感应电动势大小，等于穿过这一电路磁通量的变化率。 \[ \epsilon = -\frac{\mathrm{d}\Phi_B}{\mathrm{d}t} \] 其中 \(\epsilon\)是电动势，单位为伏特

\(\Phi_B\)是通过电路的磁通量，单位为韦伯

电动势的方向（公式中的负号）由楞次定律决定

楞次定律: 由于磁通量的改变而产生的感应电流，其方向为抗拒磁通量改变的方向。

在回路中产生感应电动势的原因是由于通过回路平面的磁通量的变化，而不是磁通量本身，即使通过回路的磁通量很大，但只要它不随时间变化，回路中依然不会产生感应电动势。

自感电动势

当电流\(I\)随时间变化时，在线圈中产生的自感电动势为 \[ \epsilon = -L\frac{\mathrm{d}I}{\mathrm{d}t} \]

magnetic flux

magnetic linkage

同名端：当两个电流分别从两个线圈的对应端子流入 ，其所 产生的磁场相互加强时，则这两个对应端子称为同名端。

Integrated Transformer Models [https://www.rfinsights.com/concepts/rfic-transformer-model/]

reference

Bogatin, E. (2018). Signal and power integrity, simplified. Prentice Hall. [pdf]

Paul, Clayton R. Inductance: Loop and Partial. Hoboken, N.J. : [Piscataway, N.J.]: Wiley ; IEEE, 2010. [pdf]

Spartaco Caniggia. Signal Integrity and Radiated Emission of High‐Speed Digital Systems. Wiley 2008

Luong, H. C., & Yin, J. (2016). Transformer-based design techniques for oscillators and frequency dividers. Springer International Publishing

ISSCC2002. Special Topic Evening Discussion Sessions SE1: Inductance: Implications and Solutions for High-Speed Digital Circuits [vSE1_Blaauw], [vSE1_Gauthier], [vSE1_Morton, [vSE1_Restle]]

Y. Massoud and Y. Ismail, "Gasping the impact of on-chip inductance," in IEEE Circuits and Devices Magazine, vol. 17, no. 4, pp. 14-21, July 2001 [https://sci-hub.se/10.1109/101.950046]

Clayton R. Paul, Partial Inductance [https://ewh.ieee.org/soc/emcs/acstrial/newsletters/summer10/PP_PartialInductance.pdf]

Cheung-Wei Lam. Common Misconceptions about Inductance & Current Return Path [https://ewh.ieee.org/r6/scv/emc/archive/022010Lam.pdf]

Randy Wolff. Signal Loop Inductance in [Pin] and [Package Model] [https://ibis.org/summits/feb10/wolff.pdf]

ANSYS Q3D Getting Started LE05. Module 5: Q3D Inductance Matrix Reduction [https://innovationspace.ansys.com/courses/wp-content/uploads/sites/5/2021/07/Q3D_GS_2020R1_EN_LE05_Ind_Matrix.pdf]

magnetic field (magnetic flux density, \(B\)), is the tesla (symbol: \(T\)), defined as one weber per square meter (\(Wb/m^2\))

Magnetic flux \(\Phi_B\) measures the total magnetic field (\(B\)) passing through a given surface area (\(A\)), representing the number of field lines penetrating that area. Measured in Webers (\(Wb\)),

Vector Calculus

3Blue1Brown, Divergence and curl: The language of Maxwell's equations, fluid flow, and more [https://youtu.be/rB83DpBJQsE]

Gradient

Divergence

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

# https://share.google/aimode/l3lNa2MRAOG8hkpOc

import numpy as np
import matplotlib.pyplot as plt

# 1. Define the grid
x = np.linspace(-6, 6, 40)
y = np.linspace(-3, 3, 20)
X, Y = np.meshgrid(x, y)

# 2. Define the vector field components F = [U, V]
U = np.sin(X)
V = np.cos(Y)

# 3. Calculate Divergence (Scalar Field)
div_F = np.cos(X) - np.sin(Y)

# 4. Create the plot
fig, ax = plt.subplots(figsize=(16, 8))

# Use 'RdBu_r' (reversed) so Positive = Red, Negative = Blue
contour = ax.contourf(X, Y, div_F, cmap='RdBu_r', levels=30, alpha=0.8)
fig.colorbar(contour, label='Divergence (Red=Source, Blue=Sink)')

# Plot Vector Field
ax.quiver(X, Y, U, V, color='black', alpha=0.9, scale=20)
ax.set_xlabel('X', fontsize=12)
ax.set_ylabel('Y', fontsize=12)
ax.set_title(r'Positive Divergence (Red) and Negative Divergence (Blue)')
plt.show()

Curl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36

# https://share.google/aimode/hWb3cR4vCBoWV4Moi

import numpy as np
import matplotlib.pyplot as plt

# 1. Setup the coordinate grid
x = np.linspace(-2, 2, 7)
y = np.linspace(-2, 2, 7)
z = np.linspace(-2, 2, 7)
X, Y, Z = np.meshgrid(x, y, z)

# 2. Define the Vector Field F and Curl(F)
U, V, W = X, Y*Z, 3*X*Z
C_U, C_V, C_W = -Y, -3*Z, np.zeros_like(Z)

# 3. Create the plot
fig = plt.figure(figsize=(16, 8), constrained_layout=True)

# Subplot 1: Vector Field
ax1 = fig.add_subplot(121, projection='3d')
ax1.quiver(X, Y, Z, U, V, W, length=0.3, normalize=True, color='royalblue')
ax1.set_title('Vector Field F')
ax1.set_xlabel('X', fontsize=16)  # Adding x-label
ax1.set_ylabel('Y', fontsize=16)  # Adding y-label
ax1.set_zlabel('Z', fontsize=16)  # Adding z-label

# Subplot 2: Curl
ax2 = fig.add_subplot(122, projection='3d')
ax2.quiver(X, Y, Z, C_U, C_V, C_W, length=0.3, normalize=True, color='crimson')
ax2.set_title('Curl(F)')
ax2.set_xlabel('X', fontsize=16)  # Adding x-label
ax2.set_ylabel('Y', fontsize=16)  # Adding y-label
ax2.set_zlabel('Z', fontsize=16)  # Adding z-label

# plt.tight_layout(pad=3.0, rect=[0, 0, 1, 0.95])
plt.show()

cross product [Google AI Mode]

Divergence Theorem (Gauss's Theorem)

surface integral -> volume integral \[ \oint_S \vec{F} \cdot d\vec{S} = \int_V (\nabla \cdot \vec{F}) dv \]

Divergence theorem is only applicable to closed surfaces

Stoke's theorem

line integral -> surface inegral \[ \oint_{c} \vec{F} \cdot \vec{\mathrm{d}l} = \int_{s} (\nabla \times \vec{F}) \cdot \vec{\mathrm{d}S} \]

Polarization \(\mathbf P\) & magnetization \(\mathbf M\)

Electric Fields in Matter & Magnetic Fields in Matter

TODO 📅

Electric Fields

Electric field intensity \(\mathbf E\) & Electric flux density \(\mathbf D\)

\[ \boxed{ \underbrace{\varepsilon_0\mathbf E}_{\text{contains free + bound charge effects}} + \underbrace{\mathbf P}_{\text{cancels bound-charge divergence}} = \underbrace{\mathbf D}_{\text{free-charge Gauss law}} }. \] Start from Gauss’s law for the actual electric field: \[ \nabla\cdot(\varepsilon_0\mathbf E) = \rho_{\text{free}}+\rho_{\text{bound}}. \] Polarization satisfies \[ \rho_{\text{bound}}=-\nabla\cdot\mathbf P. \] Therefore, \[ \nabla\cdot(\varepsilon_0\mathbf E) = \rho_{\text{free}}-\nabla\cdot\mathbf P. \] Move the polarization term to the left: \[ \nabla\cdot(\varepsilon_0\mathbf E+\mathbf P) = \rho_{\text{free}}. \] Define \[ \mathbf D=\varepsilon_0\mathbf E+\mathbf P. \]

Electric field intensity \(\mathbf E\)

Gauss's law for \(\mathbf E\) is \[ \boxed{ \nabla\cdot\mathbf E = \frac{\rho_{\text{total}}}{\varepsilon_0} } \] where \[ \rho_{\text{total}} = \rho_{\text{free}}+\rho_{\text{bound}}. \] Therefore, \(\mathbf E\) is determined by all charges:

• externally supplied free charges,
• polarization-induced bound charges.

In integral form, \[ \oint_S \mathbf E\cdot d\mathbf S = \frac{Q_{\text{total,enclosed}}}{\varepsilon_0}. \]

Electric flux density \(\mathbf D\)

Gauss's law for \(\mathbf D\) is \[ \boxed{ \nabla\cdot\mathbf D=\rho_{\text{free}} } \] or \[ \boxed{ \oint_S\mathbf D\cdot d\mathbf S = Q_{\text{free,enclosed}} } \] Thus, \(\mathbf D\) is constructed so that dielectric bound charge is absorbed into the constitutive relation \[ \mathbf D=\varepsilon_0\mathbf E+\mathbf P. \] It therefore relates directly only to free charge

Gauss's law for \(\mathbf{E}\) & Gauss's law for \(\mathbf{D}\)

\[ \boxed{ \begin{aligned} \nabla\cdot\mathbf E &=\frac{\rho_{\text{total}}}{\varepsilon_0}, \\[4pt] \nabla\cdot\mathbf D &=\rho_{\text{free}}. \end{aligned} } \]

So when a textbook writes \[ \oint_S\mathbf D\cdot d\mathbf S = \int_V\rho_v\,dV, \] the symbol \(\rho_v\) normally means \[ \boxed{\rho_v=\rho_{v,\text{free}}}. \] But when it writes \[ \oint_S\mathbf E\cdot d\mathbf S = \frac{1}{\varepsilon_0}\int_V\rho_v\,dV, \] then \(\rho_v\) means the total volume-charge density, unless the context explicitly assumes vacuum or no polarization

Magnetostatics

Magnetostatics is the study of magnetic fields in systems where the currents are steady (not changing with time)

Relationship between \(\mathbf B\) and \(\mathbf H\)

The central distinction is \[ \boxed{ \mathbf H\text{ tracks the free-current excitation, while } \mathbf B\text{ is the resulting total magnetic flux density.} } \]

The general macroscopic relation is \[ \boxed{ \mathbf B=\mu_0\left(\mathbf H+\mathbf M\right) } \] where

• \(\mathbf M\) is the magnetization density of the material,
• \(\mu_0\mathbf H\) represents the part separated from material magnetization,
• \(\mu_0\mathbf M\) is the contribution caused by magnetic dipoles inside the material.

Therefore, \[ \boxed{ \mathbf H=\frac{\mathbf B}{\mu_0}-\mathbf M }. \] In vacuum

There is no material magnetization: \[ \mathbf M=0. \] Hence, \[ \boxed{\mathbf B=\mu_0\mathbf H}. \] In a linear, isotropic material

If \[ \mathbf M=\chi_m\mathbf H, \] then \[ \mathbf B = \mu_0(1+\chi_m)\mathbf H. \] Define \[ \mu_r=1+\chi_m, \qquad \mu=\mu_0\mu_r. \] Then \[ \boxed{\mathbf B=\mu\mathbf H}. \]

Suppose a coil carries a free current \(I\).

The free current produces an \(\mathbf H\) field. If a magnetic material is placed inside the coil, the material becomes magnetized: \[ I_{\mathrm{free}} \longrightarrow \mathbf H \longrightarrow \mathbf M. \] The resulting total magnetic flux density is \[ \mathbf B=\mu_0(\mathbf H+\mathbf M). \] Therefore: \[ \boxed{ \mathbf H \text{ is related to the applied free current} } \] while \[ \boxed{ \mathbf B \text{ includes both the applied field and the material response}. } \] For the same coil current, \(\mathbf H\) may remain approximately the same, but inserting a high-permeability core can greatly increase \(\mathbf B\).

Magnetic field intensity \(\mathbf H\)

It is introduced to describe magnetic fields in matter while separating the material’s magnetization from externally supplied, or free, currents.

The Maxwell–Ampère law is \[ \boxed{ \nabla\times\mathbf H = \mathbf J_{\mathrm{free}} + \frac{\partial\mathbf D}{\partial t} } \] or, in integral form, \[ \boxed{ \oint_C\mathbf H\cdot d\mathbf l = I_{\mathrm{free}} + \int_S\frac{\partial\mathbf D}{\partial t}\cdot d\mathbf S }. \] Thus, the circulation of \(\mathbf H\) is associated with free current.

Magnetic Flux Density \(\mathbf B\)

\(\mathbf B\) determines the magnetic force on a moving charge.

It also defines magnetic flux: \[ \boxed{\Phi_B=\int_S\mathbf B\cdot d\mathbf S} \] and appears in Faraday's law: \[ \oint_C\mathbf E\cdot d\mathbf l = -\frac{d}{dt} \int_S\mathbf B\cdot d\mathbf S. \] Gauss's law for magnetism is \[ \boxed{\nabla\cdot\mathbf B=0} \] or \[ \boxed{\oint_S\mathbf B\cdot d\mathbf S=0}. \] This is universally valid because magnetic monopoles have not been observed.

Electrodynamics

Faraday's Law

Lenz's law

Ampère's law w/ Maxwell's correction

[https://youtu.be/t-W6_fn_1bI]

\[ \boxed{ \begin{aligned} \oint\mathbf H\cdot d\mathbf l &= I+\int\frac{\partial\mathbf D}{\partial t}\cdot d\mathbf S, \\[4pt] \nabla\times\mathbf B &= \mu_0\mathbf J+ \mu_0\varepsilon_0\frac{\partial\mathbf E}{\partial t} \end{aligned} } \] are equivalent in vacuum \[ \boxed{ \begin{aligned} \mathbf H,\mathbf D\text{ form} &\rightarrow \mathbf J_{\mathrm{free}},\\ \mathbf B,\mathbf E\text{ form} &\rightarrow \mathbf J_{\mathrm{total}}. \end{aligned} } \]

If \[ \mathbf B=\mu\mathbf H, \qquad \mathbf D=\varepsilon\mathbf E, \] with constant \(\mu\) and \(\varepsilon\), then \[ \boxed{ \oint_C\mathbf B\cdot d\mathbf l = \mu I_{\mathrm{free}} + \mu\varepsilon \int_S \frac{\partial\mathbf E}{\partial t}\cdot d\mathbf S } \]

Gauss’s law for magnetism

The universally valid equation is \[ \boxed{\oint_S \mathbf B\cdot d\mathbf S=0} \] or \[ \nabla\cdot\mathbf B=0. \] This expresses that magnetic field lines have no beginning or end—there are no observed magnetic monopoles

\[ \boxed{ \oint_S\mathbf B\cdot d\mathbf S=0 \text{ always, while } \oint_S\mathbf H\cdot d\mathbf S=0 \text{ only under additional conditions.} } \] i.e. \(\mathbf B=\mu\mathbf H\)

Field Inside and Outside a Current-Carrying Wire

Sources of Magnetic Fields [https://web.mit.edu/8.02t/www/802TEAL3D/visualizations/coursenotes/modules/guide09.pdf]

Displacement Current

A. Sheikholeslami, "Current Without Electrons [Circuit Intuitions]," in IEEE Solid-State Circuits Magazine, vol. 17, no. 4, pp. 8-10, Fall 2025

—, "Current Without Electric Field [Circuit Intuitions]," in IEEE Solid-State Circuits Magazine, vol. 18, no. 1, pp. 8-12, winter 2026

energy and information are carried by electric and magnetic fields (\(E\) and \(H\)) rather than by electron drift

proximity effect & skin effect

• Skin effect concentrates current near the surface of a single conductor, while proximity effect concentrates current in specific regions of multiple conductors due to their interaction
• Skin effect is caused by the conductor's own magnetic field, while proximity effect is caused by the magnetic field of a nearby conductor

proximity effect is a redistribution of electric current occurring in nearby parallel electrical conductors carrying alternating current (AC), caused by magnetic effects (eddy currents)

skin effect is the tendency of AC current flow near the surface (or "skin") of a conductor, rather than throughout its cross-section, due to the magnetic field generated by the current itself

Cause of skin effect

A main current \(I\) flowing through a conductor induces a magnetic field \(H\). If the current increases, as in this figure, the resulting increase in \(H\) induces separate, circulating eddy currents \(I_W\) which partially cancel the current flow in the center and reinforce it near the skin

Eddy current

By Lenz's law, an eddy current creates a magnetic field that opposes the change in the magnetic field that created it, and thus eddy currents react back on the source of the magnetic field

reference

Griffiths, David J. Introduction to Electrodynamics. Fifth edition. Cambridge University Press, 2024. [pdf]

David Smith, Electromagnetic Theory for Complete Idiot, 2021

邓友金. 电磁学 2022春 [http://staff.ustc.edu.cn/~yjdeng/EM2022/EM2022.html]

谢处方、饶克谨、杨显清等.《电磁场与电磁波》（第五版），高等教育出. 版社，2019.

Scott Hughes. Spring 2005 8.022: Electricity & Magnetism [https://web.mit.edu/sahughes/www/8.022/]

Aditya Varma Muppala, EE 210 - Applied Electromagnetic Theory [https://adityamuppala.github.io/teaching210/]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

# Parameters: 4th order filter with 1kHz cutoff
order = 4
w_cutoff = 2 * np.pi * 1000
w = np.logspace(2, 5, 1000) # Frequency range from 100Hz to 100kHz

# Generate Butterworth coefficients and frequency response
b_but, a_but = signal.butter(order, w_cutoff, analog=True)
w_but, h_but = signal.freqs(b_but, a_but, worN=w)

# Generate Chebyshev (1dB ripple) response
b_cheb, a_cheb = signal.cheby1(order, 1, w_cutoff, analog=True)
w_cheb, h_cheb = signal.freqs(b_cheb, a_cheb, worN=w)

# Generate Bessel response
b_bess, a_bess = signal.bessel(order, w_cutoff, analog=True)
w_bess, h_bess = signal.freqs(b_bess, a_bess, worN=w)

freqs = w / (2 * np.pi)

t = np.linspace(0, 0.005, 1000) # 5ms window

# Calculate Step Responses
t_but, y_but = signal.step((b_but, a_but), T=t)
t_cheb, y_cheb = signal.step((b_cheb, a_cheb), T=t)
t_bess, y_bess = signal.step((b_bess, a_bess), T=t)

Butterworth Filters

Butterworth Filters [https://people.eecs.ku.edu/~demarest/212/Butterworth%20Filters.pdf]

Stephen Roberts, Signal Processing & Filter Design B3 option: Lecture 2 - Frequency Selective Filters [https://www.robots.ox.ac.uk/~sjrob/Teaching/SP/l2.pdf]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42

% Parameters
wc = 100;               % Cutoff frequency in rad/s
orders = [1, 2, 6];     % Orders to compare
w = logspace(0, 3, 1000); % Frequency range (1 to 1000 rad/s)

for n = orders
    % 1. Design Analog Filter (the 's' flag is critical)
    [b, a] = butter(n, wc, 's');
    
    % 2. Calculate Complex Frequency Response
    h = freqs(b, a, w);
    
    % 3. Calculate Group Delay
    % For analog filters, we manually differentiate phase: gd = -d(phi)/dw
    phase = unwrap(angle(h));
    gd = -diff(phase) ./ diff(w);
    
    % --- Plot Magnitude (dB) ---
    subplot(3,1,1);
    semilogx(w, 20*log10(abs(h)), 'LineWidth', 2, 'DisplayName', ['n=' num2str(n)]);
    hold on;
    
    % --- Plot Phase (Degrees) ---
    subplot(3,1,2);
    semilogx(w, phase * 180/pi, 'LineWidth', 2);
    hold on;
    
    % --- Plot Group Delay ---
    subplot(3,1,3);
    semilogx(w(1:end-1), gd, 'LineWidth', 2);
    hold on;
end

% Formatting
subplot(3,1,1)
ylabel('Magnitude (dB)'); grid on; xline(wc, '--r', 'HandleVisibility','off'); legend('show');
ylim([-60 5]);

subplot(3,1,2); ylabel('Phase (deg)'); grid on; xline(wc, '--r', 'HandleVisibility','off');

subplot(3,1,3); ylabel('Group Delay (sec)'); xlabel('Frequency (rad/s)'); 
grid on; xline(wc, '--r', 'HandleVisibility','off');

Bessel Filters

Stephen Roberts, Signal Processing & Filter Design B3 option: Lecture 3 - Transient Response and Transforms [https://www.robots.ox.ac.uk/~sjrob/Teaching/SP/l3.pdf]

Bessel filter is often called Bessel–Thomson filter or simply Thomson filter

besself: Bessel analog filter design [b,a] = besself(n, Wo)

the transfer function coefficients of an \(n\)th-order lowpass analog Bessel filter, where Wo is the angular frequency up to which the filter's group delay is approximately constant. Larger values of n produce a group delay that better approximates a constant up to Wo.

scipy.signal.bessel(N, Wn, btype='low', analog=True, output='ba', norm='phase')

norm='phase' — The filter is normalized such that the phase response reaches its midpoint at angular (e.g. rad/s) frequency Wn

This is the default, and matches MATLAB's implementation.

1
2
3
4
5

octave:8> [b,a] = besself(5,1)
b = 1
a =

   1.0000   3.8107   6.7767   6.8864   3.9363   1.0000

1
2
3
4
5

b_bess, a_bess = signal.bessel(5, 1, btype='low', analog=True, norm='phase')
print(b_bess, a_bess)

# [1.] [1.         3.81070121 6.77667372 6.88636765 3.93628343 1.        ]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

order = 5
cutoff_rad = 1

b_bess, a_bess = signal.bessel(order, cutoff_rad, btype='low', analog=True, norm='phase')

# Frequency response computation (Linear scale)
w = np.linspace(0, 8 * cutoff_rad, 1000) # Frequency range 0 to 2*fc

# Frequency response for phase
w_bess, h_bess = signal.freqs(b_bess, a_bess, worN=w)

# Phase Response (unwrapped to avoid 360 degree jumps)
phase_bess = np.unwrap(np.angle(h_bess))

plt.figure(figsize=(10, 5))
plt.plot(w_bess, np.degrees(phase_bess), label='Bessel', linewidth=2, linestyle='-')
plt.axvline(cutoff_rad, color='red', linestyle=':', label='Cutoff')
plt.title('5th Order Bessel Filter Phase Response w/ Wn=1')
plt.ylabel('Phase [degrees]'); plt.xlabel('Frequency [rad/s]')
plt.grid(True); plt.legend()

1
2
3
4
5
6
7

Wn = 1;

[bb, ab] = butter(2, Wn, 's');
[bs, as] = besself(2, Wn);

pb = roots(ab);
ps = roots(as);

RC LPF

Kwantae Kim, Integrated Analog Systems D - Lecture 02 (Continuous-Time Filters) [https://youtu.be/B7-kr5zV3NA]

—, Integrated Analog Systems D - Lecture 03 (Continuous-Time Filters) [https://youtu.be/6GdDiwaKDZw]

—, Integrated Analog Systems D - Lecture 05 (Continuous-Time Filters) [https://youtu.be/LHhEK1RlC6w]

[Gist link]

2nd-Order RC LPF

Loading effect & limitation

Phase Magin with damping Factor \(\zeta\)

\[ \boxed{\phi_\text{PM}\approx 100\cdot \zeta} \]

General 2nd-Order RC LPF

Laplace Transform

Stability Analysis

Active Filters

Kwantae Kim, Integrated Analog Systems D - Lecture 05 (Continuous-Time Filters) [https://youtu.be/LHhEK1RlC6w]

[Gist link]

TODO 📅

Single-Pole LPF Algorithms

Neil Robertson, Model a Sigma-Delta DAC Plus RC Filter [https://www.dsprelated.com/showarticle/1642.php]

Jason Sachs, Ten Little Algorithms, Part 2: The Single-Pole Low-Pass Filter [https://www.embeddedrelated.com/showarticle/779.php]

—. Return of the Delta-Sigma Modulators, Part 1: Modulation [https://www.dsprelated.com/showarticle/1517/return-of-the-delta-sigma-modulators-part-1-modulation]

• Derivatives Approximation (\(H_p(s)=\frac{1}{s\tau +1}\))

  \[\begin{align} H_p(z)&=\frac{\frac{T_s}{T_s+\tau}}{1+(\frac{T_s}{T_s+\tau}-1)z^{-1}}\tag{EQ-0}\\ H_p(z)&=\frac{\frac{T_s}{\tau}}{1+(\frac{T_s}{\tau}-1)z^{-1}}\tag{EQ-1} \end{align}\]

• Matched z-Transform (Root Matching) \[ H_p(z)=\frac{1-e^{-T_s/\tau}}{1-e^{-T_s/\tau}z^{-1}}\tag{EQ-2} \] EQ-2 is connected with EQ-1 by \(1 - e^{-\Delta t/\tau} \approx \frac{\Delta t}{\tau}\)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

import numpy as np
import matplotlib.pyplot as plt
import scipy.signal


dt=0.002; tau=0.05

T=1
t = np.arange(0,T+1e-5,dt)
x = 1-np.abs(4*t-1)
x[4*t>2] = 0.95; x[t>0.75] = 0.02

alpha_derv = dt / (dt + tau)
alpha_derv_approx = dt / tau

yfilt_derv = scipy.signal.lfilter([alpha_derv], [1, alpha_derv - 1], x)
yfilt_derv_approx = scipy.signal.lfilter([alpha_derv_approx], [1, alpha_derv_approx - 1], x)

a1 = -np.exp(-dt/tau)
b0 = [1 + a1]
a = [1, a1]
y_filt_match = scipy.signal.lfilter(b0, a, x)

plt.figure(figsize=(20,10))
plt.plot(t, x, color='k', linewidth=3, label='x')
plt.plot(t, yfilt_derv, color='red', linewidth=3, label=r'$H_p(z)=\frac{\frac{T_s}{T_s+\tau}}{1+(\frac{T_s}{T_s+\tau}-1)z^{-1}}$')
plt.plot(t, yfilt_derv_approx, color='green', marker='D', linestyle='dashed',markersize=4, linewidth=3, label=r'$H_p(z)=\frac{\frac{T_s}{\tau}}{1+(\frac{T_s}{\tau}-1)z^{-1}}$')
plt.plot(t, y_filt_match, color='m', marker='x', linestyle='dashed',markersize=4, linewidth=3, label=r'$H_p(z)=\frac{1-e^{-T_s/\tau}}{1-e^{-T_s/\tau}z^{-1}}$')

plt.legend(loc='upper right', fontsize=20)
plt.grid(which='both');plt.xlabel('Time (s)');plt.show()