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



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