[求助]请教根据声压信号求A计权声压值的MATLAB程序

wangwenting

我自己写了一个，但是求出来跟声级计得到的信号差距还蛮大的
有没有哪位能给我提供一个根据声压信号求A计权声压值的MATLAB程序啊
不胜感激！！！

wangwenting

贴个我自己写的程序，根据oct3bank写的，但是算出来不太对
    [p,ff]=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 [p,f] = 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
   [B,A] = 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 [2]).
[Bu,Au] = oct3dsgn(ff(24),Fs,N);         % Upper 1/3-oct. band in last octave.
[Bc,Ac] = oct3dsgn(ff(23),Fs,N);         % Center 1/3-oct. band in last octave.
[Bl,Al] = 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([0 28 ax(3) ax(4)])
  set(gca,'XTick',[2:3:27]);                 % 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

function [B,A] = 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 [1].
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;
[B,A] = butter(N,[W1,W2]);