wangwenting 发表于 2009-3-16 10:03

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

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

yelv123 发表于 2009-3-16 10:41

自己可以把程序贴出来,让大家看看,可以帮你找出错误。。对于求程序

wangwenting 发表于 2009-3-18 16:21

回复 沙发 yelv123 的帖子

程序如下,请教各位帮忙看看问题在哪儿
谢啦!

    =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]
查看完整版本: [求助]请教根据声压信号求A计权声压值的MATLAB程序