Filling in the area between two curves, between specific points in MATLAB

1
votes

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

This is the code I used to generate the plot, where f=frequency and Zm,Z represent Z-scores of the two conditions:

plot(f,Zm,f,Z,'LineWidth',2)
xlim([0 100]);
xlabel('Frequency (Hz)');
ylabel('Normalized Power');

I believe I need to use the trapz function to calculate the area and the fill function to fill in the space, but I am unsure as to how to use these functions to perform calculations between specific frequencies.

To further complicate things, I want to only shade regions where Zm < Z for frequencies 8-30Hz and where Zm > Z for frequencies 70-100Hz.

This is the plot in question:

plot

2
matlab
Relevant: blogs.mathworks.com/graphics/2015/10/13/fill-between . Also, I've made some edits to your question, feel free to change it if you disagree with any of it. – Andras Deak
@AndrasDeak: That was the perfect reference I needed. Worked great! – SaraA
Glad to hear it:) If your final solution significantly differs from the answer of Arzeik, you should consider adding your own as another answer (especially if there's anything non-trivial about it). – Andras Deak
@AndrasDeak: I ended using a combination of the suggestions (including the blog post you posted) and I posted it below in an answer. I am still unsure as to how to make the calculation of the area though. I only want the area for Zm < Z for low frequencies and Z > Zm for higher frequencies. Should I ask this in a separate question? – SaraA
@SaraA yes. Please ask it as a separate question, as you managed to solve the problem outlined in your current question. Also please accept your own answer in that case (there's a 24h delay on self-accepts, so you'll have to do that tomorrow) – Adriaan

2 Answers

0
votes

For calculating the area, something like this may work:

r1=(f>=8).*(f<=30);
r2=(f>=70).*(f<=100);
area=trapz(f(r1),Zm(r1)-Z(r1))+trapz(f(r2),Zm(r2)-Z(r2));

For filling the first region, you can try:

f1=f(r1);
Zm1=Zm(r1);
Z1=Z(r1);
for i=1:length(f1)-1
    fill([f1(i) f1(i+1) f1(i+1) f1(i)],[Z1(i) Z1(i+1) Zm1(i+1) Zm1(i)],'blue');
end

And something similar for the second one.

If you want to fill only the region where Zm is less than Z, you can use the same code but adding an if inside the for loop. Something like:

for i=1:length(f1)-1
    if Zm1(i)<Z1(i) && Zm1(i+1)<Z1(i+1)
        fill([f1(i) f1(i+1) f1(i+1) f1(i)],[Z1(i) Z1(i+1) Zm1(i+1) Zm1(i)],'blue');
    end
end
0
votes

This is the code I ended up using (a combination of the suggestions): 

%// Fill in area between low frequency bands
f1=f(r1)'; %'// <-- this prevents string markdown
Zm1=Zm(r1);
Z1=Z(r1);

mask = Zm1 <= Z1;
fx = [f1(mask), fliplr(f1(mask))];
fy = [Zm1(mask),fliplr(Z1(mask))];
hold on
fill_color = [.040 .040 1];
fh = fill(fx,fy,fill_color);
hold off
%fh.EdgeColor = 'none';

delete(fh)

hold on
output = [];
%// Calculate t in range [0 1]
calct = @(n) (n(3,1)-n(2,1))/(n(3,1)-n(2,1)-n(3,2)+n(2,2));
%// Generate interpolated X and Y values
interp = @(t,n) n(:,1) + t*(n(:,2) - n(:,1));

for i=1:length(f1)
%// If y2 is below y1, then we don't need to add this point.
if Z1(i) <= Zm1(i)
    %// But first, if that wasn't true for the previous point, then add the
    %// crossing.
    if i>1 && Z1(i-1) > Zm1(i-1)
        neighborhood = [f1(i-1), f1(i); Zm1(i-1), Zm1(i); Z1(i-1), Z1(i)];
        t = calct(neighborhood);
        output(:,end+1) = interp(t,neighborhood);
    end
else
%// Otherwise y2 is above y1, and we do need to add this point. But first
%// ...
    %// ... if that wasn't true for the previous point, then add the
    %// crossing.
    if i>1 && Z1(i-1) <= Zm1(i-1)
        neighborhood = [f1(i-1), f1(i); Zm1(i-1), Zm1(i); Z1(i-1), Z1(i)];
        t = calct(neighborhood);
        output(:,end+1) = interp(t,neighborhood);
    end

    %// add this point.
    output(:,end+1) = [f1(i); Z1(i); Zm1(i)];
end
end

xout = output(1,:);
topout = output(2,:);
botout = output(3,:);
fh = fill([xout fliplr(xout)],[botout fliplr(topout)],fill_color);
fh.EdgeColor = 'none';
uistack(fh,'bottom')
hold off

I repeated this for high frequencies.

The resulting plot: enter image description here