148 views (last 30 days)

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.

%

% See also OCT3DSGN, OCT3SPEC, OCTDSGN, OCTSPEC.

% Author: Christophe Couvreur, Faculte Polytechnique de Mons (Belgium)

% couvreur@thor.fpms.ac.be

% 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

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.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.