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

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?