KNOW-HOW

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;

reference

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

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

Deterministic Jitter

j_Djpp can be calculated by PSD,too

1
2
3
4
5
6
7
8
9

fck = 38.4e6;
Nfft = 15000;
fres = fck/Nfft;
psddBc = -99.3343;
psBc = psddBc + 10*log10(fres);	% psd -> ps; 
phrad2 = 10^(psBc/10);
phrms = sqrt(phrad2);
Jrms = phrms/2/pi*1/fck;
Jpp = 2*sqrt(2)*Jrms;

1
2
3

Jpp =

   6.4038e-12

For DJ, we usually use peak to peak value

BTW, the psd value at half of fundamental frequency (\(f_s/2\)) is duty cycle distortion due to the NMOS/PMOS imbalance, because of rising only data

Random Jitter

RJ can be accurately and efficiently measured using PSS/Pnoise or HB/HBnoise.

Note that the transient noise can also be used to compute RJ;

However, the computation cost is typically very high, and the accuracy is lesser as compared to PSS/Pnoise and HB/HBnoise.

Since RJ follows a Gaussian distribution, it can be fully characterized using its Root-Mean-Squared value (RMS) or the standard deviation value (\(\sigma\))

The Peak-to-Peak value of RJ (\(\text{RJ}_{\text{p-p}}\)) can be calculated under certain observation conditions \[ \text{RJ}_{\text{p-p}}\equiv K \ast \text{RJ}_{\text{RMS}} \] Here, \(K\) is a constant determined by the BER specification of the system given in the following Table

1
2
3
4

K = 14.698;
Ks = K/2;
p = normcdf([-Ks Ks]);
BER = 1 - (p(2)-p(1));

1
2
3

BER =

   1.9962e-13

Total Jitter

\[ \text{TJ}_{\text{p-p}}\equiv \text{DJ}_{\text{p-p}} + \text{RJ}_{\text{p-p}}(\text{BER}) \]

In the psd of TJ, the spur is DJ and floor is RJ

Phase Noise to Jitter

The phase noise is traditionally defined as the ratio of the power of the signal in 1Hz bandwidth at offset \(f\) from the carrier \(P\), divided by the power of the carrier \[ \ell (f) = \frac {S_v'(f_0+f)}{P} \] where \(S_v'\) is is one-sided voltage PSD and \(f \geqslant 0\)

Under narrow angle assumption \[ S_{\varphi}(f)= \frac {S_v'(f_0+f)}{P} \] where \(\forall f\in \left[-\infty +\infty\right]\)

Using the Wiener-Khinchin theorem, it is possible to easily derive the variance of the absolute jitter(\(J_{ee}\))via integration of the corresponding PSD \[ J_{ee,rms}^2 = \int S_{J_{ee}}(f)df \]

And we know the relationship between absolute jitter and excess phase is \[ J_{ee}=\frac {\varphi}{\omega_0} \] Considering that phase noise is normally symmetrical about the zero frequency, multiplied by two is shown as below \[ J_{ee,rms} = \frac{\sqrt{2\int_{0}^{+\infty}\ell(f)df}}{\omega_0} \] where phase noise is in linear units not in logarithmic ones.

Because the unit of phase noise in Spectre-RF is logarithmic unit (dBc), we have to convert the unit before applying the above equation \[ \ell[linear] = 10^{\frac {\ell [dBc/Hz]}{10}} \] The complete equation using the simulation result of Spectre-RF Pnoise is \[ J_{ee,rms} = \frac{\sqrt{2\int_{0}^{+\infty}10^{\frac {\ell [dBc/Hz]}{10}}df}}{\omega_0} \]

The above equation has been verified for sampled pnoise, i.e. J_ee and Edge Phase Noise.

• For pnoise-sampled(jitter), Direct Plot Form - Function: Jee:Integration Limits can calculate it conveniently
• But for pnoise-timeaveage, you have to use the below equation to get RMS jitter.

One example, integrate to \(\frac{f_{osc}}{2}\) and \(f_{osc} = 16GHz\)

Of course, it apply to conventional pnoise simulation.

On the other hand, output rms voltage noise, \(V_{out,rms}\) divied by slope should be close to \(J_{ee,rms}\) \[ J_{ee,rms} = \frac {V_{out,rms}}{slope} \]

Pulse Width Jitter (PWJ)

[Spectre Tech Tips: Measuring Noise in Digital Circuits]

Pnoise sampled: Edge Delay mode measures the noise defined by two edges. Both edges are defined by a threshold voltage and rising or falling edges, which measures the noise of the pulse itself and direct plot calculate the variation of the pulse width

reference

Article (20500632) Title: How to simulate Random and Deterministic Jitters URL: https://support.cadence.com/apex/ArticleAttachmentPortal?id=a1O3w000009fiXeEAI

Spectre Tech Tips: Measuring Noise in Digital Circuits - Analog/Custom Design - Cadence Blogs - Cadence Community https://community.cadence.com/cadence_blogs_8/b/cic/posts/s . . .

J. Kim et al., "A 112 Gb/s PAM-4 56 Gb/s NRZ Reconfigurable Transmitter With Three-Tap FFE in 10-nm FinFET," in IEEE Journal of Solid-State Circuits, vol. 54, no. 1, pp. 29-42, Jan. 2019, doi: 10.1109/JSSC.2018.2874040

— et al., "A 224-Gb/s DAC-Based PAM-4 Quarter-Rate Transmitter With 8-Tap FFE in 10-nm FinFET," in IEEE Journal of Solid-State Circuits, vol. 57, no. 1, pp. 6-20, Jan. 2022, doi: 10.1109/JSSC.2021.3108969

DC offset

Performing FFT to a signal with a large DC offset would often result in a big impulse around frequency 0 Hz, thus masking out the signals of interests with relatively small amplitude.

One method to remove DC offset from the original signal before performing FFT

• Subtracting the Mean of Original Signal

You can also not filter the input, but set zero to the zero frequency point for FFT result.

Nyquist component

If we go back to the definition of the DFT \[ X(N/2)=\sum_{n=0}^{N-1}x[n]e^{-j2\pi (N/2)n/2}=\sum_{n=0}^{N-1}x[n]e^{-j\pi n}=\sum_{n=0}^{N-1}x[n](-1)^n \] which is a real number.

The discrete function \[ x[n]=\cos(\pi n) \] is always \((-1)^n\) for integer \(n\)

One general sinusoid at Nyquist and has phase shift \(\theta\), this is \(T=2\) and \(T_s=1\)

\[\begin{align} x[n] &= A \cos(\pi n + \theta) \\ &= A \big( \cos(\pi n) \cos(\theta) - \sin(\pi n) \sin(\theta) \big) \\ &= \big(A\cos(\theta)\big) \cos(\pi n) + \big(-A\sin(\theta)\big) \sin(\pi n) \\ &= \big(A\cos(\theta)\big) (-1)^n + \big(-A\sin(\theta)\big) \cdot 0 \\ &= B \cdot (-1)^n \\ \end{align}\]

Where \(A\cos(\theta)=B\).

Moreover \(B \cdot (-1)^n = B\cdot \cos(\pi n)\), then \[ B\cdot \cos(\pi n) = A \cdot \cos(\pi n + \theta) \] We can NOT distinguish one from another.

In other words, you CAN'T infer the signal from \(X(\frac{N}{2})\) \[\begin{align} X(k)\frac{1}{N}e^{j 2 \pi \frac{nk}{N}}\bigg|_{k=\frac{N}{2}} &= \frac{X\left(\frac{N}{2} \right)}{N}(-1)^n \\ &= \frac{X\left(\frac{N}{2} \right)}{N}\cos(\pi n) \\ &= \frac{X\left(\frac{N}{2} \right)}{N}\left( \cos(\pi n) - \beta \sin(\pi n) \right) \\ &= \frac{X\left(\frac{N}{2} \right)}{N}\sqrt{1+\beta^2}\left(\frac{1}{\sqrt{1+\beta^2}} \cos(\pi n) - \frac{\beta}{\sqrt{1+\beta^2}} \sin(\pi n) \right) \\ &= \frac{X\left(\frac{N}{2} \right)}{N} \frac{1}{\cos(\theta)}\left(\cos(\theta) \cos(\pi n) - \sin(\theta) \sin(\pi n) \right) \\ &= \frac{X\left(\frac{N}{2} \right)}{N} \frac{1}{\cos(\theta)} \cos(\pi n+\theta) \end{align}\]

where \(\beta \in \mathbb{R}\) and you wouldn't know it because \(\sin(\pi n)=0 \quad \forall n \in \mathbb{Z}\)

For example, if \(\theta=0\) \[ X(k)\frac{1}{N}e^{j 2 \pi \frac{nk}{N}}\bigg|_{k=\frac{N}{2}}= \frac{X\left(\frac{N}{2} \right)}{N} \cos(\pi n) \] However, if \(\theta=\frac{\pi}{3}\) \[ X(k)\frac{1}{N}e^{j 2 \pi \frac{nk}{N}}\bigg|_{k=\frac{N}{2}}= \frac{X\left(\frac{N}{2} \right)}{N}\cdot 2 \cos(\pi n+\frac{\pi}{3}) \]

That sort of ambiguity is the reason for the strict inequality of the sampling theorem's condition.

Duty Cycle Distortion

[Prof. Tony Chan Carusone, Low-Jitter CMOS Clock Distribution]

Both edges is used for clock waveform to evaluate the duty cycle distortion,

Assuming TIE is [0 0.2 0 0.2 0 0.2 ...], then subtract DC offset, we get [-0.1 0.1 -0.1 0.1 ...], shown as below

The amplitude is manifested in FFT amplitude spectrum, i.e. Nyquist component, which is 0.1 in follow figure

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

N = 32;
n = (1:N);
x = 0.1*(-1).^n;
figure(1)
stem(n-1, x);
X = fft(x)/N;
Xshift = fftshift(X);
fprintf("nyquist component: %.2f\n", Xshift(1));
magXshift = abs(Xshift);
ph = phase(Xshift)/pi*180;
figure(2)
subplot(3, 1, 1)
fx = (-N/2:N/2-1);
stem(fx, magXshift);
xlabel('Freq');
ylabel('|X(k)|');
title('mag of DFT');
grid on;

subplot(3, 1, 2)
stem(fx, ph);
xlabel('Freq');
ylabel('\angle X(k)(^oC)');
title('phase of DFT');
grid on;

%% inverse dft
ninv = (0:32-1);
xinv = Xshift(1)*cos(pi*ninv);
subplot(3, 1, 3);
hold on;
stem(ninv, x, "filled", 'r');
stem(ninv, xinv,'bd-.');
ninfer = (0:0.1:32+1);
xinfer1 = Xshift(1)*cos(pi*ninfer); % theta = 0
xinfer2 = Xshift(1)*2*cos(pi*ninfer+pi/3); % theta = pi/3
plot(ninfer, xinfer1, 'm--');
plot(ninfer, xinfer2, 'c--');
hold off;
legend('original', 'IDFT', '\theta=0', '\theta=\pi/3');
xlabel('time');
ylabel('V');
title('sample points');
grid on;

Single-Sided Amplitude Spectrum

DC and Nyquist frequency of FFT left over

P1(2:end-1) = 2*P1(2:end-1);

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

Fs = 1000;            % Sampling frequency
T = 1/Fs;             % Sampling period
L = 1500;             % Length of signal
t = (0:L-1)*T;        % Time vector
S = 0.7*sin(2*pi*50*t) + sin(2*pi*120*t);
X = S + 2*randn(size(t));
figure(1)
plot(1000*t(1:50),X(1:50))
title('Signal Corrupted with Zero-Mean Random Noise')
xlabel('t (milliseconds)')
ylabel('X(t)')

figure(2)
Y = fft(X);
P2 = abs(Y/L);      %!!! two-sided spectrum P2.
P1 = P2(1:L/2+1);   %!!! single-sided spectrum P1
P1(2:end-1) = 2*P1(2:end-1);	% exclude DC and Nyquist freqency
f = Fs*(0:(L/2))/L;
figure(2)
plot(f,P1)
title('Single-Sided Amplitude Spectrum of X(t)')
xlabel('f (Hz)')
ylabel('|P1(f)|')

Alternative View

The direct current (DC) bin (\(k=0\)) and the bin at \(k=N/2\), i.e., the bin that corresponds to the Nyquist frequency are purely real and unique.

sinusoidal waveform with \(10Hz\), amplitude 1 is \(cos(2\pi f_c t)\). The plot is shown as below with sampling frequency is \(20Hz\)

Amplitude and Phase spectrum, sampled with \(f_s=20\) Hz

The FFT magnitude of \(10Hz\) is 1 and its phase is 0 as shown as above, which proves the DFT and IDFT.

Caution: the power of FFT is related to samples (DFT Parseval's theorem), which may not be the power of continuous signal. The average power of samples is ([1 -1 1 -1 -1 1 ...]) is 1, that of corresponding continuous signal is \(\frac{1}{2}\).

Power spectrum derived from FFT provide information of samples, i.e. 1

Moreover, average power of sample [1 -1 1 -1 1 ...] is same with DC [1 1 1 1 ...].

code

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

fc = 10;
fs = 2*fc;
fov = 64*fc;

ts = (0:1/fs:26);
tov = (0:1/fov:26);

ys = cos(2*pi*fc*ts);
yov = cos(2*pi*fc*tov);
%% waveform
stem(ts, ys)
hold on;
plot(tov, yov);
legend('sample', 'waveform')
ylim([-1.5 1.5])
grid on;
xlabel('time(s)');
ylabel('mag(V)')

nfft = 256;
X = fftshift(fft(ys, nfft))/nfft;
f = (-nfft/2:nfft/2-1)*fs/nfft;
magX = abs(X);
phsX = atan2(imag(X),real(X));
%% fft spectrum
figure(2)
subplot(2, 1, 1);
stem(f, magX);
xlabel('Frequency(Hz)');
ylabel('mag')
xlim([min(f)-1 max(f)+1])
title('Amplitude spectrum')
subplot(2, 1, 2)
plot(f, phsX);
xlabel('Frequency(Hz)');
ylabel('phase (rad)')
xlim([min(f)-1 max(f)+1])
title('Phase spectrum')

%% power spectrum
yssq_sum_avg = sum(ys(1:nfft).^2)/nfft;
specsq_sum_avg = sum(abs(X).^2);

reference

OriginLab, How to Remove DC Offset before Performing FFT URL: http://blog.originlab.com/how-to-remove-dc-offset-before-performing-fft

How to remove DC component in FFT? URL: https://www.mathworks.com/matlabcentral/answers/712808-how-to-remove-dc-component-in-fft#answer_594373

Analyzing a signal that contains frequency content at Fs/2 doesn't seem to work unless there is a phase shift URL: https://dsp.stackexchange.com/a/59807/59253

Nyquist–Shannon sampling theorem, Critical frequency URL: https://en.wikipedia.org/wiki/Nyquist%E2%80%93Shannon_sampling_theorem#Critical_frequency

Why remove energy at Nyquist before ifft? URL: https://dsp.stackexchange.com/a/22851/59253

1

emacs --no-site-file --load path/to/verilog-mode.el --batch filename.v -f verilog-auto-save-compile

CAUTION: filename.v is overwrite by command

verilog-mode.el

1
2
3
4

/*AUTOINPUT*/
/*AUTOWIRE*/
/*AUTOINST*/
/*AUTO_TEMPLATE*/

-f verilog-batch-auto

For use with --batch, perform automatic expansions as a stand-alone tool. This sets up the appropriate Verilog mode environment, updates automatics with M-x verilog-auto on all command-line files, and saves the buffers. For proper results, multiple filenames need to be passed on the command line in bottom-up order.

-f verilog-auto-save-compile

Update automatics with M-x verilog-auto, save the buffer, and compile

Emacs

--no-site-file

Another file for site-customization is site-start.el. Emacs loads this before the user's init file (.emacs, .emacs.el or .emacs.d/.emacs.d). You can inhibit the loading of this file with the option --no-site-file

--batch

The command-line option --batch causes Emacs to run noninteractively. The idea is that you specify Lisp programs to run; when they are finished, Emacs should exit.

--load, -l FILE, load Emacs Lisp FILE using the load function;

--funcall, -f FUNC, call Emacs Lisp function FUNC with no arguments

-f FUNC

--funcall, -f FUNC, call Emacs Lisp function FUNC with no arguments

--load, -l FILE

--load, -l FILE, load Emacs Lisp FILE using the load function

Verilog-mode is a standard part of GNU Emacs as of 22.2.

multiple directories

AUTOINST only search in the file's directory default.

You can append below verilog-library-directories for multiple directories search

1
2
3

// Local Variables:
// verilog-library-directories:("." "subdir" "subdir2")
// End:

reference

Emacs Online Documentation https://doc.endlessparentheses.com/

Emacs verilog-mode 的使用 URL: https://www.wenhui.space/docs/02-emacs/verilog_mode_useguide/

always@( * )

always@( * ) blocks are used to describe Combinational Logic, or Logic Gates. Only = (blocking) assignments should be used in an always@( * ) block.

Latch Inference

If you DON'T assign every element that can be assigned inside an always@( * ) block every time that always@( * ) block is executed, a latch will be inferred for that element

The approaches to avoid latch generation:

• set default values
• proper use of the else statement, and other flow constructs

without default values

latch is generated

RTL

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

module TOP (
	input wire Trigger,
	input wire Pass,
	output reg A,
	output reg C
);
	always @(*) begin
		A = 1'b0;
		if (Trigger) begin
			A = Pass;
			C = Pass;
		end
	end
endmodule

synthesized netlist

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

/////////////////////////////////////////////////////////////
// Created by: Synopsys DC Ultra(TM) in wire load mode
// Version   : S-2021.06-SP5
// Date      : Mon May  9 17:09:18 2022
/////////////////////////////////////////////////////////////


module TOP ( Trigger, Pass, A, C );
  input Trigger, Pass;
  output A, C;


  lanhq1 C_reg ( .E(Trigger), .D(Pass), .Q(C) );
  an02d0 U3 ( .A1(Pass), .A2(Trigger), .Z(A) );
endmodule

add default value

Default values are an easy way to avoid latch generation

RTL

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

module TOP (
	input wire Trigger,
	input wire Pass,
	output reg A,
	output reg C
);
	always @(*) begin
		A = 1'b0;
        	C = 1'b1;
		if (Trigger) begin
			A = Pass;
			C = Pass;
		end
	end

synthesized netlist

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

/////////////////////////////////////////////////////////////
// Created by: Synopsys DC Ultra(TM) in wire load mode
// Version   : S-2021.06-SP5
// Date      : Mon May  9 17:12:47 2022
/////////////////////////////////////////////////////////////


module TOP ( Trigger, Pass, A, C );
  input Trigger, Pass;
  output A, C;


  nd12d0 U5 ( .A1(Pass), .A2(Trigger), .ZN(C) );
  an02d0 U6 ( .A1(Pass), .A2(Trigger), .Z(A) );
endmodule

if evaluation

signed number cast to unsigned automatically before evaluating

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

// tb.v
module tb; 
  reg signed [1:0] datasg;
  reg [1:0] dataug;

  initial begin
    datasg = 2'b11;
    dataug = 2'b11;

    $display("datasg(%%d): %d", datasg);
    $display("dataug(%%d): %d", dataug);

    if (datasg)
      $display("datasg is OK");
    if (dataug)
      $display("dataug is OK");

    $finish();
  end 
endmodule

1
2
3
4
5
6
7
8
9

$ vlog tb.v
$ vsim -c -do "run;exit" work.tb
# Loading work.tb(fast)
# run
# datasg(%d): -1
# dataug(%d): 3
# datasg is OK
# dataug is OK
# ** Note: $finish    : tb.v(16)

reference

UC Berkeley CS150 Lec #20: Finite State Machines Slides

Lee WF. Learning from VLSI Design Experience [electronic Resource] / by Weng Fook Lee. 1st ed. 2019. Springer International Publishing; 2019. doi:10.1007/978-3-030-03238-8

• A glitch is an unwanted pulse at the output of a combinational logic network – a momentary change in an output that should not have changed
• A circuit with the potential for a glitch is said to have a hazard
• In other words a hazard is something intrinsic about a circuit; a circuit with hazard may or may not have a glitch depending on input patterns and the electric characteristics of the circuit.

When do circuits have hazards ?

Hazards are potential unwanted transients that occur in the output when different paths from input to output have different propagation delays

Types of Hazards (on an output)

static 1-hazard, static 0-hazard, dynamic hazard

Hazard's Concern

• Hazards do not hurt synchronous circuits
• Hazards Kill Asynchronous Circuits
• Glitches Increase Power Consumption

referece

CPE166/EEE 270 Advanced Logic Design-Digital Design: Time Behavior of Combinational Networks: https://www.csus.edu/indiv/p/pangj/166/f/sram/Handout_Hazard.pdf

John Knight, ELEC3500 Glitches and Hazards in Digital Circuits http://www.doe.carleton.ca/~shams/ELEC3500/hazards.pdf

• Arithmetic operations have a potential to run into a condition known as overflow.
• Overflow occurs with respect to the size of the data type that must accommodate the result.
• Overflow indicates that the result was too large or too small to fit in the original data type.

Overflow when adding unsigned number

When two unsigned numbers are added, overflow occurs if

• there is a carry out of the leftmost bit.

Overflow when adding signed numbers

When two signed 2's complement numbers are added, overflow is detected if:

1. both operands are positive and the result is negative, or
2. both operands are negative and the result is positive.

Notice that when operands have opposite signs, their sum will never overflow. Therefore, overflow can only occur when the operands have the same sign.

carryout information ISN'T needed to detect overflow/underflow for signed number addition

EXTBIT:MSB

extended 1bit and msb bit can be used to detect overflow or underflow

1
2
3
4
5
6
7
8
9
10

reg signed  [1:0]      acc_inc;
reg signed  [10-1:0]   acc;
wire signed [10  :0]   acc_w;  // extend 1b for saturation
wire signed [10-1:0]   acc_stat;

assign acc_w = acc + acc_inc;   // signed arithmetic

assign acc_stat = (acc_w[10-1 +: 2] == 2'b01) ? {1'b0, {(10-1){1'b1}}} : // up saturation
    (acc_w[10-1 +: 2] == 2'b10) ? {1'b1, {(10-1){1'b0}}}  :   			 // down saturation
    acc_w[10-1:0];

2'b01 : overflow, up saturation

2'b10: underflow, down saturation

2's complement negative number

1. Flip all bits
2. Add 1.

N-bit signed number \[ A = -M_{N-1}2^{N-1}+\sum_{k=0}^{N-2}M_k2^k \] Flip all bits \[\begin{align} A_{flip} &= -(1-M_{N-1})2^{N-1} +\sum_{k=0}^{N-2}(1-M_k)2^k \\ &= M_{N-1}2^{N-1}-\sum_{k=0}^{N-2}M_k2^k -2^{N-1}+\sum_{k=0}^{N-2}2^k \\ &= M_{N-1}2^{N-1}-\sum_{k=0}^{N-2}M_k2^k -1 \end{align}\]

Add 1 \[\begin{align} A_- &= A_{flip}+1 \\ &= M_{N-1}2^{N-1}-\sum_{k=0}^{N-2}M_k2^k \\ &= -A \end{align}\]

reference

Overflow Detection: http://www.c-jump.com/CIS77/CPU/Overflow/lecture.html

Sign Extension

1. Calculate the necessary minimum width of the sum so that it contains all input possibilities
2. Extend the inputs' sign bits to the width of the answer
3. Add as usual
4. Ignore bits that ripple to the left of the answer's MSB

1. signed

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

module tb;
	reg signed [3:0] inA;
	reg signed [3:0] inB;
    reg signed [4:0] outSumSg;	// signed result
    reg [4:0] outSumUs;			// unsigned result

	initial begin
		inA = 4'b0101;
		inB = 4'b1111;
		outSumUs = inA + inB;
		outSumSg = inA + inB;

		$display("signed out(%%d): %0d", outSumSg);
		$display("signed out(%%b): %b", outSumSg);

		$display("unsigned out(%%d): %0d", outSumUs);
		$display("unsigned out(%%b): %b", outSumUs);
	end
endmodule

2. mixed

1
2

reg signed [3:0] inA;
reg [3:0] inB;

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

module tb;
	reg [3:0] inA;
	reg signed [3:0] inB;
	reg signed [4:0] outSumSg;
	reg [4:0] outSumUs;

	initial begin
		inA = 4'b0101;
		inB = 4'b1111;
		outSumUs = inA + inB;
		outSumSg = inA + inB;

		$display("signed out(%%d): %0d", outSumSg);
		$display("signed out(%%b): %b", outSumSg);

		$display("unsigned out(%%d): %0d", outSumUs);
		$display("unsigned out(%%b): %b", outSumUs);
	end
endmodule

xcelium

1
2
3
4
5
6
7

xcelium> run
signed out(%d): -12
signed out(%b): 10100
unsigned out(%d): 20
unsigned out(%b): 10100
xmsim: *W,RNQUIE: Simulation is complete.
xcelium> exit

vcs

1
2
3
4
5
6

Compiler version S-2021.09-SP2-1_Full64; Runtime version S-2021.09-SP2-1_Full64;  May  7 17:24 2022
signed out(%d): -12
signed out(%b): 10100
unsigned out(%d): 20
unsigned out(%b): 10100
           V C S   S i m u l a t i o n   R e p o r t

observation

When signed and unsigned is mixed, the result is by default unsigned.

Prepend to operands with 0s instead of extending sign, even though the operands is signed

LHS DONT affect how the simulator operate on the operands but what the results represent, signed or unsigned

Therefore, although outSumUs is declared as signed, its result is unsigned

subtraction example

In logic arithmetic, addition and subtraction are commonly used for digital design. Subtraction is similar to addition except that the subtracted number is 2's complement. By using 2's complement for the subtracted number, both addition and subtraction can be unified to using addition only.

operands are signed

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

module tb;
	reg signed [3:0] inA;
	reg signed [3:0] inB;
	reg signed [4:0] outSubSg;
	reg [4:0] outSubUs;

	initial begin
		inA = 4'b1111;
		inB = 4'b0010;
		outSubUs = inA - inB;
		outSubSg = inA - inB;

		$display("signed out(%%d): %0d", outSubSg);
		$display("signed out(%%b): %b", outSubSg);

		$display("unsigned out(%%d): %0d", outSubUs);
		$display("unsigned out(%%b): %b", outSubUs);
	end
endmodule

1
2
3
4
5
6

Compiler version S-2021.09-SP2-1_Full64; Runtime version S-2021.09-SP2-1_Full64;  May  7 17:46 2022
signed out(%d): -3
signed out(%b): 11101
unsigned out(%d): 29
unsigned out(%b): 11101
           V C S   S i m u l a t i o n   R e p o r t

1
2
3
4
5
6

xcelium> run
signed out(%d): -3
signed out(%b): 11101
unsigned out(%d): 29
unsigned out(%b): 11101
xmsim: *W,RNQUIE: Simulation is complete.

operands are mixed

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

module tb;
	reg signed [3:0] inA;
	reg [3:0] inB;
	reg signed [4:0] outSubSg;
	reg [4:0] outSubUs;

	initial begin
		inA = 4'b1111;
		inB = 4'b0010;
		outSubUs = inA - inB;
		outSubSg = inA - inB;

		$display("signed out(%%d): %0d", outSubSg);
		$display("signed out(%%b): %b", outSubSg);

		$display("unsigned out(%%d): %0d", outSubUs);
		$display("unsigned out(%%b): %b", outSubUs);
	end
endmodule

1
2
3
4
5
6

Compiler version S-2021.09-SP2-1_Full64; Runtime version S-2021.09-SP2-1_Full64;  May  7 17:50 2022
signed out(%d): 13
signed out(%b): 01101
unsigned out(%d): 13
unsigned out(%b): 01101
           V C S   S i m u l a t i o n   R e p o r t 

1
2
3
4
5
6
7

xcelium> run
signed out(%d): 13
signed out(%b): 01101
unsigned out(%d): 13
unsigned out(%b): 01101
xmsim: *W,RNQUIE: Simulation is complete.
xcelium> exit

Danger Sign

https://projectf.io/posts/numbers-in-verilog/

Verilog has a nasty habit of treating everything as unsigned unless all variables in an expression are signed. To add insult to injury, most tools won’t warn you if signed values are being ignored.

If you take one thing away from this post:

Never mix signed and unsigned variables in one expression!

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

module signed_tb ();
    logic [7:0] x;  // 'x' is unsigned
    logic signed [7:0] y;  // 'y' is signed
    logic signed [7:0] x1, y1;
    logic signed [3:0] move;

    always_comb begin
        x1 = x + move;  // !? DANGER: 'x' is unsigned but 'move' is signed
        y1 = y + move;
    end

    initial begin
        #10
        $display("Coordinates (7,7):");
        x = 8'd7;
        y = 8'd7;
        #10
        $display("x : %b  %d", x, x);
        $display("y : %b  %d", y, y);

        #10
        $display("Move +4:");
        move = 4'sd4;  // signed positive value
        #10
        $display("x1: %b  %d  *LOOKS OK*", x1, x1);
        $display("y1: %b  %d", y1, y1);

        #10
        $display("Move -4:");
        move = -4'sd4;  // signed negative value
        #10
        $display("x1: %b  %d  *SURPRISE*", x1, x1);     // 0000_0111 + {0000}_1100 = 0001_0011
        $display("y1: %b  %d", y1, y1);                 // 0000_0111 + {1111}_1100 = 0000_0011
    end
endmodule

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

Chronologic VCS simulator copyright 1991-2021
Contains Synopsys proprietary information.
Compiler version S-2021.09-SP2-2_Full64; Runtime version S-2021.09-SP2-2_Full64;  Nov 19 11:02 2022                                     
Coordinates (7,7):
x : 00000111    7
y : 00000111            7
Move +4:
x1: 00001011           11  *LOOKS OK*
y1: 00001011           11
Move -4:
x1: 00010011           19  *SURPRISE*
y1: 00000011            3
           V C S   S i m u l a t i o n   R e p o r t
Time: 60
CPU Time:      0.260 seconds;       Data structure size:   0.0Mb

reference

Lee WF. Learning from VLSI Design Experience [electronic Resource] / by Weng Fook Lee. 1st ed. 2019. Springer International Publishing; 2019. doi:10.1007/978-3-030-03238-8

Bevan Baas, VLSI Digital Signal Processing, EEC 281 - VLSI Digital Signal Processing

With implict sign extension, the implementation of signed arithmetic is DIFFERENT from that of unsigned. Otherwise, their implementations are same.

The implementations manifest the RTL's behaviour correctly

add without implicit sign extension

unsigned

rtl

1
2
3
4
5
6
7

module TOP (
	input wire [2:0] data0
	,input wire [2:0] data1
	,output wire [2:0] result
);
	assign result = data0 + data1;
endmodule

synthesized netlist

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

/////////////////////////////////////////////////////////////
// Created by: Synopsys DC Ultra(TM) in wire load mode
// Version   : S-2021.06-SP5
// Date      : Sat May  7 11:43:27 2022
/////////////////////////////////////////////////////////////


module TOP ( data0, data1, result );
  input [2:0] data0;
  input [2:0] data1;
  output [2:0] result;
  wire   n4, n5, n6;

  an02d0 U6 ( .A1(data0[0]), .A2(data1[0]), .Z(n5) );
  nr02d0 U7 ( .A1(data0[0]), .A2(data1[0]), .ZN(n4) );
  nr02d0 U8 ( .A1(n5), .A2(n4), .ZN(result[0]) );
  ad01d0 U9 ( .A(data1[1]), .B(data0[1]), .CI(n5), .CO(n6), .S(result[1]) );
  xr03d1 U10 ( .A1(n6), .A2(data0[2]), .A3(data1[2]), .Z(result[2]) );
endmodule

vcs compile with -v /path/to/lib.v

signed

rtl

1
2
3
4
5
6
7

module TOP (
	input wire signed [2:0] data0
	,input wire signed [2:0] data1
	,output wire signed [2:0] result
);
	assign result = data0 + data1;
endmodule

synthesized netlist

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

/////////////////////////////////////////////////////////////
// Created by: Synopsys DC Ultra(TM) in wire load mode
// Version   : S-2021.06-SP5
// Date      : Sat May  7 11:48:54 2022
/////////////////////////////////////////////////////////////


module TOP ( data0, data1, result );
  input [2:0] data0;
  input [2:0] data1;
  output [2:0] result;
  wire   n4, n5, n6;

  an02d0 U6 ( .A1(data0[0]), .A2(data1[0]), .Z(n5) );
  nr02d0 U7 ( .A1(data0[0]), .A2(data1[0]), .ZN(n4) );
  nr02d0 U8 ( .A1(n5), .A2(n4), .ZN(result[0]) );
  ad01d0 U9 ( .A(data1[1]), .B(data0[1]), .CI(n5), .CO(n6), .S(result[1]) );
  xr03d1 U10 ( .A1(n6), .A2(data0[2]), .A3(data1[2]), .Z(result[2]) );
endmodule

add WITH implicit sign extension

unsigned with 0 extension

rtl

1
2
3
4
5
6
7

module TOP (
    input wire [2:0] data0		// 3 bit unsigned
    ,input wire [1:0] data1		// 2 bit unsigned
    ,output wire [2:0] result	// 3 bit unsigned
);
	assign result = data0 + data1;
endmodule

synthesized netlist

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

/////////////////////////////////////////////////////////////
// Created by: Synopsys DC Ultra(TM) in wire load mode
// Version   : S-2021.06-SP5
// Date      : Sat May  7 12:15:58 2022
/////////////////////////////////////////////////////////////


module TOP ( data0, data1, result );
  input [2:0] data0;
  input [1:0] data1;
  output [2:0] result;
  wire   n4, n5, n6;

  an02d0 U6 ( .A1(data1[0]), .A2(data0[0]), .Z(n6) );
  ad01d0 U7 ( .A(data1[1]), .B(data0[1]), .CI(n6), .CO(n4), .S(result[1]) );
  xr02d1 U8 ( .A1(data0[2]), .A2(n4), .Z(result[2]) );
  nr02d0 U9 ( .A1(data1[0]), .A2(data0[0]), .ZN(n5) );
  nr02d0 U10 ( .A1(n6), .A2(n5), .ZN(result[0]) );
endmodule

signed with implicit sign extension

rtl

1
2
3
4
5
6
7

module TOP (
	input wire signed [2:0] data0
	,input wire signed [1:0] data1
	,output wire signed [2:0] result
);
	assign result = data0 + data1;
endmodule

synthesized netlist

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

/////////////////////////////////////////////////////////////
// Created by: Synopsys DC Ultra(TM) in wire load mode
// Version   : S-2021.06-SP5
// Date      : Sat May  7 12:21:51 2022
/////////////////////////////////////////////////////////////


module TOP ( data0, data1, result );
  input [2:0] data0;
  input [1:0] data1;
  output [2:0] result;
  wire   n6, n7, n8, n9, n10;

  nd02d0 U9 ( .A1(data1[0]), .A2(data0[0]), .ZN(n10) );
  inv0d0 U10 ( .I(n10), .ZN(n9) );
  nr02d0 U11 ( .A1(data0[1]), .A2(data1[1]), .ZN(n7) );
  aor221d1 U12 ( .B1(n9), .B2(data1[1]), .C1(n10), .C2(data0[1]), .A(n7), .Z(
        n6) );
  xn02d1 U13 ( .A1(data0[2]), .A2(n6), .ZN(result[2]) );
  ora21d1 U14 ( .B1(data1[0]), .B2(data0[0]), .A(n10), .Z(result[0]) );
  aor21d1 U15 ( .B1(data1[1]), .B2(data0[1]), .A(n7), .Z(n8) );
  mx02d0 U16 ( .I0(n10), .I1(n9), .S(n8), .Z(result[1]) );
endmodule

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

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([1 2 3 4 5 6])

ans =

     4     5     6     1     2     3

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

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

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

Neil Robertson, Evaluate Window Functions for the Discrete Fourier Transform

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

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

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

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