How to conduct octave analysis on frequency domain siganal?

123 views (last 30 days)
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?

Accepted Answer

Mathieu NOE
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
Bryan Wilson on 18 Jul 2022
Edited: Bryan Wilson on 18 Jul 2022
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
Mathieu NOE on 18 Jul 2022
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

Sign in to comment.

More Answers (0)

Categories

Find more on Measurements and Spatial Audio in Help Center and File Exchange

Products


Release

R2019b

Community Treasure Hunt

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

Start Hunting!