How to conduct octave analysis on frequency domain siganal?

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?

Answer 1

[Mathieu NOE]

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

[Csanad Levente Balogh]

Thank you! This is an exelent code. I am still a bit confused about the values of the frequency bands though. After a procedure, I calculate a complex-valued spectrum of a force . Doesn't octave analysis repsesent some kind of energy distribution? So for the spectrum variable in your code, should the power spectrum () be calculated, or is is enought to just average the absolute values of the original spectrum, or maybe something different from both? I actually found a lot of tutorials on how to calculate the central frequences and the frequency limits, but not a lot on how to get the rigth values for the octave bands.

[Mathieu NOE]

Octave analysis have its roots in early analog signal processing . Before the modern era of digital processing, signal could be decomposed into octave and later 1/3 octave bands by using analog bandpass filters followed by a form of "rms" computation circuitry . So you could look at RMS values for some frequency bands

Why not follow the same protocol and ensure that each 1/3 octave correspond to the rms value of the signal / spectrum in that frequency range.

It is then necessary to compute a PSD of your data , convert to power by integrating this PSD over the frequency between each pair of low / upper frequency bounds and take the square root of that to get the rms value.

Power Spectral Density Estimates Using FFT - MATLAB & Simulink (mathworks.com)

[Csanad Levente Balogh]

Thank you!

[Mathieu NOE]

My pleasure !

[Bryan Wilson]

I'm confused. Isn't the octave and 1/3 octave levels the SUM of the signal amplitudes in linear units within the frequency band, rather than the average?

[Mathieu NOE]

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