How to conduct octave analysis on frequency domain siganal?
79 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
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
on 18 Jul 2022
hello @Bryan Wilson
you are 100% right !! my bad !
this is the corrected code :
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1069470/image.png)
% 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
Find more on Measurements and Spatial Audio in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!