Real time octave analysis: one third octave band

Hi,

I'm trying to use a library called: Real time octave analysis.

Link: https://www.mathworks.com/matlabcentral/fileexchange/3076-real-time-octave-analysis?focused=6777887&tab=function

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. 
%           
%    See also OCT3DSGN, OCT3SPEC, OCTDSGN, OCTSPEC.

% Author: Christophe Couvreur, Faculte Polytechnique de Mons (Belgium)
%         [e-mail]
% Last modification: Aug. 23, 1997, 10:30pm.

% References: 
%    [1] ANSI S1.1-1986 (ASA 65-1986): Specifications for
%        Octave-Band and Fractional-Octave-Band Analog and
%        Digital Filters, 1993.
%    [2] 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 [2]). 
[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

Answer 1

[Gabriele Bunkheila]

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]

Perfect and clear! Thank you again!