The Power Spectrum

Neil RobertsonOctober 8, 2016

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$$

This article is available in PDF format for easy printing

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 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:

$$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πfcnTs), with A = sqrt(2), fc = 5 Hz,  fs = 1/Ts = 32 Hz, and N = 32.  The power into 1 ohm of the analog version of this sinusoid is A2/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πf1nTs) + A* sin(2πf2nTs), with A = sqrt(2), f1 = 125 Hz, f2 = 135.25 Hz, fs = 1/Ts = 1024 Hz, and N = 2048. 

The power into 1 ohm of the analog version of each sinusoidal component is A2/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 f1 and 2 samples on either side; likewise for f2.    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 (2nd 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 (2nd 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:

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.

Sign up
or Sign in