# The Power Spectrum

Often, when calculating the spectrum of a sampled signal, we are interested in relative powers, and we don’t care about the absolute accuracy of the y axis. However, when the sampled signal represents an analog signal, we sometimes need an accurate picture of the analog signal’s power in the frequency domain. This post shows how to calculate an accurate power spectrum.

Parseval’s theorem [1,2] is a property of the Discrete Fourier Transform (DFT) that states:

$$\sum_{n=0}^{N-1} |x(n)|^2 = \frac1N\sum_{k=0}^{N-1}|X(k)|^2$$

where

- x(n) is a time sequence (real or complex)
- X(k) is the discrete Fourier transform of x(n)
- n = 0:N-1 is the time index of x
- k = 0:N-1 is the frequency index of X

Multiply both sides of the above by 1/N:

$$\frac1N\sum_{n=0}^{N-1} |x(n)|^2 = \frac1{N^2}\sum_{k=0}^{N-1}|X(k)|^2$$

|x(n)|^{2} is the instantaneous power of a sample of
the time signal into a 1-ohm load. So the
left side of the equation is just the average power of the signal over the N
samples. Thus, for a 1-ohm load,

$$P_{av}= \frac1{N^2}\sum_{k=0}^{N-1}|X(k)|^2 \quad{watts}$$

For the each DFT bin, we can say:

$$P_{bin}(k)= \frac1{N^2}|X(k)|^2, \quad{k=0:N-1} \quad{watts/bin}\quad{(1)}$$

This is the power spectrum of the signal. For a load resistance R, just divide equation 1 by R.

Note that X(k) is the two-sided spectrum. If x(n) is real, then X(k) is symmetric about f_{s}/2, with each side containing half of the power. In that case, we can choose to keep just the one-sided spectrum, and multiply P_{bin} by 2:

$$P_{bin}(k) \quad{one-sided} = \frac2{N^2}|X(k)|^2, \quad{k=0:N/2-1} \quad{watts/bin}\quad{(2)}$$

Here are a couple of examples that compute the power spectra for sinusoidal signals.

**Example 1**. Power spectrum of a sinusoid with frequency at FFT bin center.

Let x = A*sin(2πf_{c}nT_{s}), with A = sqrt(2), fc = 5 Hz, f_{s} = 1/T_{s} = 32 Hz, and N = 32. The power into 1 ohm of the analog version of this sinusoid is A^{2}/2 = 1 watt.

Here is the Matlab code to compute the power spectrum:

% power_spectrum1.m fs= 32; % Hz sample rate N= 32; fc= 5; % Hz freq of sine wave n= 0:N-1; t= n/fs; % s A= sqrt(2); % Volts x= A*sin(2*pi*fc*t); % time signal % Compute Power Spectrum k= 0:N-1; % freq index f= k*fs/N; % Hz X= fft(x); Pbin= abs(X).^2/N^2; % watts/bin power spectrum % Plot 2-sided spectrum stem(f,Pbin),grid axis([0 fs 0 .6]); xlabel('Hz'),ylabel('Watts'),figure Pdisplay= 2*Pbin(1:N/2); % one-sided spectrum watts/bin % Plot one-sided spectrum stem(f(1:N/2),Pdisplay),grid axis([0 fs/2 0 1.1]); xlabel('Hz'),ylabel('Watts')

The two-sided and one-sided spectra for this (very simple)
example are shown below. As expected, the
one-sided plots shows a power of 1 watt at f = 5 Hz.

**Figure 1.**Two-sided power spectrum

**Figure 2.** One-sided power spectrum

**Example 2.** Power spectrum of two sinusoids – one at FFT bin center and one between bins.

Let x = A*sin(2πf_{1}nT_{s}) + A* sin(2πf_{2}nT_{s}), with A = sqrt(2), f_{1} = 125 Hz, f2 = 135.25 Hz, f_{s} = 1/T_{s} = 1024 Hz, and N = 2048.

The power into 1 ohm of the analog version of each sinusoidal component is A^{2}/2 = 1 watt. The component at 125 Hz falls on a bin center, while the component at 135.25 Hz falls between the bins (bin spacing is 1024 Hz/2048 = 0.5 Hz).

To prevent spectral leakage of the component at 135.25 Hz, we apply a window [3,4,5] to the time sequence. To preserve power accuracy, we need to normalize the window function’s power to 1. Given window function w(n), we multiply by a constant C such that:

$$\frac1{N}\sum_{n=0}^{N-1}(Cw(n))^2 = 1$$

Solving for C, we get:

$$C= \sqrt{ \frac N{\sum_{n=0}^{N-1}(w(n))^2} }$$

Here is the Matlab code to compute the power spectrum. It includes a normalized Hanning window.

% power_spectrum2.m 10/7/16 nr % Plot the power spectrum of two sines with amplitude of sqrt(2) volts. % P = A^2/2 = 1 watt. % one sine frequency is in center of bin and one is between bins % use hanning window fs= 1024; % Hz sample frequency N= 2048; f1= 125; % Hz f1 falls on bin-center f2= 135.25; % Hz f2 falls between bins n= 0:N-1; % time index t= n/fs; % s A= sqrt(2); % Volts x= A*sin(2*pi*f1*t) + A*sin(2*pi*f2*t); % time signal w= hanning(N)'; % hanning window w= w.*sqrt(N/sum(w.^2)); % normalize window for P= 1; xw= x.*w; % multiply signal by window % Compute Power Spectrum k= 0:N/2 -1; % freq index of one-sided spectrum f= k*fs/N; % Hz X= fft(xw); Pbin= 2*abs(X(1:N/2)).^2/N^2; % W/bin factor of 2 for one-sided spectrum PbindB= 10*log10(abs(Pbin)); plot(f,PbindB,'.-'),grid axis([115 145 -50 5]) xlabel('Hz'),ylabel('dBW')

The one-sided spectrum plotted around the signal frequencies is shown in Figure 3. Why are the peaks not at 0 dBW (1 watt)? Windowing in time causes each of our two spectral components to be convolved with the spectrum of the window. The power is thus spread out over the window’s bandwidth, so the peak is less than the total power. To compute the power of each component, we have to sum the power over the window’s bandwidth. In the code below, we have summed the power of the frequency sample at f_{1} and 2 samples on either side; likewise for f_{2}. The result is plotted in Figure 4 as the red circles, which show 0 dBW at each frequency.

We should also note that the scalloping effect of the Hanning window is visible – the peak at 135.25 Hz is 1.42 dB below that at 125 MHz. Although this scalloping has not harmed the power calculation, it is visually somewhat confusing. Note the scalloping could be reduced by zero-padding the windowed time signal, but this has not been done here [6].

Finally, these examples have used sinusoids, which are deterministic signals. To display the spectra of random signals, we would need to add FFT-averaging to the techniques discussed here.

% Post-Processing % Compute power of components at f1 and f2 k1= round(f1*N/fs); % index of f1 P1= sum(Pbin(k1-2:k1+2)) % W sum of spectrum components around f1 P1dB= 10*log10(P1); % dBW k2= round(f2*N/fs); % index of f2 P2= sum(Pbin(k2-2:k2+2)) % W sum of spectrum components around f2 P2dB= 10*log10(P2); % dBW plot(f,PbindB,'.-',[f1 f2],[P1dB P2dB],'or'),grid axis([115 145 -50 5]) xlabel('Hz'),ylabel('dBW')

**Figure 3.** Power Spectrum of sinusoids using windowing.

Peaks are lower than the actual power of 0 dBW.

**Figure 4.** Power Spectrum of sinusoids using windowing.

Red circles = Power of each spectral component

**References**

1. Sanjit K. Mitra, __Digital Signal Processing (2 ^{nd} Ed)__, McGraw-Hill, 2001, Table 3.5.

2. https://en.wikipedia.org/wiki/Parseval%27s_theorem, subsection *Notation Used in Engineering and Physics.*

3. Richard G. Lyons, __Understanding Digital Signal Processing (2 ^{nd} Ed)__ , Prentice Hall, 2004, Section 3.9.

4. https://en.wikipedia.org/wiki/Window_function

5. Segio Rapuano and Fred J. Harris, “An Introduction to FFT and Time Domain Windows”, IEEE Instrumentation and Measurement Magazine, Dec 2007, p. 32.

6. Lyons, sections 3.10, 3.11.

Oct 8, 2016 Neil Robertson

**Previous post by Neil Robertson:**

Digital PLL's -- Part 2

**Next post by Neil Robertson:**

Matlab Code to Synthesize Multiplierless FIR Filters

- Comments
- Write a Comment Select to add a comment

To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.

Registering will allow you to participate to the forums on ALL the related sites and give you access to all pdf downloads.