The ESD, PSD and MS Spectrum functions

4.4 The ESD, PSD and MS Spectrum functions

While the PSD $P (f)$ targets power signals, the ESD $G (f)$ is a convenient tool for analyzing energy signals, informing how the total energy $E$ is distributed over frequency. Both are “density” functions, with the independent variable $f \in ℝ$ . The mean-square (MS) spectrum $S_{ms} [k]$ is a discrete-frequency function with the independent variable $k \in ℤ$ . We will discuss each of the three functions in this section.

4.4.1 Energy spectral density (ESD)

To understand the ESD definition, it is useful to recall the Parseval theorem. The (Parseval) Theorem 2 for block transforms can be extended to continuous and discrete-time signals. For continuous-time, together with Eq. (1.22), one has

E = \int_{- \infty}^{\infty} | x (t) |^{2} d t = \int_{- \infty}^{\infty} | X (f) |^{2} d f,

(4.14)

stating that the energy is the same in both time and frequency domains. This suggests defining

G (f) = | X (f) |^{2}

(4.15)

as the ESD in continuous-time. Therefore, $G (f)$ describes how the signal energy $E$ is spread over frequency and its integral is the energy $E$ . The adopted unit for the ESD is joules/Hz.

Example 4.3. ESD of a real-valued exponential. Consider the energy signal $x (t) = e^{- at} u (t)$ in volts, where $a = 0.9$ . Its Fourier transform is $X (f) = 1 ∕ (j 2 πf + a)$ and the ESD in joules/Hz is

G (f) = | X (f) |^{2} = \frac{1}{{(2 πf)}^{2} + a^{2}} .

The signal energy is $E = \int_{- \infty}^{\infty} | x (t) |^{2} d t = \int_{- \infty}^{\infty} G (f) d f \approx 0.55556$ . Listing 4.6 illustrates how $E$ can be estimated in time-domain and also using $G (f)$ . In this case, the range $] - \infty, \infty [$ is simulated as $[0, 30]$ s and $[- 1000, 1000]$ Hz in time and frequency-domain, respectively.

Listing 4.6: MatlabOctaveCodeSnippets/snip_frequency_exp_esd.m

1%% Generate time-domain signal 
2Ts=0.0001; %defines resolution in time 
3t=0:Ts:30; %time axis 
4a=0.9; %constant 
5x=exp(-a*t); %signal 
6subplot(211), plot(t,x) 
7xlabel('t (s)'); ylabel('x(t) amplitude') 
8%% Generate ESD 
9f=linspace(-2,2,1000); %frequency axis 
10Gf= 1 ./ ((2*pi*f).^2 + a^2); %ESD 
11subplot(212) 
12plot(t,x), plot(f,Gf) 
13xlabel('f (Hz)'); ylabel('ESD (J/Hz)') 
14%% Confirm Parseval theorem 
15energy_equation = 1/(2*a); %theoretical 
16energy_time = sum(abs(x).^2)*Ts; %integration 
17df=0.0001; %defines resolution in frequency 
18min_f = -1000; %range of interest 
19max_f = 1000; 
20f=min_f:df:max_f; %frequency axis 
21Gf= 1 ./ ((2*pi*f).^2 + a^2); %recalculate ESD 
22energy_frequency = sum(Gf)*df; %integration 
23disp(['Theoretical energy =' num2str(energy_equation) ' Joules']) 
24disp(['Integration in time =' num2str(energy_time) ' Joules']) 
25disp(['Integration in frequency =' num2str(energy_frequency) ' J']) 
26%% Energy in given frequency band 
27min_f = 0.2; max_f = 1; %range of interest 
28f=min_f:df:max_f; %new frequency axis 
29Gf= 1 ./ ((2*pi*f).^2 + a^2); %recalculate ESD 
30energy_band = 2*sum(Gf)*df; %integration using 2 times 
31disp(['Energy in band =' num2str(energy_band) ' Joules'])

Listing 4.6 also illustrates how $G (f)$ can be conveniently used to estimate the energy within a given frequency band. In this example, the energy within the range from 0.2 to 1 Hz, considering both negative and positive frequencies, is $E_{r} = 2 \int_{0.2}^{1} G (f) d f \approx 0.16954$ J. $□$

Example 4.4. ESD of a sinc in time-domain. Now consider the energy signal is $x (t) = 400 sinc (100 t)$ in volts. According to Eq. (A.56), one needs to consider that $2 F = 100$ , such that $F = 50$ Hz. With this assumption, Eq. (A.56) can be used together with the linearity property of the Fourier transform. The linearity takes care of a factor of 4 that multiplies Eq. (A.56), such that one has $X (f) = 4$ for $| f | \leq 50$ Hz and 0 otherwise. Hence, the ESD of $x (t)$ is $G (f) = 16$ for $| f | \leq 50$ Hz and 0 otherwise. In this case, using the ESD instead of $x (t)$ , it is easy to conclude that the total signal energy is $E = 16 \times 100 = 1600$ joules. $□$

4.4.2 Advanced: Units of ESD when angular frequencies are adopted

Similar to the discussion in Section 2.5.4 regarding the Fourier transform $X (ω)$ with $ω$ in rad/s, there is a subtle issue when using radians per second instead of Hertz for the ESD. In this case, the factor $2 π$ needs to be taken in account such that $G (ω) ∕ (2 π) = | X (ω) |^{2} ∕ (2 π)$ has the unit of joules/(rad/s). The energy over a band of frequencies specified in rad/s can then be calculated by integrating $| X (ω) |^{2} ∕ (2 π)$ over this band. For instance, considering a band $ω \in [- \infty, \infty]$ the energy $E$ is:

E = \int_{- \infty}^{\infty} | x (t) |^{2} d t = \frac{1}{2 π} \int_{- \infty}^{\infty} | X (ω) |^{2} d ω = \frac{1}{2 π} \int_{- \infty}^{\infty} G (ω) d ω,

(4.16)

where $G (ω) = G (2 πf)$ and the unit of $G (ω) ∕ (2 π)$ is joules/(rad/s).

Example 4.5. The unit of $G (ω) ∕ (2 π)$ is joules/(rad/s). For instance, take the sinc signal of Example 4.4, which has a flat ESD given by $G (f) = 16$ J/Hz within the range $[- 50, 50]$ Hz and zero otherwise. The representation $G (ω)$ of this ESD in rad/s has the same constant value $G (ω) = 16$ but now over the range $[- 100 π, 100 π]$ rad/s. The value $G (ω) ∕ (2 π) = 16 ∕ (2 π)$ can be interpreted as a density in J/(rad/s). Within the range $[- 100 π, 100 π]$ , one has $(16 ∕ (2 π)) \times 200 π = 1600$ J. $□$

In continuous-time spectral analysis, the linear frequency in Hertz is more convenient and angular frequencies are seldom adopted. But this is not an option in discrete-time processing, given that $Ω$ is an angle (assumed in radians) and does not have a counterpart in Hertz. Hence, the ESD in discrete-time

G (e^{j Ω}) = | X (e^{j Ω}) |^{2}

(4.17)

needs to be normalized by the factor $2 π$ to be interpreted as $G (e^{j Ω}) ∕ (2 π)$ in units of joules per radians. Using this definition, one can properly interpret a version of the Parseval theorem in discrete-time:

E = \sum_{n = - \infty}^{\infty} | x [n] |^{2} = \frac{1}{2 π} \int_{< 2 π >} | X (e^{j Ω}) |^{2} d Ω .

(4.18)

Table 4.2 summarizes the discussed ESD functions.

Table 4.2: ESD functions.

E

is the total energy and the column “Variable” indicates the units and symbols used for the independent variable of each function. The units of

G (f)

G (ω) ∕ (2 π)

and

G (e^{j Ω}) ∕ (2 π)

are J/Hz, J/(rad/s) and J/rad, respectively.


Time	Variable	ESD definition	ESD main property

Continuous-time	$f$ (Hz)	$G (f) = \| X (f) \|^{2}$	$E = \int_{- \infty}^{\infty} G (f) d f$

Continuous-time	$ω$ (rad/s)	$G (ω) = \| X (ω) \|^{2}$	$E = \frac{1}{2 π} \int_{- \infty}^{\infty} G (ω) d ω$

Discrete-time	$Ω$ (rad)	$G (e^{j Ω}) = \| X (e^{j Ω}) \|^{2}$	$E = \frac{1}{2 π} \int_{< 2 π >} G (e^{j Ω}) d Ω$

4.4.3 Power spectral density (PSD)

In most cases, the signal under analysis has infinite energy, such as a deterministic power signal (e. g. a periodic signal) or realizations of a stationary random process. Therefore, the main interest in spectral analysis relies not on the ESD but on the PSD.

Main property of a PSD

Noticing from Eq. (4.14) and (4.18) that squared magnitudes provide the energy distribution over frequency, to obtain the power distribution one can intuitively consider dividing the ESD by “time”, via a normalization factor that converts energy into power.

The PSD has the important property that the average power $P$ can be obtained in continous-time by

P_{c} = \int_{- \infty}^{\infty} S (f) d f,

(4.19)

with $S (f)$ in watts/Hz.

Similar to the reasoning associated to Eq. (4.16), the version of Eq. (4.19) when the continuous-time PSD $S (ω)$ is a function of $ω = 2 πf$ in radians/s is

P_{c} = \frac{1}{2 π} \int_{- \infty}^{\infty} S (ω) d ω,

(4.20)

and $S (ω) ∕ (2 π)$ is interpreted in units of watts per rad/s.

The PSD definition in discrete-time is

P_{d} = \frac{1}{2 π} \int_{< 2 π >} S (e^{j Ω}) d Ω,

(4.21)

where $S (e^{j Ω}) ∕ (2 π)$ is interpreted in units of watts per radians.

In all three equations above, the PSDs describe the distribution of power over frequency. The frequency can be linear $f$ in Hz, angular frequency $ω$ in radians/s or discrete-time angular frequency $Ω$ , which is an angle in radians.

PSD definitions and the focus on random signals

As discussed in Section 1.6, there are many categories of signals. Figure 4.24 illustrates some of them. Different categories of signals require distinct definitions of PSD. To simplify the discussion, some spectral analysis textbooks³ choose to emphasize discrete-time random signals, and define the PSD for this kind of signals. Moreover, the spectral estimation problem is defined in these textbooks⁴ as:

From a finite-length record ${y [1], \dots, y [N]}$ of a second-order stationary random process, find an estimate $Ŝ (e^{j Ω})$ of its power spectral density $S (e^{j Ω})$ for $Ω \in [- π, π]$ .

In other words, the context of this problem is restricted to discrete-time power signals with duration of $N$ samples. Here we also discuss other cases, such as: a) continuous-time PSDs, b) when the signal is deterministic and c) the theoretical PSD expressions for infinite-duration signals.

Figure 4.24: Important categories of signals. Some of these features need to be taken in account when defining the PSD function.

It makes complete sense to focus the study of PSD estimation on random signals, because they are important in many applications, such as in digital communications. One useful model for these signals is the wide-sense stationary (WSS) random process with autocorrelation $R (τ)$ (see definitions in Appendix A.20). In many spectral analysis problems,⁵ besides being WSS, the random process is assumed to be autocorrelation ergodic.

The PSD $S (f)$ for a continuous-time power signal $x (t)$ corresponding to a realization of a WSS stochastic process $X (t)$ is defined as

S (f) \overset{def}{=} lim_{T \to \infty} \frac{𝔼 [| X_{T} (f) |^{2}]}{T},

(4.22)

where $X_{T} (f)$ is the Fourier transform of a truncated (windowed) version of $x (t)$ with duration $T$ . In other words, $X_{T} (f) = F {x_{T} (t)}$ where $x_{T} (t) = x (t)$ for $- T ∕ 2 \leq t \leq T ∕ 2$ or zero otherwise.

An alternative to obtaining a PSD via Eq. (4.22) is using the Wiener-Khinchin theorem,⁶ which states that the power spectral density (PSD) $S (f)$ of a WSS process is the Fourier transform of the corresponding autocorrelation function, i. e.

S (f) = F {R (τ)} .

(4.23)

Similar expressions exist in discrete-time and even for deterministic signals.

For example, in discrete-time processing, the PSD is the DTFT of the autocorrelation function:

S (e^{j Ω}) = ∑_{ℓ = - \infty}^{\infty} R [ℓ] e^{- j Ω ℓ},

(4.24)

where $ℓ \in ℤ$ is the lag, and the inverse DTFT gives

R [ℓ] = \frac{1}{2 π} \int_{- π}^{π} S (e^{j Ω}) e^{j Ω ℓ} d Ω .

(4.25)

Similar to Eq. (4.22), the PSD $S (e^{j Ω})$ for a discrete-time power signal $x [n]$ corresponding to a realization of a WSS stochastic process $X [n]$ is defined as

S (e^{j Ω}) \overset{def}{=} lim_{N \to \infty} \frac{1}{N} 𝔼 [{| \sum_{n = 0}^{N - 1} x [n] e^{- j Ω n} |}^{2}] = lim_{N \to \infty} \frac{1}{N} 𝔼 [{| X_{N} (e^{j Ω}) |}^{2}] .

(4.26)

where $X_{N} (e^{j Ω})$ is the DTFT of $x_{N} [n]$ , a truncated version of $x [n]$ obtained via a rectangular window of $N$ non-zero samples. The unit of $S (e^{j Ω}) ∕ (2 π)$ is watts per radians.

In practice, there is a finite number of realizations of $X [n]$ and often only one realization $x [n]$ is available. Fortunately, ergodicity of the autocorrelation can be assumed in many cases and the ensemble averages substituted by averages taken over time (see Appendix A.19). Besides, the number of samples (the duration of $x [n]$ ) is often limited to a given value that can be determined, for example, by the time over which the process can be considered stationary. For example, in speech analysis applications, it is typically assumed the process of vowel production is quasi-stationary over segments with durations from 40 to 80 ms. With limited-duration signals, the challenge for spectral analysis is to obtain accurate estimates, as discussed in this chapter.

For both Eq. (4.22) and Eq. (4.26), windows other than the rectangular can be used when the signal under analysis has a short duration. But a rectangular window with infinite duration is the adequate window for the PSD definition.

Table 4.3 summarizes the discussed PSD functions.

Table 4.3: PSD functions.

P

is the average power and the column “Variable” indicates the units and symbols used for the independent variable of each function.

F_{f}

and

F_{ω}

denote the Fourier transform in Hertz and rad/s, respectively. The units of

S (f)

S (ω) ∕ (2 π)

and

S (e^{j Ω}) ∕ (2 π)

are W/Hz, W/(rad/s) and W/rad, respectively.


Time	Variable	Definition	Main property

Continuous-time	$f$ (Hz)	$S (f) = F_{f} {R (τ)}$	$P = \int_{- \infty}^{\infty} S (f) d f$

Continuous-time	$ω$ (rad/s)	$S (ω) = F_{ω} {R (τ)}$	$P = \frac{1}{2 π} \int_{- \infty}^{\infty} S (ω) d ω$

Discrete-time	$Ω$ (rad)	$S (e^{j Ω}) = F {R [ℓ]}$	$P = \frac{1}{2 π} \int_{< 2 π >} S (e^{j Ω}) d Ω$

PSD of deterministic and periodic signals

The PSD definition that is adopted for continuous-time deterministic signals (and does not require the expected value used in Eq. (4.22)) is

S (f) = lim_{T \to \infty} \frac{| X_{T} (f) |^{2}}{T},

which is basically the ESD normalized by the time interval $T$ .

A special case of deterministic signals are the periodic ones. Assuming a continuous-time signal $x (t)$ with period $T_{0}$ , its PSD is

S (f) = \sum_{k = - \infty}^{\infty} | c_{k} |^{2} δ (f - k F_{0}),

(4.27)

where $c_{k}$ are the Fourier Series coefficients and $F_{0} = 1 ∕ T_{0}$ is the fundamental frequency. Similarly, for $S (ω)$ , the expression is

S (ω) = 2 π \sum_{k = - \infty}^{\infty} | c_{k} |^{2} δ (f - k ω_{0}),

(4.28)

where $ω_{0} = (2 π) ∕ T_{0}$ rad/s.

In summary, the PSD $S (f)$ of a deterministic (non-random) periodic signal is composed by impulses with areas determined by the squared magnitude of Fourier series coefficients.

Example 4.6. PSD of a continuous-time sinusoid. If $x (t) = A cos (2 π f_{c} t)$ , then $S (f) = \frac{A^{2}}{4} [δ (f + f_{c}) + δ (f - f_{c})]$ and the average power is $\int_{- \infty}^{\infty} S (f) d f = A^{2} ∕ 2$ . $□$

When considering a discrete-time periodic signal $x [n]$ , its PSD $S (e^{j Ω})$ can be obtained by first considering an expression $S' (e^{j Ω})$ for the frequency range $[0, 2 π [$ :

S' (e^{j Ω}) = 2 π ∑_{k = 0}^{N_{0} - 1} | \tilde{X} [k] |^{2} δ (Ω - k Ω_{0}),

(4.29)

where $N_{0}$ is the period, $Ω_{0} = (2 π) ∕ N_{0}$ is the fundamental frequency, and $\tilde{X} [k]$ the DTFS of $x [n]$ . Finally, the PSD $S (e^{j Ω})$ is simply the periodic repetition of $S' (e^{j Ω})$ :

S (e^{j Ω}) = ∑_{p = - \infty}^{\infty} S' (e^{j (Ω + p 2 π)}) .

(4.30)

Example 4.7. PSD of a discrete-time sinusoid. If $x [n] = A cos (Ω_{1} n)$ (assume $Ω_{1}$ obeys Eq. (1.21) to have $x [n]$ periodic), then $S' (e^{j Ω}) = \frac{π A^{2}}{2} [δ (Ω + Ω_{1}) + δ (Ω - Ω_{1})]$ and

S (e^{j Ω}) = ∑_{p = - \infty}^{\infty} \frac{π A^{2}}{2} [δ (Ω + Ω_{1} + p 2 π) + δ (Ω - Ω_{1} + p 2 π)]

provides its PSD, which has period $2 π$ , as expected. $□$

4.4.4 Advanced: Fourier modulation theorem applied to PSDs

If $S_{x} (f)$ is the PSD of a WSS random process $X (t)$ , the PSD $S_{y} (f)$ of a new process $Y (t) = X (t) e^{j 2 π f_{c} t}$ is

S_{y} (f) = S_{x} (f - f_{c}) .

(4.31)

Similarly, if $Y (t) = X (t) cos (2 π f_{c} t)$ , then

S_{y} (f) = \frac{1}{4} [S_{x} (f + f_{c}) + S_{x} (f - f_{c})] .

(4.32)

To observe why Eq. (4.31) (and, consequently, Eq. (4.32)) is true, recall the Wiener-Khinchin theorem of Eq. (4.23) and the autocorrelation definition of Eq. (1.50). Generally speaking, when modifying a WSS process $X (t)$ , the effect on its PSD can be obtained by checking how the modification affects its autocorrelation, and then relating this to frequency domain using Eq. (4.23). For example, multiplying $X (t)$ by a scalar $α$ corresponds to scaling its autocorrelation $R_{x} (τ)$ by $α^{2}$ and leads to a PSD $α^{2} S_{x} (f)$ given the linearity of the Fourier transform.

According to this reasoning, a proof sketch of Eq. (4.31) follows:

\begin{array}{l} S_{y} (f) & = F {R_{y} (τ)} \\ = F {𝔼 [Y (t + τ) Y^{*} (t)]} \\ = F {𝔼 [X (t + τ) e^{j 2 π f_{c} (t + τ)} X^{*} (t) e^{- j 2 π f_{c} t}]} \\ = F {e^{j 2 π f_{c} τ} 𝔼 [X (t + τ) X^{*} (t)]} \\ = F {e^{j 2 π f_{c} τ} R_{x} (τ)} \\ = S_{x} (f - f_{c}) . \end{array}

Eq. (4.32) can be obtained by decomposing the cosine $cos (2 π f_{c} t) = 1 ∕ 2 (e^{j 2 π f_{c} t} + e^{- j 2 π f_{c} t})$ into two complex exponentials and taking in account that the factor $α = 1 ∕ 2$ leads to the $1 ∕ 4$ in the PSD expression. Eq. (4.32) allows to observe that the result of a signal multiplied by a cosine of unitary amplitude has half of the original signal power.

4.4.5 Mean-square (MS) spectrum

The PSD is very useful especially when dealing with random signals and the power of the signal over a frequency range is obtained by integrating the PSD over that range. However, in some cases it is desired to use a function that allows to directly infer the average power of sinusoid components of a periodic signal, without the integration step. In these cases, the so-called MS or power spectrum is more convenient.⁷ However, in applications characterized by a mixed signal in which there is a deterministic signal of interest that is contaminated by random noise (such as sinusoids contaminated by AWGN), the PSD representation is often more convenient than the MS spectrum.

While in continuous-time the PSD $S (f)$ unit is watts/Hz, the mean-square spectrum is given directly in watts. Assuming a signal with fundamental period $N_{0}$ samples, the mean-square spectrum corresponds to the squared magnitude of the corresponding DTFS:

S_{ms} [k] = | DTFS {x [n]} |^{2},

(4.33)

with the property

\sum_{k = 0}^{N_{0} - 1} S_{ms} [k] = P,

(4.34)

where $P$ is the signal power.

Recall from the discussion associated to Eq. (2.49) that, if $x [n]$ is periodic with fundamental period $N_{0}$ , its DTFS $\tilde{X} [k]$ can be obtained with an $N_{0}$ -point FFT:

\tilde{X} [k] = \frac{FFT {x [n]}}{N_{0}} .

(4.35)

One often uses Eq. (4.35) even for a non-periodic $x [n]$ , but then the result spectrum must be properly interpreted: as if the signal were a periodic version $\sum_{p = - \infty}^{\infty} x_{N} [n - pN]$ of the windowed version $x_{N} [n]$ of $x [n]$ using $N$ samples.

If the FFT size $N$ is chosen as $N = N_{0}$ , the $k$ -th FFT value in Eq. (4.35) corresponds exactly to the $k$ -th DTFS coefficient, that is, they both represent the same frequency $k (2 π ∕ N) = k (2 π ∕ N_{0})$ . In case $N \neq N_{0}$ , Eq. (4.37) is still a valid way to obtain the average power $P$ , but the $k$ -th bin frequency $Ω_{k}$ must be interpreted according to the FFT grid as $Ω_{k} = k (2 π ∕ N)$ .

In general, an estimate $Ŝ_{ms} [k]$ of the MS spectrum of a discrete-time signal $x [n]$ can be obtained from its $N$ -length windowed version $x_{N} [n]$ with

Ŝ_{ms} [k] = {| \frac{FFT {x_{N} [n]}}{N} |}^{2}

(4.36)

and

\sum_{k = 0}^{N - 1} Ŝ_{ms} [k] = P .

(4.37)

The “hat” in $Ŝ_{ms} [k]$ indicates that in general, Eq. (4.36) is an estimate of the true MS spectrum.

³ For instance, see [Hay09, SM05].

⁴ See [SM05], page 13.

⁵ See [Hay09], page 391, Chapter 8.

⁶ See, e. g. [ url3wie] and note that a realization $x (t)$ of a WSS process $X (t)$ is not square integrable and does not have a Fourier transform.

⁷ Matlab (but not Octave) has a msspectrum spectral estimator and the private functions welch.m and computeperiodogram.m can be studied.