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?
0 个评论
采纳的回答
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
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
2022-7-18
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
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Octave 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!