Real time octave analysis: one third octave band

67 views (last 30 days)
Nycholas Maia on 4 Jan 2018
Commented: Nycholas Maia on 8 Jan 2018
Hi,
I'm trying to use a library called: Real time octave analysis.
The original source code do your calculations using a frequency vector F that have only 18 frequency numbers inside it.
Now I want to expand this code to make it possible to do your calculations with a new vector F that have 29 frequency numbers inside it.
This 29 frequency numbers correspond to 20Hz until 20.000Hz using one third octave band.
I did some little modifications inside the function oct3bank(x), but there are some lines of code that I can't understand very well and I need some help to adapt theses lines to my propouse.
Could you help me? I think that my problems will be in the lines:
Line 51 - % MODIFIED FOR LOOP VALUES: OK?
Line 58 - % HOW TO CHANGE THIS VALUES [15 14 13] TO MATCH WITH THE NEW VECTOR F?
Line 80 - % HOW TO PLOT ALL 29 NEW BANDS?
Thanks, Nycholas Maia
function [p,f] = oct3bank(x)
% OCT3BANK Simple one-third-octave filter bank.
% OCT3BANK(X) plots one-third-octave power spectra of signal vector X.
% Implementation based on ANSI S1.11-1986 Order-3 filters.
% Sampling frequency Fs = 44100 Hz. Restricted one-third-octave-band
% range (from 100 Hz to 5000 Hz). RMS power is computed in each band
% and expressed in dB with 1 as reference level.
%
% [P,F] = OCT3BANK(X) returns two length-18 row-vectors with
% the RMS power (in dB) in P and the corresponding preferred labeling
% frequencies (ANSI S1.6-1984) in F.
%
% Author: Christophe Couvreur, Faculte Polytechnique de Mons (Belgium)
% couvreur@thor.fpms.ac.be
% Last modification: Aug. 23, 1997, 10:30pm.
% References:
%  ANSI S1.1-1986 (ASA 65-1986): Specifications for
% Octave-Band and Fractional-Octave-Band Analog and
% Digital Filters, 1993.
%  S. J. Orfanidis, Introduction to Signal Processing,
% Prentice Hall, Englewood Cliffs, 1996.
Fs = 44100; % Sampling Frequency
N = 3; % Order of analysis filters.
% Original F Vector:
%F = [ 100 125 160, 200 250 315, 400 500 630, 800 1000 1250,1600 2000 2500, 3150 4000 5000 ]; % Preferred labeling freq.
% MODIFIED F VECTOR: OK!
F = [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];
% Original ff vector:
%ff = (1000).*((2^(1/3)).^[-10:7]); % Exact center freq.
% MODIFIED FF VECTOR: OK!
ff = (1000).*((2^(1/3)).^(-17:13));
P = zeros(1,length(F));
m = length(x);
% Design filters and compute RMS powers in 1/3-oct. bands
% 5000 Hz band to 1600 Hz band, direct implementation of filters.
% Original for loop values:
%for i = 18:-1:13
% MODIFIED FOR LOOP VALUES: OK?
for i = 25:-1:20
[B,A] = oct3dsgn(ff(i),Fs,N);
y = filter(B,A,x);
P(i) = sum(y.^2)/m;
end
% HOW TO CHANGE THIS VALUES [15 14 13] TO MATCH WITH THE NEW VECTOR F?
% 1250 Hz to 100 Hz, multirate filter implementation (see ).
[Bu,Au] = oct3dsgn(ff(15),Fs,N); % Upper 1/3-oct. band in last octave.
[Bc,Ac] = oct3dsgn(ff(14),Fs,N); % Center 1/3-oct. band in last octave.
[Bl,Al] = oct3dsgn(ff(13),Fs,N); % Lower 1/3-oct. band in last octave.
for j = 3:-1:0
x = decimate(x,2);
m = length(x);
y = filter(Bu,Au,x);
P(j*3+3) = sum(y.^2)/m;
y = filter(Bc,Ac,x);
P(j*3+2) = sum(y.^2)/m;
y = filter(Bl,Al,x);
P(j*3+1) = sum(y.^2)/m;
end
% Convert to decibels.
Pref = 1; % Reference level for dB scale.
idx = (P>0);
P(idx) = 10*log10(P(idx)/Pref);
P(~idx) = NaN*ones(sum(~idx),1);
% HOW TO PLOT ALL 29 NEW BANDS?
% Generate the plot
if (nargout == 0)
bar(P);
ax = axis;
axis([0 19 ax(3) ax(4)])
set(gca,'XTick',(2:3:length(F))); % Label frequency axis on octaves.
set(gca,'XTickLabel',F(2:3:length(F)));
xlabel('Frequency band [Hz]'); ylabel('Power [dB]');
title('One-third-octave spectrum')
% Set up output parameters
elseif (nargout == 1)
p = P;
elseif (nargout == 2)
p = P;
f = F;
end

Gabriele Bunkheila on 8 Jan 2018
Hi Nycholas,
I won't comment on the specifics of this function, which was written quite a while ago, but I thought I'd provide a pointer to the reference documentation page of the octaveFilter object in case it could help you or others.
Within that page you'll find a reference to the method getANSICenterFrequencies, e.g. to be used as follows to get the desired vector of frequencies that you mention:
allANSIFrequencies = getANSICenterFrequencies(octaveFilter('Bandwidth','1/3 octave'));
freqSubset = allANSIFrequencies(8:end-1);
Further down on that same page you'll find other relevant code snippets.
For example, to design an array of 1/3 octave filters centered around those frequencies:
for i = 1:length(freqSubset)
octaveFilterBank{i} = octaveFilter(freqSubset(i),'FilterOrder',12);
end
More code is provided in the same page for plotting the whole filterbank (see also this other question of yours on MATLAB Answers) or for estimating the energy of a given signal in the desired bands.
Thanks,
Gabriele.
Nycholas Maia on 8 Jan 2018
Perfect and clear! Thank you again!