Area Under the Graph

Hi,

I have these multiple plots in a figure and I am trying to calculate the area under the curve (for different frequency ranges) for each of them for comparison pursposes.

For Instance,

I want to determine the area under the curve for 0-50Hz, 50-100Hz, 100-150Hz, and so on for all the flow rates. To get area under the curve, I thought of assigning the horizontal reference line at the lowest point (around -92dB) in this case, which would be constant for all the cases. I know I can use trapz function for this but I am confused with application of limits in my case. Any suggestions?

Best wishes,

Awais

Answer 1

[Star Strider]

I am not certain what you want.

I would instead use the cumtrapz function for the entire data set, then find the frequencies you want.

Example —

A = rand(4, 60);                                                                % Amplitude Array (Assumes Row Vectors)
F = linspace(0, 600, 60);                                                       % Frequency Vector

area_mtx = cumtrapz(F, A, 2);                                                   % Cumulative Integral
frqs = 0 : 50 : 600;                                                            % Frequency Bands
for k1 = 1:numel(frqs)-1
    idxrng = [find(F >= frqs(k1),1,'first') find(F <  frqs(k1+1),1,'last')];    % Index Range
    freq_band(:,k1) = diff(area_mtx(:,idxrng),[],2);                            % Frequency Segment Integral
end 

It would be best to subtract any constant offset area from the ‘area_mtx’ array first, before the loop.

Experiment to get the result you want.

[Awais Yousaf]

Thank you.

I am a little bit confused with what is happening at the end with indexing.

I am after the areas, let's assume A, B, C, D and so on (image attached), with fixed frequency intervals, for amplitudes from multiple data series.

[Star Strider]

My pleasure.

What were the results when you ran my code with your data?

[Star Strider]

My code required a bit of adjusting to work with your data. (It always helps to have your data from the outset.)

This will measure the integrated area of each segment from the ‘refl’ refernce level (here set at -100 dB) to your data, the plot it in the defined segments. Extend the upper limit of the ‘frqs’ vector if you want more than segments between 0 and 600 Hz. (In your initial post, you limited them to 600 Hz.)

The Code —

d = importdata('50mlmin.txt');
F = d.data(:,1);
A = d.data(:,2);

ofst = -100;                                                                    % Reference Level For Area Calculations
area_mtx = cumtrapz(F, ofst-A);                                                 % Cumulative Integral
frqs = 0 : 50 : 600;                                                            % Frequency Bands
for k1 = 1:numel(frqs)-1
    idxrng = [find(F >= frqs(k1),1,'first') find(F <  frqs(k1+1),1,'last')];    % Index Range
    freq_band(k1,:) = diff(area_mtx(idxrng,:));                                 % Frequency Segment Integral
end 

figure
plot(F, A)
hold on
plot([frqs; frqs], ones(2, numel(frqs)).*ylim', '--k')
hold off
grid
xlim([0 600])
cstr = compose('\\bf%6.1f', freq_band);
% cstr = sprintfc('\\bf%6.1f', freq_band);                                        % Use This If You Do Not Have The ‘compose’ Function
text(frqs(1:end-1)+25, -98*ones(1,numel(frqs)-1), cstr, 'HorizontalAlignment','left', 'VerticalAlignment','middle', 'Rotation', 90)  

The Plot —

I plotted this for the first of your files. Since the second one is in the same format, my code should work for it as well. (I am using R2018b. This should work in all relatively recent releases.)

[Awais Yousaf]

Thanks a lot for that, I will try this now.

So in this case, am I correct to assume that the area measured in each 50Hz segment is the area between the curve and the -100dB reference line?

[Star Strider]

My pleasure.

‘So in this case, am I correct to assume that the area measured in each 50Hz segment is the area between the curve and the -100dB reference line?’

Yes. You can of course set the reference value to whatever you want.

This code will draw it as a reference line, if you need to do that:

figure
plot(F, A)
hold on
plot([frqs; frqs], ones(2, numel(frqs)).*ylim', '--k')
plot(xlim, [1 1]*refl, '-r')
hold off
grid
xlim([0 600])
cstr = compose('\\bf%6.1f', freq_band);
% cstr = sprintfc('\\bf%6.1f', freq_band);                                        % Use This If You Do Not Have The ‘compose’ Function
text(frqs(1:end-1)+25, -98*ones(1,numel(frqs)-1), cstr, 'HorizontalAlignment','left', 'VerticalAlignment','middle', 'Rotation', 90)

The rest of my code is unchanged.

[Awais Yousaf]

The "freq_band" output that is generated by the loop for each flow rate (e.g. 50mlmin.txt file), is there a way to run a loop in a way that I can save that particular output for different flow rates within one place/matrix to plot later?

Example attached.

[Star Strider]

Try this (complete code, with the plots, with a loop to read and analyze each file):

filenames = {'50mlmin.txt'; '100mlmin.txt'};
for k2 = 1:numel(filenames)
    D = dlmread(filenames{k2}, '\t', 1, 0);
    A = D(:,2);
    F = D(:,1);
    
    % A = rand(4, 60);                                                                % Amplitude Array (Assumes Row Vectors)
    % F = linspace(0, 600, 60);                                                       % Frequency Vector
    area_mtx = cumtrapz(F, A, 1);                                                   % Cumulative Integral
    frqs = 0 : 50 : 600;                                                            % Frequency Bands
    for k1 = 1:numel(frqs)-1
        idxrng = [find(F >= frqs(k1),1,'first') find(F <  frqs(k1+1),1,'last')];    % Index Range
        freq_band(k1,:) = diff(area_mtx(idxrng),[],1);                            % Frequency Segment Integral
    end
    freq_band_save{k2,:} = {filenames{k2} freq_band};
    
    refl = -50;                                                 % Reference Line For Integration
    figure
    plot(F, A)
    hold on
    plot([frqs; frqs], ones(2, numel(frqs)).*ylim', '--k')
    plot(xlim, [1 1]*refl, '-r')
    hold off
    grid
    xlim([0 600])
    cstr = compose('\\bf%6.1f', freq_band);
    % cstr = sprintfc('\\bf%6.1f', freq_band);                                        % Use This If You Do Not Have The ‘compose’ Function
    text(frqs(1:end-1)+25, -98*ones(1,numel(frqs)-1), cstr, 'HorizontalAlignment','left', 'VerticalAlignment','middle', 'Rotation', 90)
end

The ‘freq_band_save’ cell array has the frequency band data for each file. The first column in each row is the file name, and the second is the ‘frequency_band’ vector for that file.

You can do the same for any of the other variables my code calculates. If you want to avoid recalculating them, consider using the save (link) function to store whatever variables you want in a .mat file. You can then create a table such as the one you attached.

It has been a while since I originally wrote this. I had to refresh my memory.

Answer 2

[Awais Yousaf]

I tried to run the code for my case, it did not run.

I get "Index exceeds matrix dimensions" error. Before this, I tried it with random vectors and I am not sure about what the loop generates, and what I should be looking at from the loop output.

I have attached the text files for two of the series that I am plotting (figure attached).

I am trying to get the areas (A, B, C, D, E, F, ....) under each of the these series plots for 50Hz intervals, to a common baseline (let's say -100db).