Real time octave analysis: one third octave band

82 次查看（过去 30 天）
Nycholas Maia
Nycholas Maia 2018-1-4
评论： Nycholas Maia 2018-1-8
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)
% 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
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.

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Measurements and Spatial Audio 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by