# How to conduct octave analysis on frequency domain siganal?

123 views (last 30 days)

Show older comments

Csanad Levente Balogh
on 8 Mar 2021

Commented: Mathieu NOE
on 18 Jul 2022

Hi guys,

Maybe it is a silly question, but I want to calculate 1/3 octave analysis on an already frequency domain signal stored as a vector in MATLAB. I know how to calculate the center frequencies and the upper and lower frequuency bounds. What I'm not sure about is the value of the octave bands. Should I just average the values of the spectrum of the signal between lower and upper bounds? How does this work?

##### 0 Comments

### Accepted Answer

Mathieu NOE
on 8 Mar 2021

hello

this is a code that does what you are looking for

now this simply do a linear average of the spectrum amplitude between each lower / upper frequency bound

if your spectrum is given in dB , you have to convert first back to linear scale and do the conversion,then do the dB conversion (with averaging) of the 1/3 octave spectrum

% conversion fft narrow band spectrum to 1/3 octave

%% dummy data

freq = linspace(100,1000,100);

spectrum = 25+5*randn(size(freq));

[fto,sTO] = conversion2TO(freq,spectrum);

figure(1), plot(freq,spectrum,'b',fto,sTO,'*-r');

legend('narrow band FFT spectrum','1/3 octave band spectrum');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [fto,sTO] = conversion2TO(freq,spectrum)

% normalized 1/3 octave center freqs

fref = [10, 12.5 16 20, 25 31.5 40, 50 63 80, 100 125 160, 200 250 315, ...

400 500 630, 800 1000 1250, 1600 2000 2500, 3150 4000 5000, ...

6300 8000 10000, 12500 16000 20000 ];

ff = (1000).*((2^(1/3)).^[-20:13]); % Exact center freq.

a = sqrt(2^(1/3)); %

f_lower_bound = ff./a;

f_higher_bound = ff.*a;

ind1 = find (f_higher_bound>min(freq)); ind1 = ind1(1); % indice of first value of "f_higher_bound" above "min(freq)"

ind2 = find (f_lower_bound<max(freq)); ind2 = ind2(end); % indice of last value "f_lower_bound" below "max(freq)"

ind3 = (ind1:ind2);

for ci = 1:length(ind3)

ind4 = find(freq>=f_lower_bound(ind3(ci)) & freq<=f_higher_bound(ind3(ci)));

sTO(ci) = mean(spectrum(ind4)); % 1/3 octave value = averaged value of spectrum inside 1/3 octave band

fto(ci) = fref(ind3(ci)); % valid central frequency 1/3 octave

end

end

##### 6 Comments

Bryan Wilson
on 18 Jul 2022

Edited: Bryan Wilson
on 18 Jul 2022

Mathieu NOE
on 18 Jul 2022

hello @Bryan Wilson

you are 100% right !! my bad !

this is the corrected code :

% conversion fft narrow band spectrum to 1/3 octave

clc

clearvars

%% dummy data

freq = logspace(1,4,100);

% FFT narrow band spectrum (in dB !!)

spectrum_dB = 45+5*randn(size(freq));

spectrum_dB(34) = 100;

spectrum_dB(44) = 90;

spectrum_dB(54) = 80;

[fTO,sTO_dB] = conversion2TO(freq,spectrum_dB);

figure(1),

semilogx(freq,spectrum_dB,'b',fTO,sTO_dB,'*-r');

xlabel('Freq (Hz)');

ylabel('Amplitude (dB)');

legend('narrow band FFT spectrum','1/3 octave band spectrum');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [fTO,sTO_dB] = conversion2TO(freq,spectrum_dB)

% convert back from dB to linear amplitude

spectrum = 10.^(spectrum_dB/20);

% normalized 1/3 octave center freqs

fref = [10, 12.5 16 20, 25 31.5 40, 50 63 80, 100 125 160, 200 250 315, ...

400 500 630, 800 1000 1250, 1600 2000 2500, 3150 4000 5000, ...

6300 8000 10000, 12500 16000 20000 ];

ff = (1000).*((2^(1/3)).^[-20:13]); % Exact center freq.

a = sqrt(2^(1/3)); %

f_lower_bound = ff./a;

f_higher_bound = ff.*a;

ind1 = find (f_higher_bound>min(freq)); ind1 = ind1(1); % indice of first value of "f_higher_bound" above "min(freq)"

ind2 = find (f_lower_bound<max(freq)); ind2 = ind2(end); % indice of last value "f_lower_bound" below "max(freq)"

ind3 = (ind1:ind2);

for ci = 1:length(ind3)

ind4 = (freq>=f_lower_bound(ind3(ci)) & freq<=f_higher_bound(ind3(ci)));

sTO_dB(ci) = 10*log10(sum(spectrum(ind4).^2)); % 1/3 octave value = RMS sum of spectrum inside 1/3 octave band

fTO(ci) = fref(ind3(ci)); % valid central frequency 1/3 octave

end

end

### More Answers (0)

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!