MATLAB: Calculate the area under a section of the curve

area under the curveaucpsdwelch

Hello everyone!
I am trying to extract the power of certain frequency bands from my EDA signal. To this purpose, I computed Welch's power spectral density on the signal, and then I wanted to extract the power of the frequencies of interest by calculating the area under the curve between the intervals w = [0, 0.045], w = [0.045, 0.15], and w = [0.15, 0.25]. Does anyone know how I can extract the AUC just within those windows?
I tried to use trapz( ) on just some portions of the signal, but the results look too weird to be accurate. Here is what I did:
% Total signal
EDA_auc = trapz(EDA_pxx); % total power = 0.1063
% AUC for x = [0 0.045]
% cut the signal and do trapz()
y = EDA_pxx(EDA_pxx >= 0 & EDA_pxx <= 0.045);
VLF_auc = trapz(y); % output = 0.1063
% AUC for x = [0.045 0.15]
y1 = EDA_pxx(EDA_pxx >= 0.045 & EDA_pxx <= 0.15);
LF_auc = trapz(y1); % output = 0
The code here should give you a better understanding of the portions of area I would like to calculate:
% w is frequency and is plotted on the x-axis; EDA_pxx is the result of the Welch's PSD analysis and is plotted on
% the y-axis. 
figure;
subplot(2,1,1)
plot(w, EDA_pxx)
title("Welch's Power Spectral Density")
xlabel('Normalised frequency')
ylabel('Power')
subplot(2,1,2)
plot(w, EDA_pxx);
title('Portions of interest')
xlabel('Normalised frequency')
ylabel('Power')
hold on
idx = w >= 0 & w <= 0.25;
H = area(w(idx), EDA_pxx(idx));
set(H(1),'FaceColor',[1 0.5 0]);
xlim([0 0.4])
xline(0.045, ':r')
xline(0.15, ':r')
hold off
I hope it's clear what I'm asking for.
Thank you very much in advance for any suggetsion!

Best Answer

I am not certain what you want.
Try this: 
D1 = load('EDA_pxx.mat');
EDA_pxx = D1.EDA_pxx;
D2 = load('w.mat');
w = D2.w;
ws(1,:) = [0, 0.045];
ws(2,:) = [0.045, 0.15];
ws(3,:) = [0.15, 0.25];
for k = 1:size(ws,1)
    lv = w >= ws(k,1) & w <= ws(k,2);
    wv{k} = w(lv);
    pxxv{k} = EDA_pxx(lv);
    AUC{k} = trapz(wv{k}, pxxv{k});
end
pchc = 'rgb';
figure
plot(w, EDA_pxx)
hold on
for k = 1:size(ws,1)
    patch([wv{k}; flipud(wv{k})], [pxxv{k}; zeros(size(pxxv{k}))], pchc(k))
end
hold off
grid
fprintf('Area %d = %.3E\n', [1:3; [AUC{:}]])
producing: 
Area 1 = 7.149E-05
Area 2 = 8.660E-05
Area 3 = 6.038E-05
and:
.
Related Question