This example shows how to measure the power of deterministic periodic signals. Although continuous in time, periodic deterministic signals produce discrete power spectra.

We will provide two examples of how to measure a signal's average power. The examples will use sine waves and assume a load impedance of 1 Ohm.

In general signals can be classified into three broad categories, power signals, energy signals, or neither. Deterministic signals which are made up of sinusoids is an example of power signals which have infinite energy but finite average power. Random signals also have finite average power and fall into the category of power signals. A transient signal is an example of energy signals which start and end with zero amplitude. There are still other signals which can't be characterized as either a power or energy signal.

In our first example we'll estimate the average power of a sinusoidal signal with a peak amplitude of 1 volt and a frequency component at 128Hz.

Fs = 1024; t = 0:1/Fs:1-(1/Fs); A = 1; % Vpeak F1 = 128; % Hz x = A*sin(2*pi*t*F1);

Let's look at a portion of our signal in the time domain.

idx = 1:128; plot(t(idx),x(idx)); grid; ylabel('Amplitude'); xlabel('Time (sec)'); axis tight;

The theoretical average power (mean-square) of each complex sinusoid is A^2/4, which in our example is 0.25 or -6.02dB. So, accounting for the power in the positive and negative frequencies results in an average power of (A^2/4)*2.

power_theoretical = (A^2/4)*2

power_theoretical = 0.5000

in dB the power contained in the positive frequencies only is:

10*log10(power_theoretical/2)

ans = -6.0206

To measure the signal's average power we call periodogram and specify the 'power' option.

periodogram(x, hamming(length(x)), [], Fs, 'centered', 'power'); v = axis; axis([v(1) v(2) -10 -5.5]); % Zoom in Y. hgcf = gcf; hgcf.Color = [1,1,1];

As we can see from the zoomed-in portion of the plot each complex sinusoid has an average power of roughly -6dB.

Another way to calculate a signal's average power is by "integrating" the area under the PSD curve.

figure periodogram(x, hamming(length(x)), [], Fs, 'centered', 'psd'); hgcf = gcf; hgcf.Color = [1,1,1];

One thing to notice in this plot is that the peaks of the spectrum plot do not have the same height as when we plotted the power spectrum. The reason is because when taking Power Spectral Density (PSD) measurements it's the area under the curve (which is the measure of the average power) that matters. We can verify that by calling the bandpower function which uses rectangle approximation to integrate under the curve to calculate the average power.

[Pxx_hamming, F] = periodogram(x, hamming(length(x)), [], Fs, 'psd'); power_freqdomain = bandpower(Pxx_hamming, F, 'psd')

power_freqdomain = 0.5000

Since according to Parseval's theorem the total average power in a sinusoid must be equal whether it's computed in the time domain or the frequency domain, we can verify our signal's estimated total average power by summing up the signal in the time domain.

power_timedomain = sum(abs(x).^2)/length(x)

power_timedomain = 0.5000

For the second example we'll estimate the total average power of a signal containing energy at multiple frequency components: one at DC, one at 100Hz, and another at 200Hz.

Fs = 1024; t = 0:1/Fs:1-(1/Fs); Ao = 1.5; % Vpeak @ DC A1 = 4; % Vpeak A2 = 3; % Vpeak F1 = 100; % Hz F2 = 200; % Hz x = Ao + A1*sin(2*pi*t*F1) + A2*sin(2*pi*t*F2); % Let's look at a portion of our signal. idx = 1:128; plot(t(idx),x(idx)); grid; ylabel('Amplitude'); xlabel('Time (sec)'); hgcf = gcf; hgcf.Color = [1,1,1];

Like the previous example, the theoretical average power of each complex sinusoid is A^2/4. The signal's DC average power is equal to its peak power since it's constant and therefore is given by Ao^2. Accounting for the power in the positive and negative frequencies results in a signal's total average power (sum of the average power of each harmonic component) of Ao^2 + (A1^2/4)*2 + (A2^2/4)*2.

power_theoretical = Ao^2 + (A1^2/4)*2 + (A2^2/4)*2

power_theoretical = 14.7500

By calculating the average power of each unique frequency component in dB we'll see that the theoretical results match the mean-square spectrum plot below.

[10*log10(Ao^2) 10*log10(A1^2/4) 10*log10(A2^2/4)]

ans = 3.5218 6.0206 3.5218

To measure the signal's average power we will once again use the periodogram function to calculate and plot the power spectrum of the signal.

periodogram(x, hamming(length(x)), [], Fs, 'centered', 'power'); v = axis; axis([v(1) v(2) 0 7]); % Zoom in Y hgcf = gcf; hgcf.Color = [1,1,1];

As in the first example, estimating the signal's total average power by "integrating" under the PSD curve we get:

periodogram(x, hamming(length(x)), [], Fs, 'centered', 'psd'); hgcf = gcf; hgcf.Color = [1,1,1];

Note that once again the height of the peaks of the spectral density plot at a specific frequency component may not match the ones of the plot of the power spectrum for reasons noted in the first example.

[Pxx, F] = periodogram(x, hamming(length(x)), [], Fs, 'centered', 'psd'); power_freqdomain = bandpower(Pxx, F, 'psd')

power_freqdomain = 14.7500

Again, we can verify the signal's estimated average power by invoking Parseval's theorem and summing up the signal in the time domain.

power_timedomain = sum(abs(x).^2)/length(x)

power_timedomain = 14.7500

You may have noticed that while the height of the peaks of the power and power spectral density plots are different, they differ by a constant ratio.

Pxx = periodogram(x, hamming(length(x)), [], Fs, 'centered', 'psd'); Sxx = periodogram(x, hamming(length(x)), [], Fs, 'centered', 'power'); plot(F, Sxx ./ Pxx) grid on axis tight; xlabel('Frequency (Hz)'); title('Ratio between Power Spectrum and Power Spectral Density'); ratio = mean(Sxx ./ Pxx)

ratio = 1.3638

This ratio is related to the two-sided equivalent noise bandwidth (ENBW) of the window. You can compute this ratio directly by calling ENBW with the window and its corresponding sample rate.

bw = enbw(hamming(length(x)), Fs)

bw = 1.3638

In the previous sections, power was measured from one or multiple sinusoids having a frequency that coincided with a bin. Peak power estimates are usually less accurate when the signal frequency is out of bin. To see this effect, create a sinusoid with a non-integer number of cycles over a one second period.

Fs = 1024; t = 0:1/Fs:1-(1/Fs); A = 1; F = 20.4; x = A*sin(2*pi*F*t); NFFT = length(x); power_theoretical = 10*log10((A^2/4)*2);

Create a Hamming window and a flat top window.

w1 = hamming(length(x)); w2 = flattopwin(length(x));

Compute the periodogram of `x`

using the Hamming window. Zoom in on the peak.

h1 = figure; hold on stem(F,power_theoretical,... 'BaseValue',-50); [Pxx1,f1] = periodogram(x, w1, NFFT, Fs, 'power'); plot(f1,10*log10(Pxx1)) axis([0 40 -45 0]) legend('Theoretical','Periodogram') xlabel('Frequency (Hz)') ylabel('Power (dB)') title('Periodogram Power Spectrum Estimate') grid on

The peak power estimate is below the theoretical peak, and the frequency of the peak estimate is different than the true frequency.

[Pmax,imax] = max(Pxx1); dPmax_w1 = 10*log10(Pmax) - power_theoretical dFreq = f1(imax) - F

dPmax_w1 = -1.1046 dFreq = -0.4000

**Reduce Amplitude Error with Zero-Padding**

To see why this is happening, compute the periodogram using a larger number of FFT bins.

figure hold on periodogram(x, w1, 100*NFFT, Fs, 'power'); ax = gca; ax.ColorOrderIndex = 2; plot(f1,10*log10(Pxx1),'+') stem(F,power_theoretical,... 'BaseValue',-50) axis([0 40 -40 0]) legend('NFFT = 1024','NFFT = 102400',... 'Theoretical Peak')

In the orignal periodogram, the spectral peak is located between two bins and thus the estimated peak is below the theoretical peak. Increasing the number of FFT bins gives a better picture of the spectrum, although this may be a computationally expensive way to improve peak measurements.

**Reduce Amplitude Error with a Flat Top Window**

Another way to produce a better estimate for the peak amplitude is to use a different window. Compute the periodogram of `x`

using a flat top window.

figure(h1) [Pxx,F1] = periodogram(x, w2, NFFT, Fs, 'power'); plot(F1,10*log10(Pxx)) legend('Theoretical','Hamming','Flat Top')

The flat top window is broad and flat. It produces a peak estimate closer to the theoretical value when `x`

does not contain an integer number of cycles, and hence the spectral peak does not fall exactly on a bin.

dPmax_w2 = 10*log10(max(Pxx)) - power_theoretical

dPmax_w2 = -6.2007e-04

The broader peak that the flat top window produces could be a disadvantage when trying to resolve closely spaced peaks, and the frequency of the measured peak is again different from the frequency of the theoretical peak.

**Reduce Amplitude Error with Reassigned Periodogram**

Now add the 'reassigned' flag to `periodogram`

. Periodogram reassignment uses phase information, which is normally discarded, to reassign the signal to its center of energy. The procedure can result in sharper spectral estimates. Plot the reassigned periodogram of `x`

and zoom in on the peak. Use the Hamming window and the flat top window.

[RPxx1,~,~,Fc1] = periodogram(x, w1, NFFT, Fs, 'power','reassigned'); [RPxx2,~,~,Fc2] = periodogram(x, w2, NFFT, Fs, 'power','reassigned'); figure hold on stem(F,power_theoretical,... 'BaseValue',-40); stem(Fc1,10*log10(RPxx1),'BaseValue',-50); stem(Fc2,10*log10(RPxx2),'BaseValue',-50); legend('Theoretical','Hamming Reassignment',... 'Flattop Reassignment') xlabel('Frequency (Hz)') ylabel('Power (dB)') title('Periodogram Power Spectrum Estimate') axis([19.5 21 -4 -2]) grid on

The reassigned estimates of power are closer to the theoretical value for both windows, with the flat top window producing the best peak measurement.

[RPxx1max,imax1] = max(RPxx1); [RPxx2max,imax2] = max(RPxx2); dPmax_reassign_w1 = 10*log10(RPxx1max) - power_theoretical dPmax_reassign_w2 = 10*log10(RPxx2max) - power_theoretical

dPmax_reassign_w1 = -0.0845 dPmax_reassign_w2 = -1.1131e-05

The frequency estimates are also improved using the reassigned periodogram, with the flat top window again giving the best results.

Fc1(imax1)-F Fc2(imax2)-F

ans = -0.0512 ans = 5.6552e-04

Was this topic helpful?