Power over the total bandwidth doesn't equal the sum of powers over the sub-bandwidths

I have a signal x with power 10 dB and bandwidth W Hz. I want to calculate the power from the PSD of x over [0, W], [0, W/2], [W/2, W]. The sampling frequency is fs, and the total number of samples in x is U, so that the PSD bin bandwidth is fs/U. The number of bins in W hz is WU/fs, and the number of bins in W/2 is WU/(2fs). I construct the three parts and calculated the power in each part as follows

% The signal PSD
xPSD = (2/(fs*U))*abs(fft(x)).^2;

% number of bins in W Hz
numBin1 = W/(fs/U);

% number of bins in W/2 Hz
numBin2 = (W/2)/(fs/U);

% The signal PSD in [0,W] Hz
xPSD1 = xPSD(1:numBins1+1);
% The signal PSD in [0,W/2] Hz
xPSD2 = xPSD(1:numBins2+1);
% The signal PSD in [W/2,W] Hz
xPSD3 = xPSD(numBins2+2:numBins1+1); %numBins2+2 because bin numBins2+1 was included in xPSD2

% frequency axis for xPSD1
freq1 = 0:fs/U:W;
% frequency axis for xPSD2
freq2 = 0:fs/U:W/2;
% frequency axis for xPSD3
freq3 = W/2+(fs/U):fs/U:W; % I removed the last bin included in xPSD2

%% Calculaing the power in [0,W]
pVal1 = 0
for ii=1:numBins1
    pVal1 = pVal1 + xPSD1(ii)*(fs/U);
end
%% Calculaing the power in [0,W/2]
pVal2 = 0
for ii=1:numBins2
    pVal2 = pVal2 + xPSD2(ii)*(fs/U);
end
%% Calculaing the power in [W/2,W]
pVal3 = 0
for ii=1:numBins2-1
    pVal3 = pVal3 + xPSD3(ii)*(fs/U);
end

In theory, I expected that pVal1 = pVal2+pVal3, more or less, but the difference between the two sides is not negligible (sometimes it's 1 dB). I wonder if I got the indices right, or there is something I am missing!!

Answer 1

[Eric Delgado]

You can't do those operations in "logarithm world". You should convert it to linear using db2pow (if power unit) or db2mag (if volts unit) and then you have to integrated in x axes (using trapz, for example). See example below.

ScreenSpan = 50e+6;

% Wi-fi signal generator (frequency: 2.4 GHz; bandwidth: 20 MHz; power: -10 dBm)
FreqCenter = 2.4e+9;
Bandwidth  = 20e+6;

Start   = FreqCenter - ScreenSpan/2;
Stop    = FreqCenter + ScreenSpan/2;
Span    = Stop-Start;

Samples = 1024;
RBW     = 100e+3; % resolution bandwidth (spectrum analyzer parameter)

% Freq = aCoef * idx + bCoef
aCoef = Span/(Samples-1);
bCoef = Start - aCoef;

xData = linspace(Start, Stop, Samples)'; % Hz
yData = -90*ones(Samples,1) + randn(Samples, 1); % dBm

idx1 = round((FreqCenter-Bandwidth/2 - bCoef)/aCoef);
idx2 = round((FreqCenter+Bandwidth/2 - bCoef)/aCoef);

yData(idx1:idx2) = yData(idx1:idx2)+57;
plot(xData, yData)

idx3 = round((FreqCenter-Bandwidth/2 - bCoef)/aCoef);
idx4 = round((FreqCenter+Bandwidth/2 - bCoef)/aCoef);

yData_mW_Hz = db2pow(yData(idx3:idx4))/RBW; % mW/Hz (Normalitazion)

CHANNELPOWER = pow2db(trapz(xData(idx3:idx4), yData_mW_Hz))

CHANNELPOWER = -9.8315

[MAWE]

Thanks. My mistake was that I was comparing 10*log10(pVal2)+10*log10(pVal3) with 10*log10(pVal1), which is not correct, but in the linear scale pVal1 almost equal pVal2+pVal3

PS: All my code is in the linear scale. The only place where I converted to dB is after I calculated the power in linear scale, because the original signal power was defined in dB.