How to calculate the power of a discrete signal? / Clarification on PSD estimates

Question

I want to calculate the channel power $ P_\mathrm{x}$ of a given discrete and complex signal $x[n]$ (with a length of N) in a given bandwidth $B$.

I'm aware, that I could probably apply a sharp bandpass filter with bandwidth $B$ on $x[n]$ and calculate the average power in the time domain. However, I think, it might be easier to do my calculation on the DFT $X[k]$ of my signal, where I just look at those frequency bins within $B$. Here are the different approaches, that I've read about so far, to calculate the power of a signal:

1. Calculate the average power in the time domain by $$ P_\mathrm{x} = \dfrac{1}{N}\sum_{n=1}^{N}|x[n]|^2$$
2. Use Parseval's theorem to calculate power in frequency domain by $$ P_\mathrm{x} = \dfrac{1}{N^2}\sum_{k=1}^{N}|X[k]|^2$$
3. Sum the power spectrum (modified periodogram?) according to Heinzel et al., p. 15 by $$ P_\mathrm{x} = \sum_{k=1}^{N} \dfrac{|X[k]|^2}{S_1^2} $$ with $$ S_1 = \sum_{j=1}^{N} w[j] $$ where w[j] is the j'th DFT window coefficient.
4. Estimate the PSD according to Heinzel et al., p. 16, sum the result and multiply by the DFT frequency resolution $f_{res}$ by $$ P_\mathrm{x} = \sum_{k=1}^{N} \dfrac{|X[k]|^2}{S_2 \cdot f_s} \cdot f_{res} $$ where $f_s$ is the sampling frequency, $$ f_{res} = f_s / N $$ (for an N-point DFT) and $$ S_2 = \sum_{j=1}^{N} w[j]^2 $$ where $w[j]$ is the $j$'th DFT window coefficient.

All of those approaches give the same result, as long as I use a rectangular window on the signal before applying the DFT, which means $S_1 = S_2 = N$. Unfortunately, as soon as I use a different window where $S_1 \neq S_2$, my results from 1) and 4) are still correct, while my results from 2) and 3) are off.

I understand, that 2) and 3) are mathematically the same for a rectangular window where $S_1 = N$, so my guess is, I have to scale 2) by $S_1$, as well, in case I use a non-rectangular window. However, I don't really understand why 3) is incorrect as it is supposed to compute the power spectrum (before taking the sum). Perhaps, I don't understand the definition of a "power spectrum" yet, but why doesn't the sum give me the total signal power?

Also, while I think, that I understand scaling by $S_1$ in 3), I don't understand scaling with $S_2$ in 4). Yet, 4) is correct and 3) isn't (I guess, 3) is only coincidentally correct for a rectangular window, since then it is mathematically the same as 4)). So how do I correctly interpret the power spectrum (and its sum) used in 3)? I've read, that method 3) is actually called "periodogram", which can be used to estimate the PSD. However, how can this give a PSD? I'm nowhere dividing by a frequency (which would give me the unit of a PSD, which is $V^2/Hz$), and it's quite different from the PSD calculation used in 4). I'm really confused here. Can anyone clarify?

I've created a small Matlab script to illustrate those concepts and help you understand my question better.

%% Create signal
N = 4000; % length of signal
fs = 4000; % sampling freq = 4 kHz
f1 = 1000; % signal freq = 1 kHz
ts = 1/fs;
fres = fs/N; % frequency resolution of N-point DFT 

x = zeros(N,1);
for n = 0:N-1
    x(n+1) = 10 * (cos(2*pi*f1*n*ts) + 1j * sin(2*pi*f1*n*ts));
end
x_mag = abs(x);

%% Window the signal and comoute the window sums for scaling
 window = ones(N,1);
% window = hanning(N);
% window = flattopwin(N,'periodic');
S1 = sum(window);
S2 = sum(window.^2);
x_win = x.*window;

%% Transform the signal to frequency domain
X = fft(x_win,N);
X_mag = abs(X);

%% Plot signal representations in time and frequency domain
figure; subplot(2,1,1); plot(real(x)); hold on; plot(imag(x)); hold off;
subplot(2,1,2); plot(X_mag);

%% Calculate power with built-in Matlab function as a reference
power_matlab_time_domain = bandpower(x)

%% 1) Calculate power by average of instantaneous power in time domain
power_time_domain = 1/N * sum(x_mag.^2)

%% 2) Calculate power by parseval theorem
power_parseval_freq_domain = 1/N^2 * sum(X_mag.^2)

%% 3) Calculate power from power spectrum
power_ps_freq_domain = sum(X_mag.^2/S1^2)

%% 4) Calculate power from power spectral density
power_psd_freq_domain = sum(X_mag.^2/(S2*fs)) * fres

In the end, I'd like to sum only over a range of the discrete frequencies to obtain the channel power of a frequency modulated signal, using either approach 3) or 4). Is this a valid approach?

Answer 1

Let's start with the distinction between calculating the power in a signal, and estimating the power. Calculating the power is straightforward, and you've given the discrete case in equations 1 and 2. However most measurements contain various types of noise, and it is useful to perform estimation. Periodograms are a simple way of doing this (see Understanding the Windowing Method in PSD Calculation for some good details). The basics of it are you split the signal into multiple parts, find the power spectrum of each sub-signal, then combine the multiple estimates to get the final estimate. It sacrifices some frequency resolution by using shorter sub-sequences, but removes quite a bit of the effects of noise.

This is what equations 3 and 4 give(or page 15 and 16 of the linked article, excepting a constant it seems). The first describes how to normalize the frequency of each subsequence, based on the window chosen. The second describes how to normalize the combined subsequence. Lastly, "if the desired result is a power spectral density ( PSD ) expressed in V 2 /Hz, it is obtained by dividing the power spectrum (PS) by the effective noise-equivalent bandwidth ENBW (pg 15 of Heinzel et al.)." In other words, equation 4 should give the resulting spectral density. From a practical perspective a few things such as window type and length need to be chosen based on the problem and computing power available. Then, as you said, summing the frequency bins of interest will give an estimate of the power.

Answer 2

Just notice there is an error on formula under

2. Use Parseval's theorem to calculate power in frequency domain by...

The correct formular is $\frac{1}{N}$, not $\frac{1}{N^{2}}$.