Area under FFT curve

I need to calculate the area under curve for the fundamental frequency and some harmonics in an FFT curve. I have attached the FFT curve.

For FFT I used below code:

Fs=125;

L=length(A(1250000:2500000));

n=pow2(nextpow2(L));

f=Fs*(0:1/n:1-(1/n));

df=log10(abs(fft(A(1250000:2500000).*hanning(L),n)));

I need to calculate the area under the curve for a certain bandwidth. For example, the fundamental frequency is 1.5 MHz and I need the area between 1.4 MHz and 1.6 MHz. I tried to use cumtrapz function using below code:

Int = cumtrapz(f,df);

Intv = @(a,b) max(Int(f<=b)) - min(Int(f>=a));

SegmentArea1 = Intv(1.5-.1, 1.5+.1)

SegmentArea2 = Intv(3-.1, 3+.1)

But it gives the same value for SegmentArea1 and SegmentArea2. I also tried to use trapz using:

trapz(df(find(abs(f-1.4)<0.0001):find(abs(f-1.6)<0.0001)))

but i got a negative value.

Am I using the correct functions or am I doing something wrong? Is there any other way to do it?

Answer 1

[Star Strider]

Calculate the areas of only the ‘positive frequencies’ half of the fft output, not the entire output. If the data are all positive (for example, the absolute value of the fft), negative values are usually because the independent variable (frequencies in an fft) are going from the Nyquist frequency to 0, rather than from 0 to the Nyquist frequency.

[Nasrin Akter]

Can you please take a look at my edited post? I think I didn't post my question corretly.

[Star Strider]

Integrate the absolute value of the fft result without doing any conversions on it. Convert the integrated result to decibels later if you need to, however the reason for doing that is not obvious to me.

Do not integrate the decibel-transformed fft result. That is the likely origin of the negative integration results.

[Nasrin Akter]

If I take the absolute value of the fft, I get just spikes at the desired frequencies. There is no area under the curve to calculate

[Star Strider]

There should be.

Plot the cumtrapz output of the integration of ‘abs(A)’ as a function of the frequency from 0 to 10 MHz. The integrated results should be plainly visible.

[Nasrin Akter]

Yes the result is visible. Which function should I use here? 'cumtrapz' or 'trapz'? Do you think my code is correct with 'trapz'?

[Star Strider]

Use cumtrapz. That gives you cumulative areas for each frequency. Then find a way to automatically segment the parts you are interested in, and calculate the areas from those values.

[Nasrin Akter]

If you don't mind can you help me with the code?

[Star Strider]

Adapt this example code to your data:

t = linspace(0, 5, 1E+4);                                                           % Time Vector
s = sum(sin([100; 250; 300; 700]*2*pi*t));                                          % Signal Vector

figure
plot(t, s)
grid

L = numel(t);
Ts = mean(diff(t));
Fs = 1/Ts;
Fn = Fs/2;
n = nextpow2(numel(t));
FTs = fft(s,2^n)/L;
Fv = linspace(0, 1, fix(2^n/2)+1)*Fn;
Iv = 1:numel(Fv);

figure
plot(Fv, abs(FTs(Iv))*2)
grid

A = cumtrapz(Fv, abs(FTs(Iv))*2);
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0);                                % Find Approximate Zero-Crossing Indices
lims = zci(abs(FTs(Iv))-0.085);                                                     % Integration Limits For Peaks
PkAreas = diff(A(lims));
PkFreqs = Fv(islocalmax(abs(FTs(Iv)), 'MinProminence',0.1));
FA = table(PkFreqs.',PkAreas(1:2:end).', 'VariableNames',{'Frequency','Area'})

[Nasrin Akter]

Thank you. I have one more question. Is it possible to do baseline correction on a FFT signal? If yes, how?

[Star Strider]

My pleasure.

I doubt that there is any need to do a baseline correction of the fft output. If you want to remove the D-C offset (0 Hz peak) from the original time-domain signal (in order to see the other peaks more clearly), subtract the mean of the signal from the signal before doing the fft.

Depending on what the problem is with the fft output, the detrend function may be the best option.