CODEWITHSUNDEEP
matlab
SaraA
SaraA

Reputation: 31

Calculating the area between two curves

I have a PSD plot and I am trying to calculate the area between the two curves in MATLAB for two separate ranges of frequencies (8-30 Hz and 70-100 Hz). One area calculation per frequency range.

To calculate the area between the curves for all sections in each frequency range I use:

r1=(f>=8 & f<=30);
r2=(f>=70 & f<=100);
areaLF = trapz(f(r1),Z(r1)-Zm(r1));
areaHF = trapz(f(r2),Z(r2)-Zm(r2));

f=frequency and Zm, Z represent Z-scores of the two conditions (normalized power).

I can calculate the area between the curves for all areas between the curves in each frequency range, but I only want to calculate the area in the low frequency range for when Zm < Z and when Z > Zm for higher frequencies (the shaded regions in the plot).

Does someone have a suggestion about how to do this?

This is the plot showing the areas I am trying to calculate: enter image description here

Upvotes: 1

Views: 719

Answers (1)

TroyHaskin
TroyHaskin

Reputation: 8401

If I understand your question correctly, you need to change your logical masks from being defined by f to being defined by Z and Zm:

r1= Zm < Z;
r2= Z  < Zm;
area1 = trapz(f(r1),Z(r1)-Zm(r1));
area2 = trapz(f(r2),Z(r2)-Zm(r2));

For discontiguous blocks like the example you posted, you may need to integrate over the blocks one-by-one. For that case, I threw together this code (I'm not sure if it is robust or if Matlab has a simple function to perform the task):

%   Predicate
mask    = Zm < Z;

%   Determine index windows for the contiguous blocks.
ind     = 1:numel(f);
dmask   = xor(mask(1:end-1),mask(2:end));
iSwitch = ind(dmask)+1;
nSwitch = nnz(dmask);

if mask(1) == true  % start true
    iLo = [ind(1) , iSwitch(2:2:end)];
    iHi = iSwitch(1:2:end)-1;

    if mod(nSwitch,2) == 0 %    end true
        iHi = [ iHi , ind(end)];
    end

else  % start false
    iLo = iSwitch(1:2:end);
    iHi = iSwitch(2:2:end)-1;

    if mod(nSwitch,2) ~= 0 %    end true
        iHi = [ iHi , ind(end)];
    end
end
nBlocks = numel(iLo);


%   Iterate over blocks
total = 0;
for k = 1:nBlocks;
    mask  = iLo(k):iHi(k);
    total = total + trapz(f(mask),Z(mask) - Zm(mask));
end

Upvotes: 1

Related Questions