KNOW-HOW

One will see the following phenomena:

1. A power spectrum is equal to the square of the absolute value of DFT.
2. The sum of all power spectral lines in a power spectrum is equal to the power of the input signal.
3. The integral of a PSD is equal to the power of the input signal.

power spectrum has units of \(V^2\) and power spectral density has units of \(V^2/Hz\)

Parseval's theorem is a property of the Discrete Fourier Transform (DFT) that states: \[ \sum_{n=0}^{N-1}|x(n)|^2 = \frac{1}{N}\sum_{k=0}^{N-1}|X(k)|^2 \] Multiply both sides of the above by \(1/N\): \[ \frac{1}{N}\sum_{n=0}^{N-1}|x(n)|^2 = \frac{1}{N^2}\sum_{k=0}^{N-1}|X(k)|^2 \] \(|x(n)|^2\) is instantaneous power of a sample of the time signal. So the left side of the equation is just the average power of the signal over the N samples. \[ P_{\text{av}} = \frac{1}{N^2}\sum_{k=0}^{N-1}|X(k)|^2\text{, }V^2 \] For the each DFT bin, we can say: \[ P_{\text{bin}}(k) = \frac{1}{N^2}|X(k)|^2\text{, k=0:N-1, }V^2/\text{bin} \] This is the power spectrum of the signal.

Note that \(X(k)\) is the two-sided spectrum. If \(x(n)\) is real, then \(X(k)\) is symmetric about \(fs/2\), with each side containing half of the power. In that case, we can choose to keep just the one-sided spectrum, and multiply Pbin by 2 (except DC & Nyquist):

\[ P_{\text{bin,one-sided}}(k) = \left\{ \begin{array}{cl} \frac{1}{N^2}|X(0)|^2 & : \ k = 0 \\ \frac{2}{N^2}|X(k)|^2 & : \ 1 \leq k \leq N/2-1 \\ \frac{1}{N^2}|X(N/2)|^2 & : \ k = N/2 \end{array} \right. \]

1
2
3
4
5
6
7
8

rng default
Fs = 1000;
t = 0:1/Fs:1-1/Fs;
x = cos(2*pi*100*t) + randn(size(t));
N = length(x);
xdft = fft(x);
xsq_sum_avg = sum(x.^2)/N;
specsq_sum_avg = sum(abs(xdft).^2)/N^2;

where xsq_sum_avg is same with specsq_sum_avg

For a discrete-time sequence x(n), the DFT is defined as: \[ X(k) = \sum_{n=0}^{N-1}x(n)e^{-j2\pi kn/N} \] By it definition, the DFT does NOT apply to infinite duration signals.

scaling DFT results

Tutorial on Scaling of the Discrete Fourier Transform and the Implied Physical Units of the Spectra of Time-Discrete Signals Jens Ahrens, Carl Andersson, Patrik Höstmad, Wolfgang Kropp

[https://github.com/AppliedAcousticsChalmers/scaling-of-the-dft]

Different scaling is needed to apply for amplitude spectrum, power spectrum and power spectrum density, which shown as below

\(f_s\) in Eq.(13) is sample rate rather than frequency resolution.

And Eq.(13) can be expressed as \[ \text{PSD}(k) =\frac{1}{f_{\text{res}}\cdot N\sum_{n}w^2(n)}\left|X_{\omega}(k)\right|^2 \] where \(f_{\text{res}}\) is frequency resolution

We define the following two sums for normalization purposes:

\[\begin{align} S_1 &= \sum_{j=0}^{N-1}w_j \\ S_2 &= \sum_{j=0}^{N-1}w_j^2 \end{align}\]

Given Eq.(12) and Eq.(13) \[\begin{align} \text{PS(k)} &= \text{PSD(k)} \cdot \frac{f_s \sum_{n}w^2(n)}{(\sum_{n}\omega(n))^2} \\ &= \text{PSD(k)} \cdot f_{\text{res}} \cdot \frac{N \sum_{n}w^2(n)}{(\sum_{n}w(n))^2} \\ &= \text{PSD(k)} \cdot f_{\text{res}} \cdot \frac{N S_2}{S_1^2} \\ &= \text{PSD(k)} \cdot f_{\text{res}} \cdot \text{NENBW} \\ &= \text{PSD(k)} \cdot \text{ENBW} \end{align}\]

where Normalized Equivalent Noise BandWidth is defined as \[ \text{NENBW} =\frac{N S_2}{S_1^2} \] and Effective Noise BandWidth is \[ \text{ENBW} =f_{\text{res}} \cdot \frac{N S_2}{S_1^2} \]

For Rectangular window, \(\text{ENBW} =f_{\text{res}}\)

This equivalent noise bandwidth is required when the resulting spectrum is to be expressed as spectral density (such as for noise measurements).

Eq.(13) Example:

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

Fs = 1024; 
t  = (0:1/Fs:1-1/Fs).'; 
x  = sin(2*pi*t*200);
Nx = length(x);
% Window data
w = hamming(Nx);	% hamming window
xw = x.*w; 
% Calculate power
nfft = Nx; 
X = fft(xw,nfft);

NumUniquePts = nfft/2+1;
Pxx1 = abs(X(1:NumUniquePts)).^2/sum(w.^2)/Fs;    %!!! Eq.(13)
Pxx1(2:end-1) = Pxx1(2:end-1)*2;    % twosided -> onesided
Fx1 = (0:NumUniquePts-1)*Fs/nfft; 

[Pxx2,Fx2] = pwelch(x,w,0,nfft,Fs);

plot(Fx1,10*log10(Pxx1),Fx2,10*log10(Pxx2),'r--');
legend('PSD via Eq.(13)','PSD via pwelch')

window effects

It is possible to correct both the amplitude and energy content of the windowed signal to equal the original signal. However, both corrections cannot be applied simultaneously

Amplitude correction

Linear (amplitude) Spectrum \[ X(k) = \frac{X_{\omega}(k)}{S_1} \] power spectrum \[ \text{PS}=\frac{\left| X_{\omega}(k) \right|^2}{S_1^2} \]

usage: tone amplitude

Power correction

power spectral density (PSD) \[ \text{PSD} =\frac{\left|X_{\omega}(k)\right|^2}{f_s\cdot S_2} \]

We have \(\text{PSD} = \frac{\text{PS}}{\text{ENBW}}\), where \(\text{ENBW}=\frac{N \cdot S_2}{S_1^2}f_{\text{res}}\)

linear power spectrum \[ \text{PS}_L=\frac{|X_{\omega}(k)|^2}{N\cdot S_2} \]

usage: RMS value, total power \[ \text{PS}_L(k)=\text{PSD(k)} \cdot f_{\text{res}} \]

Window Correction Factors

While a window helps reduce leakage (The window reduces the jumps at the ends of the repeated signal), the window itself distorts the data in two different ways:

• Amplitude – The amplitude of the signal is reduced

  This is due to the fact that the window removes information in the signal

• Energy – The area under the curve, or energy of the signal, is reduced

Window correction factors are used to try and compensate for the effects of applying a window to data. There are both amplitude and energy correction factors.

Only the Uniform window (rectangular window), which is equivalent to no window, has the same amplitude and energy correction factors.

\[\begin{align} X(k) &= \frac{X_\omega(k)}{\sum_n w(n)} \\ &= \frac{N}{\sum_n w(n)} \frac{X_\omega(k)}{N} \\ &= K_a \cdot \frac{X_\omega(k)}{N} \end{align}\] where \(K_a = \frac{N}{\sum_n \omega(n)}\)

In literature, Coherent power gain is defined show below, which is close related to \(K_a\) \[ \text{Coherent power gain (dB)} = 20 \; log_{10} \left( \frac{\sum_n w[n]}{N} \right) \]

With amplitude correction, by multiplying by two, the peak value of both the original and corrected spectrum match. However the energy content is not the same.

The amplitude corrected signal (red) appears to have more energy, or area under the curve, than the original signal (blue).

\[\begin{align} \text{PSD} &=\frac{1}{f_{\text{res}}\cdot N\sum_{n}w^2(n)}\left|X_{\omega}(K)\right|^2 \\ &= \frac{N}{\sum_{n}w^2(n)} \frac{\left|X_{\omega}(K)\right|^2}{f_{\text{res}}\cdot N^2} \\ &=\left|K_e \cdot \frac{X_{\omega}(K)}{N}\right|^2 \cdot \frac{1}{f_{\text{res}}} \end{align}\]

wherer \(K_e = \sqrt{\frac{N}{\sum_{n}\omega^2(n)}}\)

Multiplying the values in the spectrum by 1.63, rather than 2, makes the area under the curve the same for both the original signal (blue) and energy corrected signal (red)

hanning's correction factors:

1
2
3
4

N = 256;
w = hanning(N);
Ka = N/sum(w)
Ke = sqrt(N/sum(w.^2))

1
2
3
4
5
6
7
8

Ka =

    1.9922


Ke =

    1.6298

Amplitude,Power Correction Example

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

%% 200 Hz, amplitude 1 and 200/3 Hz amplitude 0.3
fosc = 200;
Fs = 8192;
Nx = Fs*4;
tstop = Nx/Fs;
t  = (0:1/Fs:tstop-1/Fs).'; 
x  = sin(2*pi*t*fosc) + 0.3*sin(2*pi*t*fosc/3);
Nx = length(x);

%% direct method, time domain
P0 = sum(x.^2)/Nx;
xrms0 = sqrt(P0);   % RMS value caculated by time domain
fprintf("xrms0: %.4f\n", xrms0);

%% Window data
w = hamming(Nx);	% hamming window
%w = ones(Nx, 1);	% rectangular window

S1= sum(w);
S2 = sum(w.^2);
NENBW = Nx*S2/S1^2;

xw = x.*w; 
nfft = Nx; 
X = fft(xw,nfft);   % windowed signal fft
NumUniquePts = nfft/2+1;
Fx1 = (0:NumUniquePts-1)*Fs/nfft;

%% fft - ps -> rms
Xmag = abs(X);
PS1 = Xmag.^2/S1^2;    %!!! Eq.(12)
PS1 = PS1(1:NumUniquePts);
PS1(2:end-1) = PS1(2:end-1)*2;    % twosided -> onesided
Pmax_fft = max(PS1);    % !!! amplitude OK
fprintf("Pmax_fft: %.4f\n", Pmax_fft);

% linear power spectrum
PSL_fft = Xmag.^2/nfft/S2;
PSL_fft = PSL_fft(1:NumUniquePts);
PSL_fft = PSL_fft(2:end-1)*2;
Ptot_fft1 = sum(PSL_fft);
xrms_PSL = sqrt(Ptot_fft1); % !!! power OK
fprintf("xrms_PSL: %.4f\n", xrms_PSL);

% use psd*fres to calculate total power
Pxx1 = abs(X(1:NumUniquePts)).^2/sum(w.^2)/Fs;    %!!! Eq.(13)
Pxx1(2:end-1) = Pxx1(2:end-1)*2;    % twosided -> onesided
Ptot_fft2 = sum(Pxx1*Fs/nfft);
xrms_psd = sqrt(Ptot_fft2);     % !!! power OK
fprintf("xrms_psd: %.4f\n", xrms_psd);


%% pwelch
[Pxx2,Fx2] = pwelch(x,w,0,nfft,Fs);
fres = Fx2(2)-Fx2(1);

PS2 = Pxx2*fres;
Ptot_pwelch = sum(PS2);
xrms_pwelch = sqrt(Ptot_pwelch);    % !!! amplitude OK
fprintf("xrms_pwelch: %.4f\n", xrms_pwelch);

% PS = PSD * ENBW always hold
ENBW = fres * NENBW;
PS2_pwelch = Pxx2*ENBW;
Pmax_pwelch = max(PS2_pwelch);  % !!! amplitude OK
fprintf("Pmax_pwelch: %.4f\n", Pmax_pwelch);

%% plot psd of two methods
plot(Fx1,10*log10(Pxx1),Fx2,10*log10(Pxx2),'r--');
legend('PSD via Eq.(13)','PSD via pwelch')
grid on;
xlabel('Freq (Hz)');
ylabel('dB; V^2/Hz')

1
2
3
4
5
6

xrms0: 0.7383
Pmax_fft: 0.5000
xrms_PSL: 0.7382
xrms_psd: 0.7382
xrms_pwelch: 0.7382
Pmax_pwelch: 0.5000

RMS

We may also want to know for example the RMS value of the signal, in order to know how much power the signal generates. This can be done using Parseval’s theorem.

For a periodic signal, which has a discrete spectrum, we obtain its total RMS value by summing the included signals using \[ x_{\text{rms}}=\sqrt{\sum R_{xk}^2} \] Where \(R_{xk}\) is the RMS value of each sinusoid for \(k=1,2,3,...\) The RMS value of a signal consisting of a number of sinusoids is consequently equal to the square root of the sum of the RMS values.

This result could also be explained by noting that sinusoids of different frequencies are orthogonal, and can therefore be summed like vectors (using Pythagoras’ theorem)

For a random signal we cannot interpret the spectrum in the same way. As we have stated earlier, the PSD of a random signal contains all frequencies in a particular frequency band, which makes it impossible to add the frequencies up. Instead, as the PSD is a density function, the correct interpretation is to sum the area under the PSD in a specific frequency range, which then is the square of the RMS, i.e., the mean-square value of the signal \[ x_{\text{rms}}=\sqrt{\int G_{xx}(f)df}=\sqrt{\text{area under the curve}} \] The linear spectrum, or RMS spectrum, defined by \[\begin{align} X_L(k) &= \sqrt{\text{PSD(k)} \cdot f_{\text{res}}}\\ &=\sqrt{\frac{\left|X_{\omega}(k)\right|^2}{f_{\text{res}}\cdot N\sum_{n}\omega^2(n)} \cdot f_{\text{res}}} \\ &= \sqrt{\frac{\left|X_{\omega}(k)\right|^2}{N\sum_{n}\omega^2(n)}} \\ &= \sqrt{\frac{|X_{\omega}(k)|^2}{N\cdot S_2}} \end{align}\]

The corresponding linear power spectrum or RMS power spectrum can be defined by \[\begin{align} \text{PS}_L(k)&=X_L(k)^2=\frac{|X_{\omega}(k)|^2}{S_1^2}\frac{S_1^2}{N\cdot S_2} \\ &=\text{PS}(k) \cdot \frac{S_1^2}{N\cdot S_2} \end{align}\]

So, RMS can be calculated as below \[\begin{align} P_{\text{tot}} &= \sum \text{PS}_L(k) \\ \text{RMS} &= \sqrt{P_{\text{tot}}} \end{align}\]

DFT averaging

we use \(N= 8*\text{nfft}\) time samples of \(x\) and set the number of overlapping samples to \(\text{nfft}/2 = 512\). pwelch takes the DFT of \(\text{Navg} = 15\) overlapping segments of \(x\), each of length \(\text{nfft}\), then averages the \(|X(k)|^2\) of the DFT’s.

In general, if there are an integer number of segments that cover all samples of N, we have \[ N = (N_{\text{avg}}-1)*D + \text{nfft} \] where \(D=\text{nfft}-\text{noverlap}\). For our case, with \(D = \text{nfft}/2\) and \(N/\text{nfft} = 8\), we have \[ N_{\text{avg}}=2*N/\text{nfft}-1=15 \] For a given number of time samples N, using overlapping segments lets us increase \(N_{avg}\) compared with no overlapping. In this case, overlapping of 50% increases \(N_{avg}\) from 8 to 15. Here is the Matlab code to compute the spectrum:

1
2
3
4
5
6
7
8

nfft= 1024;
N= nfft*8;                                % number of samples in signal
n= 0:N-1;
x= A*sin(2*pi*f0*n*Ts) + .1*randn(1,N);   % 1 W sinewave + noise
noverlap= nfft/2;                         % number of overlapping time samples
window= hanning(nfft);
[pxx,f]= pwelch(x,window,noverlap,nfft,fs);  % W/Hz  PSD
PdB_bin= 10*log10(pxx*fs/nfft);              % dBW/bin

DFT averaging reduces the variance \(\sigma^2\) of the noise spectrum by a factor of \(N_{avg}\), as long as noverlap is not greater than nfft/2

fftshift

The result of fft and its index is shown as below

After fftshift

1
2
3
4
5

>> fftshift([0 1 2 3 4 5 6 7])

ans =

     4     5     6     7     0     1     2     3

1
2
3
4
5
6
7
8
9

clear;
N = 100;
fs = 1000;
fx = 100;
x = cos(2*pi*fx/fs*[0:1:N-1]);
s = abs(fftshift(fft(x)))/N;
fx = linspace(0,N-1,N)*fs/N -fs/2;
plot(fx, s, 'o-', 'linewidth', 2);
grid on; grid minor;

dft and psd function in virtuoso

• dft always return

  

To compensate windowing effect, \(W(n)\), the dft output should be multiplied by \(K_a\), e.g. 1.9922 for hanning window.

• psd function has taken into account \(K_e\), postprocessing is not needed

reference

Heinzel, Gerhard, A. O. Rüdiger and Roland Schilling. "Spectrum and spectral density estimation by the Discrete Fourier transform (DFT), including a comprehensive list of window functions and some new at-top windows." (2002). URL: https://holometer.fnal.gov/GH_FFT.pdf

Demonstration-DFT-PS-PSD

mathworks, Power Spectral Density Estimates Using FFT

Lyons, Richard, “Widowing Functions Improve FFT Results”, EDN, June 1, 1998. https://www.edn.com/electronics-news/4383713/Windowing-Functions-Improve-FFT-Results-Part-I

Rapuano, Sergio, and Harris, Fred J., An Introduction to FFT and Time Domain Windows, IEEE Instrumentation and Measurement Magazine, December, 2007. https://ieeexplore.ieee.org/document/4428580

Neil Robertson, The Power Spectrum

—, Use Matlab Function pwelch to Find Power Spectral Density – or Do It Yourself

—, Evaluate Window Functions for the Discrete Fourier Transform

—. Add a Power Marker to a Power Spectral Density (PSD) Plot [https://www.dsprelated.com/showarticle/1387.php]

Mathuranathan Viswanathan, Interpret FFT, complex DFT, frequency bins & FFTShift

robert bristow-johnson, Does windowing affect Parseval's theorem?

Window Correction Factors URL: https://community.sw.siemens.com/s/article/window-correction-factors

Brandt, A, 2011. Noise and vibration analysis: signal analysis and experimental procedures. John Wiley & Sons

What should be the correct scaling for PSD calculation using fft. URL:https://dsp.stackexchange.com/a/32205

POWER SPECTRAL DENSITY CALCULATION VIA MATLAB Revision C. URL: http://www.vibrationdata.com/tutorials_alt/psd_mat.pdf

Jens Ahrens, Carl Andersson, Patrik Höstmad, Wolfgang Kropp, “Tutorial on Scaling of the Discrete Fourier Transform and the Implied Physical Units of the Spectra of Time-Discrete Signals” in 148th Convention of the AES, e-Brief 56, May 2020 [ pdf, web ].

Manolakis, D., & Ingle, V. (2011). Applied Digital Signal Processing: Theory and Practice. Cambridge: Cambridge University Press. doi:10.1017/CBO9780511835261

• Laplace transform
  • a generalization of the continuous-time Fourier transform
  • converts integro-differential equations into algebraic equations
• z-transforms
  • a generalization of the discrete-time Fourier transform
  • converts difference equations into algebraic equations

system function, transfer function: \(H(s)\)

frequency response: \(H(j\omega)\), if the ROC of \(H(s)\) includes the imaginary axis, i.e.\(s=j\omega \in \text{ROC}\)

Laplace Transform

To specify the Laplace transform of a signal, both the algebraic expression and the ROC are required. The ROC is the range of values of \(s\) for the integral of \(t\) converges

bilateral Laplace transform

where \(s\) in the ROC and \(\mathfrak{Re}\{s\}=\sigma\)

The formal evaluation of the integral for a general \(X(s)\) requires the use of contour integration in the complex plane. However, for the class of rational transforms, the inverse Laplace transform can be determined without directly evaluating eq. (9.56) by using the technique of partial fraction expansion

ROC Property

The range of values of s for which the integral in converges is referred to as the region of convergence (which we abbreviate as ROC) of the Laplace transform

i.e. no pole in RHP for stable LTI sytem

System Causality

For a causal LTI system, the impulse response is zero for \(t \lt 0\) and thus is right sided

causality implies that the ROC is to the right of the rightmost pole, but the converse is not in general true, unless the system function is rational

System Stability

The system is stable, or equivalently, that \(h(t)\) is absolutely integrable and therefore has a Fourier transform, then the ROC must include the entire \(j\omega\)-axis

all of the poles have negative real parts

Unilateral Laplace transform

analyzing causal systems and, particularly, systems specified by linear constant-coefficient differential equations with nonzero initial conditions (i.e., systems that are not initially at rest)

A particularly important difference between the properties of the unilateral and bilateral transforms is the differentiation property \(\frac{d}{dt}x(t)\)

Integration by parts for unilateral Laplace transform

in Bilateral Laplace Transform

In fact, the initial- and final-value theorems are basically unilateral transform properties, as they apply only to signals \(x(t)\) that are identically \(0\) for \(t \lt 0\).

\(z\)-Transform

The \(z\)-transform for discrete-time signals is the counterpart of the Laplace transform for continuous-time signals

where \(z=re^{j\omega}\)

The \(z\)-transform evaluated on the unit circle corresponds to the Fourier transform

ROC Property

system stability

The system is stable, or equivalently, that \(h[n]\) is absolutely summable and therefore has a Fourier transform, then the ROC must include the unit circle

\(z\)-Transform Table

[http://www.ws.binghamton.edu/Fowler/fowler%20personal%20page/EE301_files/ZT%20Tables_rev2.pdf]

Unilateral \(z\)-Transform

The time shifting property is different in the unilateral case because the lower limit in the unilateral transform definition is fixed at zero, \(x[n-n_0]\)

bilateral \(z\)-transform

unilateral \(z\)-transform

Initial rest condition

Initial Value Theorem & Final Value Theorem

Laplace Transform

Two valuable Laplace transform theorem

• Initial Value Theorem, which states that it is always possible to determine the initial value of the time function \(f(t)\) from its Laplace transform \[ \lim _{s\to \infty}sF(s) = f(0^+) \]

• Final Value Theorem allows us to compute the constant steady-state value of a time function given its Laplace transform \[ \lim _{s\to 0}sF(s) = f(\infty) \]

  If \(f(t)\) is step response, then \(f(0^+) = H(\infty)\) and \(f(\infty) = H(0)\), where \(H(s)\) is transfer function

\(z\)-transform

Initial Value Theorem \[ f[0]=\lim_{z\to\infty}F(z) \] final value theorem \[ \lim_{n\to\infty}f[n]=\lim_{z\to1}(z-1)F(z) \]

Coert Vonk. Initial/final value proofs [https://coertvonk.com/physics/lfz-transforms/z/initial-final-value-proofs-31543]

\(s\)- and \(z\)-Domains Connection

discrete-time systems also can be analyzed by means of the Laplace transform — the \(z\)-transform is the Laplace transform in disguise and that discrete-time systems can be analyzed as if they were continuous-time systems.

A continuous-time system with transfer function \(H(s)\) that is identical in structure to the discrete-time system \(H[z]\) except that the delays in \(H[z]\) are replaced by elements that delay continuous-time signals

delay element	Time domain	Transform
\(z\)-transform	\(\delta[n-1]\)	\(z^{-1}\)
Laplace transform	\(\delta(t-T)\)	\(e^{-sT}\)

Therefore, we obtain

\[\begin{align} \sum_{n=0}^{\infty}y[n]z^{-n} &=\sum_{n=0}^{\infty}x[n]z^{-n}\cdot\sum_{n=0}^{\infty}h[n]z^{-n} \\ \sum_{n=0}^{\infty}y[n]e^{-snT} &=\sum_{n=0}^{\infty}x[n]e^{-snT}\cdot \sum_{n=0}^{\infty}h[n]e^{-snT} \end{align}\]

The \(z\)-transform of a sequence can be considered to be the Laplace transform of impulses sampled train with a change of variable \(z = e^{sT}\) or \(s = \frac{1}{T}\ln z\)

Note that the transformation \(z = e^{sT}\) transforms the imaginary axis in the \(s\) plane (\(s = j\omega\)) into a unit circle in the \(z\) plane (\(z = e^{sT} = e^{j\omega T}\), or \(|z| = 1\))

Note: \(\bar{h}(t)\) is the impulse sampled version of \(h(t)\)

Note \(\bar{h}(t)\) is impulse sampled signal, whose CTFT is scaled by \(\frac{1}{T}\) of continuous signal \(h(t)\), \(\overline{H}[e^{sT}]=\overline{H}(e^{sT})\) is the approximation of continuous time system response, for example summation \(\frac{1}{1-z^{-1}}\) \[ \frac{1}{1-z^{-1}} = \frac{1}{1-e^{-sT}} \approx \frac{1}{j\omega \cdot T} \] And we know transform of integral \(u(t)\) is \(\frac{1}{s}\), as expected there is ratio \(T\)

Staszewski, Robert Bogdan, and Poras T. Balsara. All-digital frequency synthesizer in deep-submicron CMOS. John Wiley & Sons, 2006

D. Sundararajan. Signals and Systems A Practical Approach Second Edition, 2023

Frequency Response of Discrete-Time LTI Systems

Eugenio Schuster, Lehigh University, Discrete-Time Signal Analysis in the Frequency-Domain [https://www.lehigh.edu/~eus204/teaching/ME450_SIRC/notes/DTS_Tutorial.pdf]

impulse invariance

A straightforward and useful relationship between the continuous-time impulse response \(h_c(t)\) and the discrete-time impulse response \(h[n]\)

\[ h[n] = Th_c(nT) \]

and \(T\) is chosen such that

\[ H_c(j\omega)=0, \space\space |\hat{\omega}| \ge \frac{\pi}{T} \]

When \(h[n]\) and \(h_c(t)\) are related through the above equation, i.e., the impulse response of the discrete-time system is a scaled, sampled version of \(h_c(t)\), the discrete-time system is said to be an impulse-invariant version of the continuous-time system

we have \[ H(e^{j\hat{\omega}}) = H_c\left(j\frac{\hat{\omega}}{T}\right),\space\space |\hat{\omega}| \lt \pi \]

• \(h[n] = Th_c(nT)\) and \(T\to 0\)

  guarantees \(y_c(nT) = y_r(nT)\), i.e. output equivalence only at the sampling instants

• \(H_c(j\Omega)\) is bandlimited and \(T \lt 1/2f_h\)

  guarantees \(y_c(t) = y_r(t)\)

The scaling of \(T\) can alternatively be thought of as a normalization of the time domain, that is average impulse response shall be same i.e., \(h[n]\times 1 = h(nT)\times T\)

note that the impulse-invariant transform is not appropriate for high-pass responses, because of aliasing errors

⭐ bandlimited

⭐ bandlimited

Scale the Result if impulse not multiplied with \(T_s\)

[https://www.mathworks.com/matlabcentral/answers/1582704-how-to-convolve-signal-and-transfer-function#comment_1825159]

[https://www.google.com.hk/search?q=matlab+continuous+transfer+function+in+time+domain+convolution]

1
2
3
4
5
6
7
8
9

t = linspace(0, 10);
dt = t(2) - t(1);
%u = @(t) heaviside(t); % by default, heavside(0) = 1/2, but at t = 0, we should have u(0) = 1 by definition of unit step
u = @(t) (t>=0);
Vin = @(t) u(t)-u(t-3);
h = @(t) ( exp(-2*t)/6 ) + ( 3/16*exp(-t) ) - ( exp(-5*t)/48 );
y = conv(Vin(t), h(t))*dt;  % must scale by dt
plot(t,y(1:numel(t)),'-x')
legend('wt - Laplace','wt - convolution integral','wt - convolution sum')

Transfer function & sampled impulse response

continuous-time filter designs to discrete-time designs through techniques such as impulse invariance

useful functions

• using fft

  The outputs of the DFT are samples of the DTFT

• using freqz

  modeling as FIR filter, and the impulse response sequence of an FIR filter is the same as the sequence of filter coefficients, we can express the frequency response in terms of either the filter coefficients or the impulse response

  fft is used in freqz internally

freqz method is straightforward, which apply impulse invariance criteria. Though fft is used for signal processing mostly,

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

clear all;
clc;

%% continuous system
s = tf('s');
h = 2*pi/(2*pi+s); 	% First order lowpass filter with 3-dB frequency 1Hz
[mag, phs, wout] = bode(h);
fct = wout(:)/2/pi;
Hct_dB = 20*log10(mag(:));


fstep = 0.01;           % freq resolution
fnyqst = 1000;
Ts = 1/(2*fnyqst);
Fs = 1/Ts;              % sampling freq
Ns = ceil(Fs/fstep);    % samping points
fstep = Fs/Ns;          % update fstep
t = (0:Ns-1)*Ts;        % sampling time points

y = impulse(h, t);      % impulse resp

%% modelling as discrete system
Y = fft(y);                 % dft
Hfft = Y * Ts;              % !!! multiply Ts
Hfft_dB = 20*log10(abs(Hfft(1:Ns/2+1)));
ffft = (1:Ns/2+1)*fstep - fstep;


[Hfir, ffir] = freqz(y, 1, [], 1/Ts);   % modelling as FIR
Hfir = Hfir * Ts;           % !!! multiply Ts
Hfir_dB = 20*log10(abs(Hfir));

%% plot
semilogx(fct, Hct_dB, '-ok', ffft, Hfft_dB, 'r.-', ffir, Hfir_dB, 'b--');
legend('bode(s)', 'fft', 'FIR model')
xlabel('Freq(Hz)');
ylabel('dB');
xlim([1e-2 1e2]);
grid on;
title('frequency response of different methods');

!!! only multiply Ts

Neil Robertson. The Discrete Fourier Transform as a Frequency Response [https://www.dsprelated.com/showarticle/1498.php]

H(k) by using complex exponentials is indeed identical to the DFT of h(n)

H(f) by DFT -> fft

H(k) by using complex exponential -> freqz

FIR Equalization

Frequency Response

\[ z = e^{j\omega T_s} \]

Unit impulse

filter coefficients are [-0.131, 0.595, -0.274] and sampling period is 100ps

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

%% Frequency response
w = [-0.131, 0.595, -0.274];
Ts = 100e-12;
[mag, w] = freqz(w, 1, [], 1/Ts);
subplot(2, 1, 1)
plot(w/1e9, 20*log10(abs(mag)));
xlabel('Freq(GHz)');
ylabel('dB');
grid on;
title('frequency response');

%% unit impulse response from transfer function
subplot(2, 1, 2)
z = tf('z', Ts);
h = -0.131 + 0.595*z^(-1) -0.274*z^(-2);
[y, t] = impulse(h);
stem(t*1e10, y*Ts);  % !!! y*Ts is essential
grid on;
title("unit impulse response");
xlabel('Time(\times 100ps)');
ylabel('mag');

impulse:

For discrete-time systems, the impulse response is the response to a unit area pulse of length Ts and height 1/Ts, where Ts is the sample time of the system. (This pulse approaches \(\delta(t)\) as Ts approaches zero.)

Scale output:

Multiply impulse output with sample period Ts in order to correct 1/Ts height of impulse function.

PSD transformation

If we have power spectrum or power spectrum density of both edge's absolute jitter (\(x(n)\)) , \(P_{\text{xx}}\)

Then 1UI jitter is \(x_{\text{1UI}}(n)=x(n)-x(n-1)\), and Period jitter is \(x_{\text{Period}}(n)=x(n)-x(n-2)\), which can be modeled as FIR filter, \(H(\omega) = 1-z^{-k}\), i.e. \(k=1\) for 1UI jitter and \(k=2\) Period jitter \[\begin{align} P_{\text{xx}}'(\omega) &= P_{\text{xx}}(\omega) \cdot \left| 1-z^{-k} \right|^2 \\ &= P_{\text{xx}}(\omega) \cdot \left| 1-(e^{j\omega T_s})^{-k} \right|^2 \\ &= P_{\text{xx}}(\omega) \cdot \left| 1-e^{-j\omega T_s k} \right|^2 \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

clear all
close all
clc

xf_fs = 0:0.01:0.5;
k = 1;
H_1UI = 1 - exp(-1i*2*pi*xf_fs);
HH_1UI = abs(H_1UI).^2;
subplot(2, 1, 1);
plot(xf_fs, HH_1UI);
grid on;
xlabel('Freq');
ylabel('|H|^2')
title('Weight for 1UI jitter');

k = 2;
H_period = 1 - exp(-1i*2*pi*xf_fs);
HH_period = abs(H_period).^2;
subplot(2, 1, 2)
plot(xf_fs, HH_period);
grid on;
xlabel('Freq');
ylabel('|H|^2')
title('Weight for Period jitter');

\[ x(t-\Delta T)\overset{FT}{\longrightarrow} X(s)e^{-\Delta T \cdot s} \]

Bilinear Transform

The design of IIR filters (cont.) [https://ocw.mit.edu/courses/2-161-signal-processing-continuous-and-discrete-fall-2008/cc00ac6d468dc9dcf2238fc1d1a194d4_lecture_19.pdf]

[https://tttapa.github.io/Pages/Mathematics/Systems-and-Control-Theory/Digital-filters/Discretization/Bilinear-transform.html]

an algebraic transformation between the variables \(s\) and \(z\) that maps the entire imaginary \(j\Omega\)-axis in the \(s\)-plane to one revolution of the unit circle in the \(z\)-plane

\[\begin{align} z &= \frac{1+s\frac{T_s}{2}}{1-s\frac{T_s}{2}} \\ s &= \frac{2}{T_s}\cdot \frac{1-z^{-1}}{1+z^{-1}} \end{align}\]

where \(T_s\) is the sampling period

frequency warping:

• The bilinear transformation avoids the problem of aliasing encountered with the use of impulse invariance, because it maps the entire imaginary axis of the \(s\)-plane onto the unit circle in the \(z\)-plane

• impulse invariance cannot be used to map highpass continuous-time designs to high pass discrete-time designs, since highpass continuous-time filters are not bandlimited

• Due to nonlinear warping of the frequency axis introduced by the bilinear transformation, bilinear transformation applied to a continuous-time differentiator will not result in a discrete-time differentiator. However, impulse invariance can be applied to bandlimited continuous-time differentiator

  The feature of the frequency response of discrete-time differentiators is that it is linear with frequency

The simple approximation \(z=e^{sT}\approx1+sT\), the first equal come from impulse invariance essentially

Stephen Roberts. Lecture 3 - Design of Digital Filters [https://www.robots.ox.ac.uk/~sjrob/Teaching/B14_SP/b14_sp_lect3.pdf]

Design of IIR Filters From Analog Filters

Modeling a Continuous-Time System with Matlab

Because the mapping between the continuous (\(s\)-plane) and discrete domains (\(z\)-plane) cannot be done exactly, the various design methods are at best approximations

Neil Robertson. Modeling a Continuous-Time System with Matlab [https://www.dsprelated.com/showarticle/1055.php]

—. Modeling Anti-Alias Filters [https://www.dsprelated.com/showarticle/1418.php]

impinvar [https://www.mathworks.com/help/signal/ref/impinvar.html]

bilinear [https://www.mathworks.com/help/signal/ref/bilinear.html]

filter [https://www.mathworks.com/help/matlab/ref/filter.html]

---

The design of IIR filters [https://ocw.mit.edu/courses/2-161-signal-processing-continuous-and-discrete-fall-2008/14f29be83ac38f11eb1a1bd5fdb6f5ea_lecture_18.pdf]

The design of IIR filters (cont.) [https://ocw.mit.edu/courses/2-161-signal-processing-continuous-and-discrete-fall-2008/cc00ac6d468dc9dcf2238fc1d1a194d4_lecture_19.pdf]

Lecture 19: Design of IIR Filters [http://smartdata.ece.ufl.edu/eee5502/2020_fall/media/2020_eee5502_slides28.pdf]

Derivatives Approximation

Perhaps the simplest method for low-order systems is to use backward-difference approximation to continuous domain derivatives

Note this approximation is not same with impulse invariance, e.g. \(\frac{1}{s^3} \to \frac{T^3z(z+1)}{2(-1)^3}\) employing impulse invariance

\[ 1- z^{-1} = 1-e^{-j\Omega T} = 1-\cos(\omega T) + j\sin(\Omega T) \approx 1-1+j\Omega T = s T \]

That is

\[ s \approx \frac{1-z^{-1}}{T} \]

Suppose we wish to make a discrete-time filter based on a prototype first-order low-pass filter \[ H_p(s) = \frac{1}{s\tau + 1} \] The differential equation describing this filter is \[ \tau\frac{dy}{dt} + y = x \] then differential equation gives

\[ \tau\frac{y_n-y_{n-1}}{T} + y_n=x_n \] or \[ (\frac{\tau}{T}+1)y_n - \frac{\tau}{T}y_{n-1} = x_n \] The transfer function is \[ H(z) = \frac{\frac{T}{T+\tau}}{1-\frac{\tau}{T+\tau}z^{-1}} \]

or substitute \(s\) with \(\frac{1-z^{-1}}{T}\) into \(H_p(s)\) \[ H_p(z) = \frac{1}{\tau\frac{1-z^{-1}}{T} + 1} = \frac{\frac{T}{T+\tau}}{1-\frac{\tau}{T+\tau}z^{-1}}= \frac{\alpha}{1 +(\alpha -1)z^{-1}} \] where \(\alpha = \frac{T}{\tau+T}\)

Under the assumption, time constants much longer than the timestep \(\tau \gg T\) \[\begin{align} \frac{T}{T+\tau}& = \frac{T}{\tau}\cdot \frac{\tau}{T+\tau}\approx \frac{T}{\tau} \\ \frac{\tau}{T+\tau} &= \frac{\tau - T}{\tau}\cdot\frac{\tau^2}{\tau^2-T^2} \approx \frac{\tau - T}{\tau} = 1- \frac{T}{\tau} \end{align}\]

Then the first-order low pass fiter \[ H_p(z) \approx \frac{\alpha}{1 +(\alpha -1)z^{-1}} \] where \(\alpha = \frac{T}{\tau}\)

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

# https://www.dsprelated.com/showarticle/1517/return-of-the-delta-sigma-modulators-part-1-modulation
# https://www.embeddedrelated.com/showarticle/779.php
def show_dsmod_samplewave(t,x,dsmod,args=(1,),tau=0.05, R=1, fig=None,
                          return_handles=False, filter_dsmod=False):
    dt = t[1]-t[0]
    if filter_dsmod:
        x1 = dsmod(x, *args,R=R, modulate=False)
    else:
        x1 = x
    y = dsmod(x,*args,R=R)
    a = dt/tau
    yfilt = scipy.signal.lfilter([a],[1,a-1],y)
    xfilt = scipy.signal.lfilter([a],[1,a-1],x1)

Matched z-Transform (Root Matching)

The matched Z-transform method, also called the pole-zero mapping or pole-zero matching method, and abbreviated MPZ or MZT

\[ z = e^{sT} \]

1
2
3
4
5
6
7
8
9
10

% https://www.dsprelated.com/showarticle/1642.php

%II One-pole RC filter model
Ts= 1/fs;
Wc= 1/(R*C);               % rad -3 dB frequency
fc= Wc/(2*pi);             % Hz -3 dB frequency
a1= -exp(-Wc*Ts);
b0= 1 + a1;                % numerator coefficient
a= [1 a1];                 % denominator coeffs
y_filt= filter(b0,a,y);    % filter the DAC's output signal y

Impulse Invariance (impinvar)

• To extend the accurate frequency range, we would need to increase the sample rate
• impulse-invariant transform is not appropriate for high-pass responses, because of aliasing errors

Note \(h[n] = Th_c(nT)\), Multiply the analog impulse response by this gain to enable meaningful comparison (other response, like step, the amplitude correction is not needed)

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

% https://www.dsprelated.com/showarticle/1055.php
% modified by hguo, Sun Jun 22 09:21:48 AM CST 2025

% butter_3rd_order.m      6/4/17 nr
% Starting with the butterworth transfer function in s,
% Create discrete-time filter using the impulse invariance xform and compare
% its time and frequency responses to those of the continuous time filter.
% Filter fc = 1 rad/s = 0.159 Hz


% I.  Given H(s), find H(z) using the impulse-invariant transform
fs= 4;            % Hz  sample frequency
% 3rd order butterworth polynomial
num= 1;
den= [1 2 2 1];
[b,a]= impinvar(num,den,fs);         % coeffs of H(z)
%[b,a]= bilinear(num,den,fs)

% II.  Impulse Response and Step Response
% find discrete-time impulse response
Ts= 1/fs;
N= 16*fs;
n= 0:N-1;
t= n*Ts;
x= [1, zeros(1,N-1)];   % impulse
x= fs*x;              % make impulse response amplitude independent of fs
y= filter(b,a,x);    % filter the impulse  
subplot(2,2,1),plot(t,y,'.'),grid,
xlabel('seconds')

% Continuous-time Impulse response from inverse Laplace transform
% Blinchikoff and Zverev, p116
h= exp(-t) - 2/sqrt(3)*exp(-t/2).*cos(sqrt(3)/2*t + pi/6); 
subplot(2,2,2),plot(t,y,'.',t,h),grid
title('Impulse Response'),legend('discrete-time filter response', 'continuous-time filter response'),xlabel('seconds')

% find discrete-time step response
x= ones(1,N);         % step
y= filter(b,a,x);    % filter the step         
% Continuous-time step response. Blinchikoff and Zverev, p116
t= t+Ts/2;                 % offset t to to align step responses
h= 1 - exp(-t) - 2/sqrt(3)*exp(-t/2).*sin(sqrt(3)/2*t); 
e= h-y;                    % error of discrete-time response
subplot(2,2,3),plot(t,y,'.',t,h),grid
title('Step Response'),legend('discrete-time filter response', 'continuous-time filter response'),xlabel('seconds')
% discrete-time response error
subplot(2,2,4),plot(t,e),grid
title('Error of discrete-time step response'),xlabel('seconds')

• To use the filter function with the b coefficients from an FIR filter, use y = filter(b,1,x).

Bilinear Transformation (bilinear)

The Filter Wizard Substack: "Countdown to s-to-z…" [https://filterwizard.substack.com/p/countdown-to-s-to-z]

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

a2 = 9.751405; a1 = 11.5677; a0 = 1.290113;
b2 = 0.118385; b1 = 30.27418; b0 = 0.502194;

fs = 1e3; % sampling rate 1KHz

num0 = (4*fs^2*a2 + 2*fs*a1 + a0)/(4*fs^2*b2 + 2*fs*b1 + b0);
num1 = (-8*fs^2*a2 + 2*a0)/(4*fs^2*b2 + 2*fs*b1 + b0);
num2 = (4*fs^2*a2 - 2*fs*a1 + a0)/(4*fs^2*b2 + 2*fs*b1 + b0);
den1 = (-8*fs^2*b2 + 2*b0)/(4*fs^2*b2 + 2*fs*b1 + b0);
den2 = (4*fs^2*b2 - 2*fs*b1 + b0)/(4*fs^2*b2 + 2*fs*b1 + b0);

s = tf('s');
H1_s = (a2*s^2 + a1*s + a0) / (b2*s^2 + b1*s + b0);

z = tf('z', 1/fs);
H1_z = (num0 + num1*z^(-1) + num2*z^(-2))/(1 + den1*z^(-1) + den2*z^(-2));

%% visualization
options = bodeoptions;
options.FreqUnits = 'Hz';

freq_range = logspace(log10(1e-3), log10(500), 1000);

[mag_s,phase_s,fout_s] = bode(H1_s, freq_range);
[mag_z,phase_z,fout_z] = bode(H1_z, freq_range);

subplot(2,1,1)
plot(fout_s(:), mag_s(:), 'r', LineWidth=4)
hold on
plot(fout_z(:), mag_z(:), 'b--', LineWidth=3)
title('Magnitude $H_1(s)$ \& $H_1(z)$ with $f_s=1$KHz',Interpreter='latex', FontSize=16)
legend('$|H_1(s)|$', '|$H_1(z)|$', fontsize=14, Location='southeast',Interpreter='latex')
xlabel('Hz', FontSize=12); ylabel('mag', FontSize=12); grid on;

subplot(2,1,2)
plot(fout_s(:), phase_s(:), 'r', LineWidth=4)
hold on
plot(fout_z(:), phase_z(:), 'b--', LineWidth=3)
title('Phase $H_1(s)$ \& $H_1(z)$ with $f_s=1$KHz',Interpreter='latex', FontSize=16)
legend('$\angle H_1(s)$', '$\angle H_1(z)$', fontsize=14, Location='southeast',Interpreter='latex')
xlabel('Hz', FontSize=12); ylabel('phase', FontSize=12); grid on;

disp(mag_s(1,1:10))
disp(mag_z(1,1:10))
%{
    2.5644    2.5643    2.5641    2.5640    2.5639    2.5637    2.5636    2.5635    2.5633    2.5632

    2.5644    2.5643    2.5641    2.5640    2.5639    2.5637    2.5636    2.5635    2.5633    2.5632
%}

bandwidth from pulse response

L.W. Couch, Digital and Analog Communication Systems, 8th Edition, Pearson, 2013.

Generating Basic signals – Rectangular Pulse and Power Spectral Density using FFT

Convolution Property of the Fourier Transform \[ x(t)*h(t)\longleftrightarrow X(\omega)H(\omega) \] pulse response can be obtained by convolve impulse response with UI length rectangular \[ H(\omega) = \frac{Y_{\text{pulse}}(\omega)}{X_{\text{rect}}(\omega)} = \frac{Y_{\text{pulse}}(\omega)}{\text{sinc}(\omega)} \]

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

% Convolution Property of the Fourier Transform
% pulse(t) = h(t) * rect(t)
% -> fourier transform
% PULSE = H * RECT
% FT(RECT) = sinc
% H = PULSE/RECT = PULSE/sinc
xx = pi*ui.*w(1:plt_num);
y_sinc = ui.*sin(xx)./xx;
y_sinc(1) = y_sinc(2);
y_sinc = y_sinc/y_sinc(1);  % we dont care the absoulte gain
h_ban1 = abs(h(1:plt_num))./abs(y_sinc);

h_dB = 20*log10(abs(h_ban1));
[hmin, Index] = min(abs(h_dB +3 ));
f_3dB = w(Index);
f_3dB = f_3dB/1e9;

sinc function

Notice that the complete definition of \(\operatorname{sinc}\) on \(\mathbb R\) is \[ \operatorname{Sa}(x)=\operatorname{sinc}(x) = \begin{cases} \frac{\sin x}{x} & x\ne 0, \\ 1, & x = 0, \end{cases} \] which is continuous.

To approach to real spectrum of continuous rectangular waveform, \(\text{NFFT}\) has to be big enough.

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

clear;
close all;
clc;

fs=500; %sampling frequency
Ts = 1/fs;
T=0.2; %width of the rectangule pulse in seconds

figure(1)
hold on;
t=-0.5:1/fs:0.5; %time base
x=rectpuls(t,T); %generating the square wave
sum(x>0.5)
stem(t,x,'--k');
plot(t, x, 'b.-')
xstart = T/2-Ts/2;
xend = -T/2-Ts/2;
plot([xstart, xstart], [-1, 1.02], 'r--');
plot([xend, xend], [-1, 1.02], 'r--');
plot([-0.11, 0.11], [1, 1], 'r--')
hold off;
grid on;
title(['Rectangular Pulse width=', num2str(T),'s']);
xlabel('Time(s)');
ylabel('Amplitude');
xlim([-0.12, 0.12]);
ylim([-0.05, 1.05]);

figure(2)
hold on
Titer = [1, 4, 16, 32, 64, 128];
color = {'k', 'g', 'r', 'm', 'c', 'b'};
formatSpec = 'NFFT=%d';
for k=1:length(Titer)
    t=-Titer(k):Ts:Titer(k); %time base
    x=rectpuls(t,T); %generating the square wave

    L=length(x);
    Np = nextpow2(L)-1;
    NFFT = 2^Np;
    X = fftshift(fft(x,NFFT)); %FFT with FFTshift for both negative & positive frequencies
    f = fs*(-NFFT/2:NFFT/2-1)/NFFT; %Frequency Vector
    Xc = abs(X)*Ts; % continuous signal spectrum
    plot(f, Xc, "Color",color{k}, 'LineWidth',2);
    legend_str{k} = [num2str(NFFT, formatSpec)];
end
hold off
legend(legend_str);
title('Magnitude of X_c (original continuous signal)');
xlabel('Frequency (Hz)')
ylabel('Magnitude |X(f)|');
grid on;

Stability check in z-domain

Unit circle criterion; Jury’s stability criterion

reference

Alan V. Oppenheim, Alan S. Willsky, and S. Hamid Nawab. 1996. Signals & systems (2nd ed.)

Alan V Oppenheim, Ronald W. Schafer. 2010. Discrete-Time Signal Processing, 3rd edition

B.P. Lathi, Roger Green. Linear Systems and Signals (The Oxford Series in Electrical and Computer Engineering) 3rd Edition

Sam Palermo, ECEN720, Lecture 7: Equalization Introduction & TX FIR Eq

Sam Palermo, ECEN720, Lab5 –Equalization Circuits

B. Razavi, "The z-Transform for Analog Designers [The Analog Mind]," IEEE Solid-State Circuits Magazine, Volume. 12, Issue. 3, pp. 8-14, Summer 2020. [https://www.seas.ucla.edu/brweb/papers/Journals/BR_SSCM_3_2020.pdf]

Jhwan Kim, CICC 2022, ES4-4: Transmitter Design for High-speed Serial Data Communications

Mathuranathan. Digital filter design – Introduction [https://www.gaussianwaves.com/2020/02/introduction-to-digital-filter-design/]

Daniel Boschen, "Fast Track to Designing FIR Filters with Python" [https://events.gnuradio.org/event/24/contributions/598/attachments/186/485/Boschen%20FIR%20Filter%20Presentation.pdf]

Daniel Boschen, "Quick Start on Control Loops with Python" [https://events.gnuradio.org/event/24/contributions/599/attachments/187/480/Boschen%20Control%20Presentation.pdf]

T-coil

coupling factor

note CTAP is grounded

\[ L_{12} = \frac{\text{im}[Z_{diff}]}{2\pi f} \]

where \(Z_{diff} = Z_{11} - Z_{12} - Z_{21} + Z_{22}\)

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

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt


df_emx  = pd.read_csv('./tcoil-emx.csv')
L1sim = df_emx['L1sim Y'].to_numpy()/1e-12  # pH
L2sim = df_emx['L2sim Y'].to_numpy()/1e-12  # pH
Q1sim = df_emx['Q1sim Y'].to_numpy()
Q2sim = df_emx['Q2sim Y'].to_numpy()
frqsim =  df_emx['ksim X'].to_numpy()/1e9   # GHz
Ksim =  df_emx['ksim Y'].to_numpy()

df_calc = pd.read_csv('./tcoil-calc.csv')
frqcalc = df_calc.iloc[:,0].to_numpy()/1e9   # GHz
L1calc = df_calc.iloc[:,1].to_numpy()/1e-12  # pH
L2calc = df_calc.iloc[:,3].to_numpy()/1e-12  # pH
Q1calc = df_calc.iloc[:,5].to_numpy()
Q2calc = df_calc.iloc[:,7].to_numpy()
L12calc = df_calc.iloc[:,9].to_numpy()/1e-12  # pH

kcalc = []
for l1i, l2i, l12i in zip(L1calc, L2calc, L12calc):
    l1il2i = l1i * l2i
    kk = (l12i - l1i -l2i)/2/np.sqrt(l1il2i) if l1il2i > 0.0 else 0.0   # L1 * L2 < 0 is meaningless
    kk = kk if np.abs(kk) < 2.0 else 0.0                                # constrain |k| < 2
    kcalc.append(-kk)

plt.figure(figsize=(20,8))
plt.subplot(3, 2, 1)
plt.plot(frqsim, L1sim, frqcalc, L1calc, 'r--', linewidth=3); plt.legend(['L1sim', 'L1calc'], loc='upper right'); plt.grid()
plt.subplot(3, 2, 2)
plt.plot(frqsim, Q1sim, frqcalc, Q1calc, 'r--', linewidth=3); plt.legend(['Q1sim', 'Q1calc'], loc='upper right'); plt.grid()
plt.subplot(3, 2, 3)
plt.plot(frqsim, L2sim, frqcalc, L2calc, 'r--', linewidth=3); plt.legend(['L2sim', 'L2calc'], loc='upper right'); plt.grid()
plt.subplot(3, 2, 4)
plt.plot(frqsim, Q2sim, frqcalc, Q2calc, 'r--', linewidth=3); plt.legend(['Q2sim', 'Q2calc'], loc='upper right'); plt.grid()
plt.subplot(3, 1, 3)
plt.plot(frqsim, Ksim, frqcalc, kcalc, 'r--', linewidth=3); plt.legend(['Ksim', 'kcalc'], loc='upper right'); plt.grid()

plt.show()

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

(get_k (lambda (l12 l1122) (letseq ((xvec (drGetWaveformXVec l12)) (n (drVectorLength xvec)) (l12v (drGetWaveformYVec l12)) (l1122v (drGetWaveformYVec l1122)) (resultv (drCreateVec 'double n))) (do ((i 0 i+1)) ((i >= n)) (letseq ((l12i (drGetElem l12v i)) (l1122i (drGetElem l1122v i)) (kk (if (l1122i > 0.0) l12i/(sqrt l1122i) 0.0)) (k (if ((abs kk) < 2.0) kk 0.0))) (drSetElem resultv i k))) (drCreateWaveform xvec resultv))))) (EMX_plot_aux bgui wid what 3 '("Inductance" "Q" "k") '("Henry" "" "") (lambda (ys) (letseq ((zs (reduce ys)) (z11 (nth 0 zs)) (z12 (nth 1 zs)) (z22 (nth 2 zs)) (pi 3.14159265358979) (f (xval z11)) (l11 (imag z11)/(2*pi*f)) (q11 (imag z11)/(real z11)) (l12 (imag z12)/(2*pi*f)) (l22 (imag z22)/(2*pi*f)) (q22 (imag z22)/(real z22)) (k (get_k l12 l11*l22))) `((,l11 ,l22) (,q11 ,q22) (,k)))) '(("L1" "L2") ("Q1" "Q2") ("k")))))

T-coil vs tapped inductor

tcoil and tapped inductor share same EM simulation result, and use modelgen with different model formula.

The relationship is \[ L_{\text{sim}} = L1_{\text{sim}}+L2_{\text{sim}}+2\times k_{\text{sim}} \times \sqrt{L1_{\text{sim}}\cdot L2_{\text{sim}}} \] where \(L1_{\text{sim}}\), \(L2_{\text{sim}}\) and \(k_{\text{sim}}\) come from tcoil model result, \(L_{\text{sim}}\) comes from tapped inductor model result

\(k_{\text{sim}}\) in EMX have assumption, induce current from P1 and P2 Given Dot Convention:

Same direction : k > 0

Opposite direction : k < 0

So, the \(k_{\text{sim}}\) is negative if routing coil in same direction

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

% EMX - shield tcoil model
L1 = csvread('./L1sim.csv', 1, 0);
L2 = csvread('./L2sim.csv', 1, 0);
k = csvread('./ksim.csv', 1, 0);

% EMX - Tapped shield inductor
L = csvread('./Lsim.csv', 1, 0);

freq = L1(:, 1)/1e9;    % GHz
L1 = L1(:, 2);
L2 = L2(:, 2);
k = -k(:, 2);   % Caution: minus of EMX ksim due to same current direction

L = L(:, 2);

Lcalc = L1 + L2 + 2*k.*(L1.*L2).^0.5;

plot(freq, L*1e9, 'r', 'LineWidth', 3);
hold on;
plot(freq, Lcalc*1e9, '--b', 'LineWidth', 3);
grid on;
xlabel('Freq (GHz)');
ylabel('Inductance (nH)');
legend('Tapped inductor model', 'tcoil model calc');

LVS check

LVS issue for circuits with customized devices

• auCdl: Analog and Microwave CDL, is a netlister used for creating CDL netlist for analog circuits

• auLVS: Analog and Microwave LVS, is used for analog circuit LVS

EMX ports

plain labels

• pin layer
• uncheck Cadence pins in Advanced options

rectangle pins

• drawing layer rectangle pin and specify Access Direction as intended
• check Cadence pins in Advanced options

The rectangle pins are always selected as driven port while there are only rectangle pin whether Cadence pins checked or not.

check ports used for simulation

use GDS view - EMX

EMX Synthesis Kits

Synthesis is a capability of the EMX Pcell library and uses scalable model data pre-generated by Continuum for a specific process and metal scheme combination.

Synthesis is supported by the Pcells that are suffixed _scalable, and these Pcells have the additional fields and buttons needed for synthesis.

port order (signals)

emxform.ils

EMX device info

EMX plot function

EMX's formulation is defined in

1	/path/to/EMX/share/emx/virtuoso_ui/emxinterface/emxskill/emxform.ils

EMX import this file at Virtuoso startup, you have to relaunch Virtuoso if you change this file

Single-ended inductor

Both with and without shield apply

• port-1 impedance when port-2 short

\[ Z_1 = \frac{1}{Y_{11}} \]

• port-2 impedance when port-1 short

\[ Z_2 = \frac{1}{Y_{22}} \] Then \[\begin{align} L1 &= \frac{Im(Z_1)}{2\pi f} \\ Q1 &= \frac{Im(Z_1)}{Re(Z_1)} \\ L2 &= \frac{Im(Z_2)}{2\pi f} \\ Q2 &= \frac{Im(Z_2)}{Re(Z_2)} \end{align}\]

EMX only plot L1 and Q1

differential impedance

Y parameters to Z parameters

\[\begin{align} |Y| &= Y_{11}*Y_{22} - Y_{12}*Y_{22} \\ \begin{bmatrix} Z_{11} & Z_{12}\\ Z_{21} & Z_{22} \end{bmatrix} &= \begin{bmatrix} \frac{Y_{22}}{|Y|} & \frac{-Y_{12}}{|Y|}\\ \frac{-Y_{21}}{|Y|} & \frac{Y_{11}}{|Y|} \end{bmatrix} \end{align}\]

Then differential impedance is \[ Z_{diff} = Z_{11} - Z_{12} - Z_{21} + Z_{22} \]

similarly, Z parameters to Y parameters \[ \begin{bmatrix} Y_{11} & Y_{12}\\ Y_{21} & Y_{22} \end{bmatrix} = \begin{bmatrix} \frac{Z_{22}}{|Z|} & \frac{-Z_{12}}{|Z|}\\ \frac{-Z_{21}}{|Z|} & \frac{Z_{11}}{|Z|} \end{bmatrix} \] where \[ |Z| = Z_{11}Z_{22} - Z_{12}Z_{21} \]

Differential inductor

Both with and without shield apply

\[\begin{align} L_{diff} &= \frac{Im(Z_{diff})}{2\pi f} \\ Q_{diff} &= \frac{Im(Z_{diff})}{Re(Z_{diff})} \end{align}\]

Center-tapped inductor

\[ Y = \begin{bmatrix} Y_{11} & Y_{12} & Y_{13}\\ Y_{21} & Y_{22} & Y_{23}\\ Y_{31} & Y_{32} & Y_{33} \end{bmatrix} \]

where port order is P1 P2 CT.

1
2
3
4
5
6
7
8
9
10
11

(define (EMX_plot_ct_ind bgui wid what)
  (EMX_plot_aux bgui wid what 3
		'("Differential inductance" "Differential Q")
		'("Henry" "")
		(lambda (ys)
		  (letseq ((z (EMX_differential (nth 0 ys) (nth 1 ys) (nth 3 ys) (nth 4 ys)))
			   (f (xval z))
			   (L (imag z)/(2*3.14159265358979*f))
			   (Q (imag z)/(real z)))
		    `((,L) (,Q))))
		'(("L") ("Q"))))

Assume CT i.e. port 3 in S-parameter is grounded, (z (EMX_differential (nth 0 ys) (nth 1 ys) (nth 3 ys) (nth 4 ys))) obtain differential impedance with \(Y_{11}\), \(Y_{12}\), \(Y_{21}\) and \(Y_{22}\). \[ Y = \begin{bmatrix} Y_{11} & Y_{12}\\ Y_{21} & Y_{22} \end{bmatrix} \] Finally, differential inductance and Q are obtained, shown as below

\[\begin{align} L_{diff} &= \frac{Im(Z_{diff})}{2\pi f} \\ Q_{diff} &= \frac{Im(Z_{diff})}{Re(Z_{diff})} \end{align}\]

Single-ended cap

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

(define (EMX_plot_se_cap bgui wid what)
  (EMX_plot_aux bgui wid what 2
		'("Capacitance" "Q" "Capacitance" "Q")
		'("Farad" "" "Farad" "")
		(lambda (ys)
		  (letseq ((z1 1.0/(nth 0 ys))
			   (y12 (nth 1 ys))
			   (z2 1.0/(nth 3 ys))
			   (f (xval z1))
			   (C1 (-1.0/(imag z1))/(2*3.14159265358979*f))
			   (C12 -(imag y12)/(2*3.14159265358979*f))
			   (C2 (-1.0/(imag z2))/(2*3.14159265358979*f))
			   (Q1 -(imag z1)/(real z1))
			   (Q12 (imag y12)/(real y12))
			   (Q2 -(imag z2)/(real z2)))
		    `((,C1) (,Q1) (,C12))))
		'(("Cse") ("Qse") ("C12"))))

We define Port-1 impedance \(Z_1\), Port-2 impedance \(Z_2\)

\[\begin{align} Z_1 &= \frac {1}{Y_{11}}\\ Z_2 &= \frac {1}{Y_{22}} \end{align}\]

Then single-ended cap and Q \[\begin{align} C_1 &= -\frac{1/Im(Z_1)}{2\pi f} \\ Q_1 &= -\frac{Im(Z_1)}{Re(Z_1)} \\ C_2 &= -\frac{1/Im(Z_2)}{2\pi f} \\ Q_2 &= -\frac{Im(Z_2)}{Re(Z_2)} \\ C_{12} &= -\frac{Im(Y_{12})}{2\pi f}\\ Q_{12} &= \frac{Im(Y_{12})}{Re(Y_{12})} \end{align}\]

• Series equivalent model is used in \(C_1\), \(Q_1\), \(C_2\) and \(Q_2\)
  • \(Z_1 = R + \frac{1}{sC_1}\) and \(Z_2 = R + \frac{1}{sC_2}\)
• Parallel model is used in \(C_{12}\) and \(Q_{12}\)
  • \(Y_{12} = \frac{1}{R} + sC_{12}\)

EMX plot \(C_{se}\), \(Q_{se}\) and \(C_{12}\), i.e. \(C_1\), \(Q_1\) and \(C_{12}\)

Differential cap

1
2
3
4
5
6
7
8
9
10
11

(define (EMX_plot_diff_cap bgui wid what)
  (EMX_plot_aux bgui wid what 2
		'("Differential capacitance" "Differential Q")
		'("Farad" "")
		(lambda (ys)
		  (letseq ((z (apply EMX_differential ys))
			   (f (xval z))
			   (C (-1.0/(imag z))/(2*3.14159265358979*f))
			   (Q -(imag z)/(real z)))
		    `((,C) (,Q))))
		'(("C") ("Q"))))

First obtain differential impedance, \(Z_{diff}\) then apply series equivalent model \[\begin{align} C_{diff} &= -\frac{1/Im(Z_{diff})}{2\pi f} \\ Q_{diff} &= -\frac{Im(Z_{diff})}{Re(Z_{diff})} \end{align}\]

Tline

Open circuit impedance \(Z_o\), short circuit impedance \(Z_s\) and characteristic impedance \(Z_0\)

\[\begin{align} Z_o &= Z_{11}\\ Z_s &= \frac{1}{Y_{11}}\\ Z_0 &= \sqrt{Z_o*Z_s} \end{align}\]

propagation constant is given as \[\begin{align} \gamma &= \frac{1}{2}\log\left( \frac{Z_0+Z_s}{Z_0-Z_s} \right) \\ &= \alpha + j\beta \end{align}\] where \(\alpha\) is attenuation constant and \(\beta\) is phase constant

The relationship between these parameter and geometry of the transmission line \[\begin{align} Z_0 &= \sqrt{\frac{R+j\omega L}{G+j\omega C}} \\ \gamma &= \sqrt{(G+j\omega C)(R+j\omega L)} \end{align}\] EMX plot the real and imaginary part of \(Z_0\), \(\alpha\) and \(\beta\) of \(\gamma\)

Note EMX plot the absolute value of \(\alpha\) and \(\beta\)

Transformer

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

  (get_k
   (lambda (l12 l1122)
     (letseq ((xvec (drGetWaveformXVec l12))
       (n (drVectorLength xvec))
       (l12v (drGetWaveformYVec l12))
       (l1122v (drGetWaveformYVec l1122))
       (resultv (drCreateVec 'double n)))
(do ((i 0 i+1))
    ((i >= n))
  (letseq ((l12i (drGetElem l12v i))
	   (l1122i (drGetElem l1122v i))
	   (kk (if (l1122i > 0.0)
		   l12i/(sqrt l1122i)
		   0.0))
	   (k (if ((abs kk) < 2.0) kk 0.0)))
    (drSetElem resultv i k)))
(drCreateWaveform xvec resultv)))))
(EMX_plot_aux bgui wid what num_ports
'("Inductance" "Q" "k")
'("Henry" "" "")
(lambda (ys)
  (letseq ((zs (reduce ys))
	   (z11 (nth 0 zs))
	   (z12 (nth 1 zs))
	   (z22 (nth 2 zs))
	   (pi 3.14159265358979)
	   (f (xval z11))
	   (l11 (imag z11)/(2*pi*f))
	   (q11 (imag z11)/(real z11))
	   (l12 (imag z12)/(2*pi*f))
	   (l22 (imag z22)/(2*pi*f))
	   (q22 (imag z22)/(real z22))
	   (k (get_k l12 l11*l22)))
    `((,l11 ,l22) (,q11 ,q22) (,k))))
'(("L1" "L2") ("Q1" "Q2") ("k")))))

[IC Prophet GDSII 文件使用和仿真测试说明]

EMX autoplot

using AC simulation, and inductor's parallel model or series model

That is to say: both sp (network parameter) and ac (impedance) can be used to plot inductance, Q value.

usually EMX choose ac method

left 2 figures are used for AC simulation, \(Y_{nn}\) can be obtained conveniently

Foundary model

for single-end capicator \[\begin{align} Q_1 &= -\frac{Im(Z_1)}{Re(Z_1)} \\ &= -\frac{Im(1/Y_{11})}{Re(1/Y_{11})} \\ &= -\frac{Im(Y_{11}^*)/|Y_{11}|^2}{Re(Y_{11}^*)/|Y_{11}|^2} \\ &= \frac{Im(Y_{11})}{Re(Y_{11})} \end{align}\]

So, the EMX model and foundary model is consistent.

pad & bump

EMX process file contain M0 up to RDL-AP

PM, CB2_FC, UBM is in the chip package

PEX extract up to RDL-AP as expected

DC model

quick find routing resistance

GPDK045 metal resistor model is not consistent with its process file on sheet resistance

Tips

Process file

Process file encryption mostly for advanced nodes, like TSMC 16nm Finfet, whose process file is encrypted.

• Use --key=EMXkey in the EMX Advanced options

GDSviewer has two options

• EMX: shows the final gds sent to EMX for simulation after it has been processed by EMX
• Raw: shows the raw gds

If there are port name with the # sign, it means EMX sees a port but it is not in the signal list.

EMX Accuracy

• Edge mesh: controls layout discretization in the X-Y plane

  • For MoM capacitors, use the edge mesh to be the same as the width of the finger (for example, 0.1um).

• Thickness: controls layout discretization in the Z dimension

• 3D metals: skips all 2D assumptions about conductors and their currents and charges

  • If you set 3D metals to * then all metals are treated as 3D
    • For Inductor type structures, only thick metal needs 3D.
    • For MoM, all layers are needed.

Ports entered in Grounds will cause these nets to be grounded; these ports will not show up in the S-parameter result.

Setup Temperature

• EMX: --temperature=100

ParaView

• If check ParaView related options when ParaView is not setup properly, EMX simulation stop at Creating mesh... without waring or errors (version 6.2).

Paraview & stimulus

reference

Tips on Specifying Ports in EMX

Using 'Cadence pins' as ports with access direction in EMX simulations

EMX miscellaneous features [https://picture.iczhiku.com/resource/eetop/WyIFKleSLTRIuvCb.pdf]

Cadence Rapid Adoption Kit (RAK). Analysis of a Figure-Eight Inductor with EMX

discrete time jitter impulse response (normalized to the input jitter stimulus similar to the procedure used to represent a conventional system impulse response)

When impulsive jitter is injected into clock distribution circuits (i.e., a small incremental time delay or advance applied to an individual clock edge), it results in jitter in multiple subsequent edges in the output clock

transient noise and rms_jitter function

RJ(rms): single Edge or Both Edge?

RJ(seed): what is it?

phase noise method

Directly compare the input phase noise and output phase noise, the input waveform maybe is the PLL output or other clock distribution end point

Jitter Impulse Response & Jitter Transfer Function

assuming linear, time-invariant phase response

n =5 buffers, fclk = 10GHz

[Alphawave’s CTO, Tony Chan Carusone, High Speed Communications Part 8 – On Die CMOS Clock Distribution]

Theoretically, the DC gain of JTF of low pass filter shall be 1. Unfortunately, the gain less than 1 due to numerical error or nonlinearity

Example

Low Pass Filter

1
2
3
4
5
6
7
8
9
10
11

N = 32;
x = zeros(N,1);
x(1) = 6;
x(2) = -2;
x(3) = 0.5;
x = x/5;
figure(1)
stem(x)
Y = fft(x, N);
figure(2)
plot(abs(Y(1:N/2)));

discrete time jitter impulse response

both input and output are discrete time signal, i.e. no sampling in the input, that's why ratio \(1/T_s\) is not in the jtf

High Pass Filter

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

N = 128;
j_hp = zeros(N, 1);
j_hp(1)= 1;
j_hp(2) = 0.5;
j_hp(3) = -0.3;
j_hp(4) = 0.3;
j_hp(5) = -0.1;
jtf_hp = abs(fft(j_hp));
semilogx(jtf_hp(1:N/2+1));
xlabel('Freq');
ylabel('Jitter Amplification Factor');
grid on;

inverter chain

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

ji = 1e-12; % 1ps
data = importdata('/path/to/jir.csv');
jo = data.data(:, 2);
Ts = 31.25e-12;
Fs = 1/Ts;
jir = jo/ji;
N = 2^(nextpow2(length(jir)-1));
Y = fft(jir, N);
jtf = abs(Y(1:N/2+1));
freqs= Fs/N*(0:N/2);
plot(feqs/1e9, jtf, 'linewidth', 2);
grid on;
xlabel('Freq (GHz)');
ylabel('JTF');
title('JTF of inverter chain');

Colored Jitter Amplification

[Alphawave’s CTO, Tony Chan Carusone, High Speed Communications Part 8 – On Die CMOS Clock Distribution]

Four major noise sources are included in the modeling: Input noise, DAC quantization noise (DAC QN), DCO random noise (DCO RN), and delay line random noise (DL RN).

H. Kang et al., "A 42.7Gb/s Optical Receiver with Digital CDR in 28nm CMOS," 2023 IEEE Radio Frequency Integrated Circuits Symposium (RFIC), San Diego, CA, USA, 2023 [https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=10630516]

Dirac impulse at edge position

Full-rate JTF

Singe edge is using

Half-rate

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

ji = 10;
jir = readmatrix("~/cargo/jir.csv");
jir = jir / ji;
jir_2e = jir(:, 2);		  % both edge
jir_1e = jir(1:2:end, 2); % single edge

[mag, w] = freqz(jir_2e, 1, [], 1);
plot(w, abs(mag), LineWidth=2);
hold on
[mag, w] = freqz(jir_1e, 1, [], 0.5);
plot(w, abs(mag), LineWidth=2);
grid on;
legend("2Edge", "1Edge")
title("JTF")

fck = 1GHz

Note: Phase Noise dBc is SSB, that's why we add 10log10(2)

The calculated JTF from JIR is too small compared with Pnoise simualtion

245.1/272.4 = 89.98%

245.1/272.4 = 83.74%

Reference

Sam Palermo, ECEN 720, Lecture 13 - Forwarded Clock Deskew Circuits

B. Casper and F. O'Mahony, "Clocking Analysis, Implementation and Measurement Techniques for High-Speed Data Links-A Tutorial," in IEEE Transactions on Circuits and Systems I. [https://people.engr.tamu.edu/spalermo/ecen689/clocking_analysis_hs_links_casper_tcas1_2009.pdf]

Rhee, W. (2020). Phase-locked frequency generation and clocking : architectures and circuits for modern wireless and wireline systems. The Institution of Engineering and TechnologyMathuranathan Viswanathan, Digital Modulations using Matlab : Build Simulation Models from Scratch

Tony Chan Carusone, University of Toronto, Canada, 2022 CICC Educational Sessions "Architectural Considerations in 100+ Gbps Wireline Transceivers"

X. Mo, J. Wu, N. Wary and T. Chan Carusone, "Design Methodologies for Low-Jitter CMOS Clock Distribution," in IEEE Open Journal of the Solid-State Circuits Society, 2021 [https://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=9559395]

Y. Zhao and B. Razavi, "Phase Noise Integration Limits for Jitter Calculation," 2022 IEEE International Symposium on Circuits and Systems (ISCAS), Austin, TX, USA, 2022 [https://www.seas.ucla.edu/brweb/papers/Conferences/YZ_ISCAS_22.pdf]

Thomas Toifl. TWEPP 2012. Low-power High-Speed CMOS I/Os: Design Challenges and Solutions [https://indico.cern.ch/event/170595/contributions/266344/attachments/212179/297391/twepp_sept2012_final_v2.pdf]

Ganesh Balamurugan and Naresh Shanbhag, "Modeling and mitigation of jitter in multiGbps source-synchronous I/O links," Proceedings 21st International Conference on Computer Design, San Jose, CA, USA, 2003, pp. 254-260, doi: 10.1109/ICCD.2003 [https://shanbhag.ece.illinois.edu/publications/ganesh-ICCD2203.pdf]

Balamurugan, G. & Casper, Bryan & Jaussi, James & Mansuri, Mozhgan & O'Mahony, Frank & Kennedy, Joseph. (2009). Modeling and Analysis of High-Speed I/O Links. Advanced Packaging, IEEE Transactions on. [https://sci-hub.se/10.1109/TADVP.2008.2011366]

Jihwan Kim, ISSCC2019 F5: Design Techniques for a 112Gbs PAM-4 Transmitter

Process Variation

Global/Local Variation

Timing and RC Modeling with Process Corners

Global and Local variation by Gaussian

Local Monte-Carlo (SSG, FFG with Local Gaussian) as Signoff golden

Process Corner Model Limitations

Variation section

• Total corner (TT/SS/FF/SF/FS)
  • E.g. TTMacro_MOS_MOS_MOSCAP
• Global Corner (TTG/SSG/FFG/SFG/FSG) + Local MC
  • E.g. TTGlobalCorner_LocalMC_MOS_MOSCAP
• Local MC
  • E.g. LocalMCOnly_MOS_MOSCAP
• Global MC + Local MC (Total MC)
  • GlobalMC_LocalMC_MOS_MOSCAP

SSGNP, FFGNP:

When N/P global correlation is weak (R^2=0.15), the corner of N/PMOS balance circuit (e.g. inverter) can be tightened (3sigma -> 2.5sgma) due to the cancellation between NMOS and PMOS

SSGNP, FFGNP usually used in Digital STA

• Global variation validation with global corner
  • 3-sigma of global MC simulation is aligned with global corner
• Total variation validation with total corner
  • 3-sigma of global MC + local MC (total) simulation is aligned with total corner

Global Corner

https://community.cadence.com/cadence_technology_forums/f/custom-ic-design/20466/monte-carlo-simulation-global-local-vs-local-and-process-vs-mismatch/1365101#1365101

The "total corner" is representative of the maximum device parameter variation including local device variation effects. However, it is not a statistical corner.

The "global" corner is defined as the "total" corner minus the impact of "local variation"

Hence, if you were to examine simulation results for a parameter using a "total" and "global" corner, you would find the range of variation will be less with the "global" corner than with the "total" corner.

The "global" corner is provided for use in statistical simulations. Hence, when performing a Monte-Carlo simulation, the "global" corner is selected - NOT the "total" corner.

\[ \Delta V_{T,\sigma_{total}} = \sqrt{\Delta^2 _{T, \sigma_{global}}+\Delta^2 _{T, \sigma_{local}}} \]

RC corner

Traditional RC Corners

Metal width variation (\(\Delta W\)), Metal thickness variation (\(\Delta T\)), IMD thickness variation (\(\Delta H\))

Capacitance Dominant: C-best, C-worst

Resistance Dominant: RC-best, RC-worst

DPT effect

When using two masks per layer (Double Patterning Technology, DPT) there is an issue of mask alignment where any mis-alignment will cause layer spacing values to change, therefore changing the parasitic coupling capacitance values.

Misalignment scale and direction are not deterministic facts: coupling cap and total cap may be increased or decreased.

Five new corners are added in a DPT flow to account for RC variations accurately:

sapced-dependent side-wall dielectric constant also affect coupling cap

and CC_worst means to increase both K1 and K2

CC_best means to decrease both K1 and K2

• Setup time sign-off would use:

  Cworst_CCworst / RCworst_CCworst

• Hold time sign-off would use:

  Cbest_CCbest / RCbest_CCbest / Cworst_CCworst / RCworst_CCworst

Signoff corner

BEOL Target: typical

The recommended RC corner:

​ cworst_CCworst, cbest_CCbest, rcworst_CCworst, rcbest_CCbest and typical

The others are for pre-color RC calculation purpose

_T stands for "Tighten DPT corner"; these are less pessimistic 1.5 sigma corners

Below table is caputre of Aragio's TSMC16: LVDS datasheet

BEOL corner

Spacing variation is implicitly defined by \(\Delta W_m\).

We denote the conductor width and thickness of the layer m by \(W_m\) and \(T_m\), respectively.

Similarly, we denote the thickness of the layer's interlayer dielectric (i.e., the distance between layer m and layer m +1) by \(H_m\)

• C-based means worst and best caps
• RC-based means worst and best R in adjustment with C (RC product)

Based on experience, it was found that C-based extraction provides worst and best case over RC for internal timing paths because Capacitance dominates short wire.

However, for large design, inter-block timing paths were often worst with RC worst parasitic since R dominates for long wires.

signoff corners for setup & hold

reference

Process Variation

Eric J.-W. Fang, T5: Fundamentals of Process Monitors for Signoff-Oriented Circuit Design, 2022 IEEE International Solid-State Circuits Conference

Alvin Loke, Device and Physical Design Considerations for Circuits in FinFET Technology, ISSCC 2020 Short Course

簡報 Cln16ffcll Sr V1d0 2p1 Usage Guide URL: https://usermanual.wiki/Document/cln16ffcllsrv1d02p1usageguide.1649731847/view

Radojcic, Riko, Dan Perry and Mark Nakamoto. “Design for manufacturability for fabless manufactuers.” IEEE Solid-State Circuits Magazine 1 (2009): n. pag.

How To Reduce Implementation Headaches In FinFET Processes URL: https://semiengineering.com/how-to-reduce-implementation-headaches-in-finfet-processes/

陌上风骑驴看IC, STA | SSGNP, FFGNP. https://mp.weixin.qq.com/s/eJ8fYRJBR1E9XbfH95OUOg

陌上风骑驴看IC, STA | ssg 跟ss corner 的区别——谬误更正版 https://mp.weixin.qq.com/s?__biz=MzUzODczODg2NQ==&mid=2247486225&idx=1&sn=e9c68f6108ae6c9958d47ca0b29373ca&chksm=fad262cfcda5ebd949cc91353c7cbfaf4ba61179306f7d8e98461a4f4ca9d8a9baef5e9f2cc1&scene=178&cur_album_id=1326356275000705025#rd

The Evolution, Pitfalls, and Cargo Cult Engineering of Advanced Digital Timing Sign-off https://www.tauworkshop.com/2021/speaker_slides/christian_l.pdf

Don O'Riordan Cadence Design Systems. Recommended Spectre Monte Carlo modeling methodology [https://designers-guide.org/modeling/montecarlo.pdf]

RC corner

Modeling Sub-90nm On-Chip Variation Using Monte Carlo Method for DFM https://www.aspdac.com/aspdac2007/pdf/archive/2D-1.pdf

Double Patterning for IC Design, Extraction and Signoff https://semiwiki.com/eda/synopsys/1974-double-patterning-for-ic-design-extraction-and-signoff/

抽刀断水水更流，RC Corner不再愁：STA之RC Corner URL: http://mp.weixin.qq.com/s?__biz=MzUzODczODg2NQ==&mid=2247484115&idx=1&sn=de99f27aadf58ea316c284dad9000b7c&chksm=fad26b0dcda5e21b8c9750f738b55053f695843a66c3c202ff0ba586c738f45aa270254c3722&scene=21#wechat_redirect

一曲新词酒一杯，RC Corner继续飞: STA之RC Corner拾遗 URL:https://mp.weixin.qq.com/s?__biz=MzUzODczODg2NQ==&mid=2247484135&idx=2&sn=bddc632850bd10c32b5688fd7af46218&chksm=fad26b39cda5e22f1c3970f8c8c2e1287c9492c526c4caf02b61f61faffdf829381c392d6ea1&scene=21#wechat_redirect

且将新火试新茶，深究趁年华：STA之RC Corner再论 URL:https://mp.weixin.qq.com/s?__biz=MzUzODczODg2NQ==&mid=2247484144&idx=1&sn=059843381e77cd4008d25166db388d02&chksm=fad26b2ecda5e23816b33b3a949f34d4ca09118bad76f1089dfcb228b7da6491f423a5f4e703&cur_album_id=1326356275000705025&scene=189#wechat_redirect

LDP_IN_800_25V_DN: 1.2GHz LVDS Receiver http://aragiosolutions.com/pdf/rgo_tsmc16_18v25_lvds_product_brief_rev_1a.pdf

Parasitic extraction technologies Advanced node and 3D-IC design https://static.sw.cdn.siemens.com/siemens-disw-assets/public/81845/en-US/Siemens-SW-Parasitic-extraction-technologies-for-advanced-WP-81845-C2.pdf

New Game, New Goal Posts: A Recent History of Timing Closure https://pdfs.semanticscholar.org/9360/5ce48f9bd3b7527ae8979f41a9c7e310efa4.pdf

The Evolution, Pitfalls, and Cargo Cult Engineering of Advanced Digital Timing Sign-off https://www.tauworkshop.com/2021/speaker_slides/christian_l.pdf

T. -B. Chan, S. Dobre and A. B. Kahng, "Improved signoff methodology with tightened BEOL corners," 2014 IEEE 32nd International Conference on Computer Design (ICCD), Seoul, Korea (South), 2014, pp. 311-316, doi: 10.1109/ICCD.2014.6974699.

Chan, T. (2014). Mitigation of Variability and Reliability Margins in IC Implementation /. UC San Diego. ProQuest ID: Chan_ucsd_0033D_14269. Merritt ID: ark:/20775/bb52916761. Retrieved from https://escholarship.org/uc/item/35r1m001

Dr. Adam Temanm, Digital VLSI Design:Lecture 10: Routing https://www.eng.biu.ac.il/temanad/files/2017/02/Lecture-10-Routing.pdf

Article (20487193) Title: Setting Pegasus - LVS to Quantus av_extracted view Flow with TSMC16 packages URL: https://support.cadence.com/apex/ArticleAttachmentPortal?id=a1O0V000009MprZUAS

primary clock, generated clock and virtual clock in SDC

primary clocks

• Primary clocks should be created at input ports and output pins of black boxes.
• Never create clocks on hierarchy pins. Creating clocks on hierarchy will cause problems when reading SDF. The net timing arc becomes segmented at the hierarchy and PrimeTime will be unable to annotate the net successfully.

generated clocks

• Generated clocks are generally created for waveform modifications of a primary clock (not including simple inversions). PrimeTime does not simulate a design and thus will not derive internally generated clocks automatically - these clocks must be created by the user and applied as a constraint.
• PrimeTime caculate source latency for generated clocks if primary clock is propagated, otherwise its source latency is zero.

1
2
3
4
5

# primary clock
create_clock -period 4 [get_ports Clk]
set_clock_latency -source 2 [get_clocks Clk]
set_propagated_clock [get_clocks Clk]
create_generated_clock -divide_by 2 -name div_clk -source [get_ports Clk] FF3/Q

virtual clocks

• Are clock objects without a source
• Do not clock sequential devices within the current_design
• Serve as references of input or output delays

1
2
3
4

# create a virtual clock, vclk, for input and output delay constraints
create_clock -period 5 -name vclk
set_input_delay -max 2 -clock vclk [get_ports in1]
set_output_delay -max 1 -clock vclk [get_ports out2]

There is no network latency to calculate, even if a virtual clock is propagated.

Clock Edges used for Setup and Hold in PrimeTime

PT picks the most restrictive pair of edges for setup and for hold. It determines which edges to be used as follows:

1. Evaluate waveforms over the smallest common base period

2. For each capture edge, find the closest setup launch edge. Call these the primary pairs

3. Out of the primary pairs, pick the most restrictive setup launch and capture edges.

4. For each primary pair, draw two hold relationships:

   • Launch to (capture - 1)

   • (Launch + 1) to Capture

     From all of these hold relationships, pick the most restrictive.

PrimeTime uses the ideal clock waveform (as reported in report_clock) to determine the appropriate clock edges for inter-clock analysis.

The most restrictive setup pair is from Clk1 8ns to Clk2 9ns

The most restrictive hold pair is from Clk1 0ns to Clk2 0ns

Multicycle Paths Constraints in PrimeTime

PrimeTime® User Guide Version O-2018.06-SP4 Chapter 1: Introduction to PrimeTime Overview of Static Timing Analysis - Timing Exceptions

PrimeTime does not automatically identifies multicycle paths

specifying multicycle path for setup

1
2
3
4
5

set_multicycle_path 6 -from reg[26]/CP -to reg/D
# or
set_multicycle_path -setup 6 -from reg[26]/CP -to reg/D
# check the exception
report_exception

specifying multicycle path for hold (new data every 6 cycles)

1
2

set_multicycle_path -setup 6 -to [get_pins "*reg[*]/D"]
set_multicycle_path -hold [expr 6-1] -to [get_pins "*reg[*]/D"]

M_H stands for Hold Multiplier, M_S for Setup Multiplier. The Setup multiplier counts up with increasing clock cycles, the Hold Multiplier counts up with decreasing cycles. The origin (0) for the Hold Multiplier is always at the Setup Multiplier - 1 position.

Reporting a multicycle path with report_timing

1

report_timing -exceptions all -from *reg[26]/CP -to *reg/D

Timing Exceptions

If certain paths are not intended to operate according to the default setup and hold behavior assumed by the PrimeTime tool, you need to specify those paths as timing exceptions. Otherwise, the tool might incorrectly report those paths as having timing violations.

The PrimeTime tool lets you specify the following types of exceptions:

• False path – A path that is never sensitized due to the logic configuration, expected data sequence, or operating mode.
• Multicycle path – A path designed to take more than one clock cycle from launch to capture.
• Minimum or maximum delay path – A path that must meet a delay constraint that you explicitly specify as a time value.

package import in SystemVerilog compilation

class forward declaration [https://verificationacademy.com/forums/systemverilog/class-forward-declaration#answer-41383]

import and compile order

project

yapp_pkg.sv

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

package yapp_pkg;
import uvm_pkg::*;
`include "uvm_macros.svh"

 typedef uvm_config_db#(virtual yapp_if) yapp_vif_config;

`include "yapp_packet.sv"
`include "yapp_tx_monitor.sv"
`include "yapp_tx_sequencer.sv"
`include "yapp_tx_seqs.sv"
`include "yapp_tx_driver.sv"
`include "yapp_tx_agent.sv"
`include "yapp_env.sv"

endpackage

yapp_tx_monitor.sv

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

class yapp_tx_monitor extends uvm_monitor;
  
  // Collected Data handle
  yapp_packet pkt;

  // Count packets collected
  int num_pkt_col;

  // component macro
  `uvm_component_utils_begin(yapp_tx_monitor)
    `uvm_field_int(num_pkt_col, UVM_ALL_ON + UVM_NOCOMPARE)
  `uvm_component_utils_end
/////////////////// virtual interface ////////////////////
  virtual interface yapp_if vif;
/////////////////////////////////////////////////////////
      
  function new (string name, uvm_component parent);
    super.new(name, parent);
  endfunction : new
      ...
     
endclass : yapp_tx_monitor

yapp_tx_driver.sv

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

class yapp_tx_driver extends uvm_driver #(yapp_packet);

 ////////////////// virtual interface ////////////////////
  virtual interface yapp_if vif;
 /////////////////////////////////////////////////////
      
  int num_sent;

  // component macro
  `uvm_component_utils_begin(yapp_tx_driver)
    `uvm_field_int(num_sent, UVM_ALL_ON + UVM_NOCOMPARE)
  `uvm_component_utils_end

  // Constructor - required syntax for UVM automation and utilities
  function new (string name, uvm_component parent);
    super.new(name, parent);
  endfunction : new

  function void start_of_simulation_phase(uvm_phase phase);
    `uvm_info(get_type_name(), {"start of simulation for ", get_full_name()}, UVM_HIGH)
  endfunction : start_of_simulation_phase
      ...
     
endclass : yapp_tx_driver 

yapp_if.sv

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

interface yapp_if (input clock, input reset );
timeunit 1ns;
timeprecision 100ps;

import uvm_pkg::*;
`include "uvm_macros.svh"

////////////////// import package //////////////////////
import yapp_pkg::*;
//////////////////////////////////////////////////////    

  // Actual Signals
  logic              in_data_vld;
  logic              in_suspend;
  logic       [7:0]  in_data;

  // signal for transaction recording
  bit monstart, drvstart;

  // local storage for payload
  logic [7:0] payload_mem [0:63];
    ...
    
endinterface : yapp_if

virtual interface is in package;

yapp_pkg is imported into interface

!!! typedef uvm_config_db#(virtual yapp_if) yapp_vif_config; is forward declaration

A forward typedef declares an identifier as a type in advance of the full definition of that type

VCS compile

1

vcs -full64 -R -sverilog  -ntb_opts uvm-1.2 +UVM_TESTNAME=short_yapp_012_test  -f vcs.f  

vcs.f

1
2
3
4
5
6
7
8
9

-timescale=1ns/1ns

+incdir+../sv

../sv/yapp_if.sv
../sv/yapp_pkg.sv
./clkgen.sv
./hw_top.sv
./tb_top.sv

output ERROR

1
2
3
4
5
6

Error-[SV-LCM-PND] Package not defined
../sv/yapp_if.sv, 18
yapp_if, "yapp_pkg::"
  Package scope resolution failed. Token 'yapp_pkg' is not a package. 
  Originating module 'yapp_if'.
  Move package definition before the use of the package.

xrun compile

1

xrun -f xrun.f

xrun.f

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

-64

-uvmhome  CDNS-1.2

// incdir for include files
-incdir ../sv

+UVM_TESTNAME=short_yapp_012_test

-timescale 1ns/1ns

// compile files
../sv/yapp_if.sv
../sv/yapp_pkg.sv
./clkgen.sv
./hw_top.sv
./tb_top.sv

output Error

1
2
3
4
5
6

file: ../sv/yapp_if.sv
import yapp_pkg::*;
              |
xmvlog: *E,NOPBIND (../sv/yapp_if.sv,18|14): Package yapp_pkg could not be bound.
        interface worklib.yapp_if:sv
                errors: 1, warnings: 0

solution

place ../sv/yapp_pkg.sv before ../sv/yapp_if.sv.

In this particular example, yapp_pkg is NOT needed in interface. just delete import yapp_pkg::* is enough in ../sv/yapp_if.sv

plain example

The order of compilation unit DON'T matter

project

hw_top.sv

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

module hw_top;

  // Clock and reset signals
  logic [31:0]  clock_period;
  logic         run_clock;
  logic         clock;
  logic         reset;

  // YAPP Interface to the DUT
  yapp_if in0(clock, reset);

  // CLKGEN module generates clock
  clkgen clkgen (
    .clock(clock),
    .run_clock(1'b1),
    .clock_period(32'd10)
  );

  yapp_router dut(
    .reset(reset),
    .clock(clock),
    .error(),
    // YAPP interface signals connection
    .in_data(in0.in_data),
    .in_data_vld(in0.in_data_vld),
    .in_suspend(in0.in_suspend),
    // Output Channels
    //Channel 0   
    .data_0(),
    .data_vld_0(),
    .suspend_0(1'b0),
    //Channel 1   
    .data_1(),
    .data_vld_1(),
    .suspend_1(1'b0),
    //Channel 2   
    .data_2(),  
    .data_vld_2(),
    .suspend_2(1'b0),
    // Host Interface Signals
    .haddr(),
    .hdata(),
    .hen(),
    .hwr_rd());
    
    ...
        
endmodule

yapp_router.sv

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

module yapp_router (input clock,                              
                    input reset,                            
                    output error,

                    // Input channel
                    input [7:0] in_data,                           
                    input in_data_vld,                     
                    output in_suspend,
                    // Output Channels
                    output [7:0] data_0,  //Channel 0
                    output reg data_vld_0, 
                    input suspend_0, 
                    output [7:0] data_1,  //Channel 1
                    output reg data_vld_1, 
                    input suspend_1, 
                    output [7:0] data_2,  //Channel 2
                    output reg data_vld_2,
                    input suspend_2,
     
                    // Host Interface Signals
                    input [15:0] haddr,
                    inout [7:0] hdata,
                    input hen,
                    input hwr_rd);                            
...
endmodule    

compile

place hw_top.sv before or after yapp_router.sv doesn't matter, the compiler (xrun, vcs) can compile them successfully.

Conclusion

Always place package before DUT is preferred choice during compiling.

Transaction Level Modeling (TLM)

• Blocking methods are defined with get() or put() tasks to allow them to consume time
• Non-blocking methods are defined with try_get() or try_put() functions as they execute in zero time

Uni-Directional TLM Methods Reference

The peek() methods are similarly to the get() methods, but copy the transaction instead of removing it. The transaction is not consumed, and a subsequent get or peek operation will return the same transaction

Selected Connector and Method Options

in the connectors above, * can be replaced by port, imp, or export

All the methods for a specific connector type MUST be implemented. If you define an uvm_put connection between two compoents, then the component with the uvm_put_imp object must provide implementations of ALL three put methods, put, try_put and can_put, even if these methods are not explicitly called

TLM FIFO

The TLM FIFO is a FIFO component wrapped in get and put imp connectors. This has the benefit of data storage as well as providing implementations of the communication methods. Components connected to the TLM FIFO are in control of data transfer and can simply defined port connectors to initiate read and write operations on the FIFO

uvm_tlm_fifo

The TLM FIFO object is effectively a FIFO component instantiated between and connected to two components. The FIFO contains imp connectors for the standard TLM put and get/peek interfaces, therefore the user does not have to defineimp ports or communication methods and the FIRO takes care of data storage

The advantages are:

• The user does not need to define communication methods or imp connectors
• The FIFO provides data storage between the write (put) and read (get/peek) components
• There are a number of built-in methods for checking FIFO status

The disadvantages are:

• The user must now initiate both sides of the transfer (both get/peek and put) to complete the transaction
• Two connections must be made (both sides of the FIFO) rather than one

The put_export and get_peek_export connection objects have alternatives which provide subsets of the full connector. For example, blocking_put_export and nonblocking_put_export can replace put_export. blocking_get_export, nonblocking_get_export and get_export (as well as others) can replace get_peek_export.

built-in methods

1
2
3

class uvm_tlm_fifo #(type T=int) extends uvm_tlm_fifo_base #(T);
    ...
endclass    

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

virtual class uvm_tlm_fifo_base #(type T=int) extends uvm_component;
    uvm_put_imp #(T, this_type) put_export;
    
    uvm_get_peek_imp #(T, this_type) get_peek_export;
    
    uvm_analysis_port #(T) put_ap;
    
    uvm_analysis_port #(T) get_ap;
    
    // The following are aliases to the above put_export
    uvm_put_imp      #(T, this_type) blocking_put_export;
    uvm_put_imp      #(T, this_type) nonblocking_put_export;
    
    // The following are all aliased to the above get_peek_export, which provides
    // the superset of these interfaces.
    uvm_get_peek_imp #(T, this_type) blocking_get_export;
    uvm_get_peek_imp #(T, this_type) nonblocking_get_export;
    uvm_get_peek_imp #(T, this_type) get_export;
    
    uvm_get_peek_imp #(T, this_type) blocking_peek_export;
    uvm_get_peek_imp #(T, this_type) nonblocking_peek_export;
    uvm_get_peek_imp #(T, this_type) peek_export;
    
    uvm_get_peek_imp #(T, this_type) blocking_get_peek_export;
    uvm_get_peek_imp #(T, this_type) nonblocking_get_peek_export;
    
    function new(string name, uvm_component parent = null);
        super.new(name, parent);
        
        put_export = new("put_export", this);
        blocking_put_export     = put_export;
        nonblocking_put_export  = put_export;
        
        get_peek_export = new("get_peek_export", this);
        blocking_get_peek_export    = get_peek_export;
        nonblocking_get_peek_export = get_peek_export;
        blocking_get_export         = get_peek_export;
        nonblocking_get_export      = get_peek_export;
        get_export                  = get_peek_export;
        blocking_peek_export        = get_peek_export;
        nonblocking_peek_export     = get_peek_export;
        peek_export                 = get_peek_export;
        
        put_ap = new("put_ap", this);
        get_ap = new("get_ap", this);
        
    endfunction

Analysis FIFO

uvm_tlm_analysis_fifo

uvm_tlm_analysis_fifo is a specialization of uvm_tlm_fifo

• Intended to buffer write transactions between the UVC monitor analysis port and scoreboard

It has the following characteristics:

• Unbounded (size=0)
• analysis_export connector replaces put_export
  • Support the analysis write method

By declaring the FIFO in the scoreboard, we can get directly from the FIFO output.

However the write side connection of the FIFO to the interface UVC monitor analysis port must be made (usually) in the testbench which has visibility of both UVC and scoreboard components. The connection is made using a connect method call inside the connect phase method

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