Hi Raf, you should do this:
should i multiply power values * frequency then add up the result
BUT you want to multiply by the width in frequency, or your delta f. What is the different in Hz between adjacent DFT bins? You can use diff() on your frequency vector to get that info. Multiplying your power estimates by that delta f and summing the result approximates the integral under the PSD.
I don't know how long this vector is in samples (250 msec), but you can check yourself by just using the avgpower method on a single 250 msec sample
x = randn(1e3,1);
psdestx = psd(spectrum.periodogram,x,'Fs',1,'NFFT',length(x));
pwr = avgpower(psdestx,[1/4 1/2]);
It's probably more meaningful to take the frequency interval width into account as you suggest in your second option.
t = 0:0.001:1-0.001;
x = 2*cos(2*pi*100*t);
psdestx = psd(spectrum.periodogram,x,'Fs',1e3,'NFFT',length(x));
pwr = avgpower(psdestx);
Answer is 4/2 as you would expect, now just look in the range around 100 Hz.
pwr = avgpower(psdestx,[95 105])
Still 2, because all the power is there!