Real time octave analysis
MASSIMO MERONI
I read carefully a prevoius article "Real time octave
analysis" posted on the Mathworks website and I dowloaded
and applied the Matlab files oct3dsgn.m and oct3bank.m to
my vibration analysis.
I work on mechanical issues, I know something about
aliasing problems and maybe nothing more than this.
I changed the oct3bank.m file for my vibration analysis in
order to wide the frequencies band to 10 HZ --- 10000HZ (
the original file was with the band 100 HZ --- 5000 HZ).
I ask you politely if you can say I made some big errors
and if my octave analysis may create some fake results.
(all my m. files work correctly).
I thank you in advance, I hope you can help me.
best regards,
Massimo
this is my program:
function [p,f] = oct3bank_my(x);
narg=nargin;
pi = 3.14159265358979;
Fs =44100; % Sampling
Frequency
N = 3; % Order of analysis
filters.
F = [10 12.5 16 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]; % Preferred labeling freq.
ff = (1000).*((2^(1/3)).^[-20:13]); % Exact center
freq.
P = zeros(1,31);
if narg
m = length(x);
end;
% Design filters and compute RMS powers in 1/3-oct. bands
% 10000 Hz band to 1600 Hz band, direct implementation of
filters.
BB=[];
AA=[];
for i =32:-1:22
[B,A] = oct3dsgn(ff(i),Fs,N);
if narg
y = filter(B,A,x);
P(i) = sum(y.^2)/m;
else
BB(i,:)=B;
AA(i,:)=A;
end;
end
if narg
% 1250 Hz to 10 Hz, multirate filter implementation
(see [2]).
[Bu,Au] = oct3dsgn(ff(25),Fs,N); % Upper 1/3-
oct. band in last octave.
[Bc,Ac] = oct3dsgn(ff(24),Fs,N); % Center
1/3-oct. band in last octave.
[Bl,Al] = oct3dsgn(ff(23),Fs,N); % Lower 1/3-
oct. band in last octave.
%
for j = 6:-1:0
x = decimate(x,2);
m = length(x);
y = filter(Bu,Au,x);
P(j*3+4) = sum(y.^2)/m;
y = filter(Bc,Ac,x);
P(j*3+3) = sum(y.^2)/m;
y = filter(Bl,Al,x);
P(j*3+2) = sum(y.^2)/m;
end
x = decimate(x,2);
y = filter(Bu,Au,x);
P(1) = sum(y.^2)/m;
% Convert to decibels.
Pref = 10e-5; % Reference
level for dB scale.
idx = (P>0);
P(idx) = 10*log10(P(idx)/Pref);
P(idx) =1*(P(idx));
P(~idx) = NaN*ones(sum(~idx),1);
% Generate the plot
if (nargout == 0)
bar(P);
ax = axis;
axis([0 33 ax(3) ax(4)])
set(gca,'XTick',[2:3:28]); %
Label frequency axis on octaves.
%set(gca,'XTickLabels',F(2:3:length(F)));
% MATLAB 4.1c
set(gca,'XTickLabel',F(2:3:length(F)));
% MATLAB 5.1
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
else
BBx=BB(13:15,:); %Section to repeat
AAx=AA(13:15,:);
for ii=1:4
ind=(ii-1)*3+1:ii*3;
BB(ind,:)=BBx;
AA(ind,:)=AAx;
end;
p=BB;
f=AA;
end;
MASSIMO MERONI
I read carefully a prevoius article "Real time octave