求助:帮忙解读程序！
我是一个初学者，在理解这个程序时遇到了一些麻烦，烦请各位帮帮忙，帮我讲解一下这个程序的算法！谢谢各位！function = filtbank(x,Fs,T,arg4);
% FILTBANK One-third-octave band frequency analyser.
% = FILTBANK(X,Fs) computes the one-third-octave spectrum
% of X, assuming sampling frequency Fs (in Hz). P is a row vector
% containing RMS powers computed for each 1/3-octave band and expressed
% in dB with 1 as reference level. F contains the corresponding preferred
% labeling frequencies (standard ANSI S1.6-1984).
%
% Plots can be obtained, e.g., with BANKDISP.
% Example: (spectrum of white noise)
% X = rand(1,100000);
% = filtbank(X,44100);
% bankdisp(P,F,-40,-20);
%
% The frequency range covered is the standard 'restricted' audio range:
% from 100 Hz to 5000 Hz.
%
% = FILTBANK(X,Fs,T) computes a sequence of one-third-octave
% spectra with integration time T (in s). P is a matrix of size
% (LENGTH(X)/(T*Fs)) x LENGTH(F). If T = [] (default), the integration
% time is equal to the length of X.
%
% = FILTBANK(X,Fs,T,'extended') covers the 'extended' audio range:
% from 25 Hz to 20 000 Hz. Note that depending on the value of Fs
% the frequency range might be automatically reduced. Warnings will be
% issued if this is the case.
%
% See also BANKDISP, LEQ.
if (nargin < 2) | (nargin > 4)
error('Invalid number of arguments');
end
pi = 3.14159265358979;
Fref = [ 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 ]; % Preferred labeling freq.
Fc = 1000*((2^(1/3)).^[-16:1:13]); % Exact center freq.
N = 3; % Order of analysis filters.
extended = 0;
U = 2^(1/3);
% Integration time.
if (nargin>=3)
if isempty(T)
T = length(x);
P = zeros(1,length(Fref));
else % Convert T to number of samples.
T = floor(Fs*T);
P = zeros(floor(length(x)/T),length(Fref));
end
else
T = length(x);
P = zeros(1,length(Fref));
end
if (nargin >= 4)
if strcmp(lower(arg4),'extended')% Extended (25 Hz to 20000 Hz)
extended = 1;
end
end
% Frequency range.
if (extended) % extended (25 Hz to 20 000 Hz).
i_up = 30;
i_low = 1;
else % restricted (100 Hz to 5000 Hz).
i_up = 24;
i_low = 7;
end
% Check sampling frequencies and issue warnings.
if (Fs/2) < Fref(i_low)*1.5
error('Sampling frequency Fs too low.');
elseif (Fs/2) < Fref(i_up)*1.5
disp('Warning: frequency range must be reduced (Fs too low).');
i_up = max(find(Fc<=Fs/3));
end
% Compute 'pivot' frequency for multirate filter implementation
% All filters below Fs/20 will be implemented after a decimation.
if (Fc(i_low) > Fs/20)
i_dec = 0;
else
i_dec = max(find(Fc<=Fs/20));
end
% Design filters and compute RMS powers in 1/3-oct. bands.
% Higher frequencies, direct implementation of filters.
for i = i_up:-1:i_dec+1
= oct3dsgn(Fc(i),Fs,N);
y = filter(B,A,x);
P(:,i) = leq(y,T);
end
% Lower frequencies, decimation by series of 3 bands.
if (i_dec > 0)
= oct3dsgn(Fc(i_dec),Fs/2,N); % Upper 1/3-oct. band in last octave.
= oct3dsgn(Fc(i_dec)/U,Fs/2,N); % Center 1/3-oct. band in last octave.
= oct3dsgn(Fc(i_dec)/(U^2),Fs/2,N); % Lower 1/3-oct. band in last octave.
i = i_dec;
whilei >= i_low+2
x = decimate(x,2);
T = T/2;
y = filter(Bu,Au,x);
P(:,i) =leq(y,T);
y = filter(Bc,Ac,x);
P(:,i-1) =leq(y,T);
y = filter(Bl,Al,x);
P(:,i-2) =leq(y,T);
i = i-3;
end
if (i == (i_low+1))
x = decimate(x,2);
T = T/2;
y = filter(Bu,Au,x);
P(:,i) =leq(y,T);
y = filter(Bc,Ac,x);
P(:,i-1) =leq(y,T);
elseif (i == (i_low))
x = decimate(x,2);
T = T/2;
y = filter(Bu,Au,x);
P(:,i) =leq(y,T);
end
end
f = Fref(i_low:i_up);
p = P(:,i_low:i_up);
页:
[1]