How to conduct octave analysis on frequency domain siganal?

114 次查看(过去 30 天)
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?

采纳的回答

Mathieu NOE
Mathieu NOE 2021-3-8
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 个评论
Bryan Wilson
Bryan Wilson 2022-7-18
编辑:Bryan Wilson 2022-7-18
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 2022-7-18
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

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Octave 的更多信息

产品


版本

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by