[求助]请教根据声压信号求A计权声压值的MATLAB程序
我自己写了一个，但是求出来跟声级计得到的信号差距还蛮大的有没有哪位能给我提供一个根据声压信号求A计权声压值的MATLAB程序啊
不胜感激！！！
回复 楼主 wangwenting 的帖子
贴个我自己写的程序，根据oct3bank写的，但是算出来不太对=oct3bank(noi); %noi是从声级计中采出来的声压信号
plot(ff,p);
Cj=[-19.1 -16.1 -13.4 -10.9 -8.6 -6.6 -4.8 -3.2 -1.9 -0.8 0 0.6 1.0 1.2 1.3 1.2 1.0 0.5 -0.1 -1.1 -2.5];
Lj=zeros(1,21);
for n=1:21
Lj(n)=10^(0.1*(p(n)+Cj(n)));
end
Lp(i)=10*log10(sum(Lj));
function = oct3bank(x);
pi = 3.14159265358979;
Fs = 44100; % Sampling Frequency
N = 3; % Order of analysis filters.
F = [ 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)).^[-16:10]); % Exact center freq.
P = zeros(1,27);
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.
for i = 27:-1:22
= oct3dsgn(ff(i),Fs,N);
y = filter(B,A,x);
P(i) = sum(y.^2)/m;
end
% 1250 Hz to 100 Hz, multirate filter implementation (see ).
= oct3dsgn(ff(24),Fs,N); % Upper 1/3-oct. band in last octave.
= oct3dsgn(ff(23),Fs,N); % Center 1/3-oct. band in last octave.
= oct3dsgn(ff(22),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+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 = 0.00002; % Reference level for dB scale.
idx = (P>0);
P(idx) = 20*log10(P(idx)/Pref);
P(~idx) = NaN*ones(sum(~idx),1);
% Generate the plot
if (nargout == 0)
bar(P);
ax = axis;
axis()
set(gca,'XTick',); % 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 '); ylabel('Power ');
title('One-third-octave spectrum')
% Set up output parameters
elseif (nargout == 1)
p = P;
elseif (nargout == 2)
p = P;
f = F;
end
function = oct3dsgn(Fc,Fs,N);
if (nargin > 3) | (nargin < 2)
error('Invalide number of arguments.');
end
if (nargin == 2)
N = 3;
end
if (Fc > 0.88*(Fs/2))
error('Design not possible. Check frequencies.');
end
% Design Butterworth 2Nth-order one-third-octave filter
% Note: BUTTER is based on a bilinear transformation, as suggested in .
pi = 3.14159265358979;
f1 = Fc/(2^(1/6));
f2 = Fc*(2^(1/6));
Qr = Fc/(f2-f1);
Qd = (pi/2/N)/(sin(pi/2/N))*Qr;
alpha = (1 + sqrt(1+4*Qd^2))/2/Qd;
W1 = Fc/(Fs/2)/alpha;
W2 = Fc/(Fs/2)*alpha;
= butter(N,);
页:
[1]