Zwicker响度-Methods for calculating loudness

Too_L

已于 2024-08-09 14:07:42 修改

阅读量721

点赞数 18

文章标签：网络

于 2024-08-08 21:42:03 首次发布

本文链接：https://blog.csdn.net/qq_52501378/article/details/141035366

版权

写在前面：最近做了一些和Zwicker响度方面有关的调研，此文主要记录个人的一些日常总结，和分享相关的资料，如果只需要最终结论可以直接跳到Matlab实现中的第三节，有相关代码！~

ISO532系列标准的目的是指定基于声音物理特性的计算程序,以估计在特定聆听条件下耳科正常听力的人所感知的声音的响度和响度级别。

ISO532-1和ISO532-2规定了两种不同的响度计算方法,对于给定的声音可能会产生不同的结果。由于目前无法说明对一种或另一种方法的一般偏好,因此由用户选择最适合给定情况的方法。下面介绍每种方法的一些主要特征,以方便选择。

ISO532中的中文文档可见：ISO532-1-2017响度（译文） - 道客巴巴

Critical band 临界频带

在声学中，“Bark” 是一个用于描述听觉频率的单位。它是由德国心理声学家 Eberhard Zwicker 提出的，表示听觉频率的心理尺度。Bark 频率尺度将频率划分为若干个等感知带宽的区间，每个区间称为一个 Bark。这个尺度试图更准确地反映人类听觉系统对不同频率的感知。

具体来说，Bark 频率尺度将听觉频率从 20 Hz 到 20,000 Hz 划分为 24 个 Bark 单位。每个 Bark 单位对应一个等效矩形带宽（ERB），这意味着在这个频率范围内，人耳对声音的感知是相对均匀的。

示例：中心频率为1000Hz的临界带宽值约为160Hz,因此频率从920Hz增加到1080Hz对应于1Bark的步长。

通常将人耳可听范围内的20Hz～16kHz分成24个临界频带，用临界频带级来表示临界频带的宽度，单位为巴克（Bark）： Bark=一个临界频带的宽度。当频率f＜500Hz时，1 Bark=f/100，临界带宽几乎恒定为100Hz；当频率f＞500Hz时，1 Bark=4log(f/100)，临界带宽随中心频率的升高而增加，约为中心频率的20%。不同中心频率的临界带宽、频率边界、临界频带级，如下表所示。

相关资料可见：【听力学名词释义（15）】临界带宽临界频带

响度和响度级的计算

表示响度的物理量为宋（Sone）和方（Phon），宋和方的关系类似于声压和声压级的关系。大体上，声音(1000Hz纯音)增加10dB(声压级)，其响度约增加1倍。

宋(sone)是响度衡量(loudness scaling)实验的结果,Stevens于1936年提出用宋作为另一个衡量响度的单位。宋是由响度感受主观度量试验而来。同样的,事先须规定一个标声作为参考。约定俗成的标准是以40方的响度(40dBSPL,1kHz)为1宋,比这个标准声听起来响度大N倍的声音的响度即为N宋,如果听起来响度只只有这个标准声响度的一半,那么其响度值即为0.5宋,依此类推。实验数据表明,纯音的响度每增加10方,响度听起来会加倍,因此响度宋与方之间的换算关系可用以下的函数式来表达z:sone=2^[(phon-40)/10]。根据上述关系式，40方、60方和80方的响度分别对应1宋、4宋和16宋，实际上，上述方与宋的函数式只对40方以上的响度才成立，即响度每增加10方，以宋为单位的响度值即加倍。当响度低于40方时，宋与方之间的换算关系略有改变，此时响度增加10方，宋值增加将超过1倍。

$N$ ：以宋(sone)为单位的响度

$L_N$ ：以方(phon)为单位的响度级

下面给出它们换算的方法：

$N=2^{0.1(\frac{L_{N}}{phon}-40)}~sone, ~L_N\geq 40~phon$

$L_N= \left\{ \begin{aligned} &[40+33.22lg\frac{N}{sone}]~phon,~N\geq1~ sone\\&40{(\frac{N}{sone}+0.0005)}^{0.35} ~phon,~N<1~sone \end{aligned} \right.$

根据上述关系式，40phon、60phon和80phon的响度分别对应1宋、4宋和16宋。

人耳的等响曲线

响度级→声压级

IS0226 —2003标准H等响度曲线的定义如下,

首先给出如下定义，假设有频率为 $f$ 的纯音，它的声响度级是 $L_N$ ；

接下来我们给出三个与频率 $f$ 有关的函数，其中 $T_f$ 为听力阈值； $a_f$ 为响度感知指数； $L_u$ 为以1000Hz为标准所计算的线性传输函数的幅值。其具体值在代码中给出

接下来我们首先算出一个新的函数 $A_f$ :

$A_f=4.4710*10^{\left( -3 \right)}\left( 10^{\left( 0.025L_N \right)}-1.15 \right) +\left[ 0.4*10^{\left( \frac{T_f+L_U}{10}-9 \right)} \right] ^{a_f}$

则当前响度对应不同频率 $f$ 的声压级 $L_p$ 为:

$L_p=(\frac{10}{a_f}*lgA_f-L_u+94)~\mathrm{dB}$

即： $L_p=f(f,L_N)$

给出相关代码：

%响度计算函数
function [spl, freq] = iso226(phon)
%                /---------------------------------------\
%%%%%%%%%%%%%%%%%          TABLES FROM ISO226             %%%%%%%%%%%%%%%%%
%                \---------------------------------------/
f = [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 12500];

af = [0.532 0.506 0.480 0.455 0.432 0.409 0.387 0.367 0.349 0.330 0.315 ...
      0.301 0.288 0.276 0.267 0.259 0.253 0.250 0.246 0.244 0.243 0.243 ...
      0.243 0.242 0.242 0.245 0.254 0.271 0.301];

Lu = [-31.6 -27.2 -23.0 -19.1 -15.9 -13.0 -10.3 -8.1 -6.2 -4.5 -3.1 ...
       -2.0  -1.1  -0.4   0.0   0.3   0.5   0.0 -2.7 -4.1 -1.0  1.7 ...
        2.5   1.2  -2.1  -7.1 -11.2 -10.7  -3.1];

Tf = [ 78.5  68.7  59.5  51.1  44.0  37.5  31.5  26.5  22.1  17.9  14.4 ...
       11.4   8.6   6.2   4.4   3.0   2.2   2.4   3.5   1.7  -1.3  -4.2 ...
       -6.0  -5.4  -1.5   6.0  12.6  13.9  12.3];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    

%错误处理
if((phon < 0) || (phon > 90))
    disp('Phon value out of bounds!')
    spl = 0;
    freq = 0;
else
    Ln = phon;
    %从响度级计算声压级 
    Af=4.47E-3 * (10.^(0.025*Ln) - 1.15) + (0.4*10.^(((Tf+Lu)/10)-9 )).^af;
    Lp=((10./af).*log10(Af)) - Lu + 94;
    spl = Lp;  
    freq = f;
end

clc
clear all
phon=0:10:90;
for phon=phon
  [spl,freq]=iso226(phon);          %计算声压级
  figure(1)
  semilogx(freq,spl,'-','color','k')
  axis([20,20000,-10,130])
  title('Phon')
  xlabel('频率(Hz)')
  ylabel('声压级别(dB)')
  set(gca,'ytick',-10:10:130)
  hold on
  grid on
  box off
end
line([20,20],[-10,140])

运行效果：

声压级→响度级

现在已经有了声压级 $L_p$ ，希望得到它的声响度级 $L_N$ ；将上述函数取逆

$A_f=10^{\frac{a_f}{10}(L_p+L_u-94)}$

$L_N=40*lg\left[\left[A_f-\left[ 0.4*10^{\left( \frac{T_f+L_U}{10}-9 \right)} \right] ^{a_f}\right]*10^3/4.471+1.15\right]$

为了简化代码，定义一个 $\beta=0.4*10^{\frac{T_f+L_U}{10}-9}$ ，则：

$L_N=40*lg\left[1.15+10^3/4.471*\left[A_f-\beta^{a_f}\right]\right]$

即 $L_N=f(f,L_p)$

代码如下：

clc;clear;
Source=zeros(1,1);
Source=load("C:\Users\anker\Downloads\Bose TNC.txt");

%对原数据插值
x1=logspace(1,5,120)
% x1=[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 12500];
y1=interp1(Source(:,1),Source(:,2),x1,'linear');
y1_e=y1+90;%加90dB后的频响

%初始噪声，全频段90dB
%对三个函数都进行插值
f = [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 12500];
af = [0.532 0.506 0.480 0.455 0.432 0.409 0.387 0.367 0.349 0.330 0.315 ...
      0.301 0.288 0.276 0.267 0.259 0.253 0.250 0.246 0.244 0.243 0.243 ...
      0.243 0.242 0.242 0.245 0.254 0.271 0.301];
af1=interp1(f,af,x1,'linear');

Lu = [-31.6 -27.2 -23.0 -19.1 -15.9 -13.0 -10.3 -8.1 -6.2 -4.5 -3.1 ...
       -2.0  -1.1  -0.4   0.0   0.3   0.5   0.0 -2.7 -4.1 -1.0  1.7 ...
        2.5   1.2  -2.1  -7.1 -11.2 -10.7  -3.1];
Lu1=interp1(f,Lu,x1,'linear');

Tf = [ 78.5  68.7  59.5  51.1  44.0  37.5  31.5  26.5  22.1  17.9  14.4 ...
       11.4   8.6   6.2   4.4   3.0   2.2   2.4   3.5   1.7  -1.3  -4.2 ...
       -6.0  -5.4  -1.5   6.0  12.6  13.9  12.3];
Tf1=interp1(f,Tf,x1,'linear');

%%%%%%%%%%%%%%%%%%%白噪声90dB%%%%%%%%%%%%%

x2=x1;
y2=90.*ones(1,120);

%声压级换算响度级
A_f=zeros(1,1);
L_N=zeros(1,1);
beta=zeros(1,1);

for i =1:length(x2)
    A_f(i)=10^(af1(i)/10.*(Lu1(i)+y2(i)-94));
    beta(i)=0.4.*10^((Tf1(i)+Lu1(i))./10-9);
    L_N2(i)=40*log10(1.15+1e3./4.471.*(A_f(i)-beta(i).^(af1(i))));
end


figure(1)
yyaxis left
semilogx(x2,y2,'LineStyle','-','LineWidth',2)
axis([60,16300,30,100])
ylabel('SPL(dB)')


yyaxis right
semilogx(x2,L_N2,'LineStyle',':','LineWidth',2)
axis([60,16300,30,100])

legend({'声压级','响度级'},'Location','south')
ylabel('方(phon)')
xlabel('Frequency (Hz)')

结果

播放90dB SPL的白噪声后

1.环境中的声压级与响度级，3kHz附近实际听感响度更强

2.配到降噪耳机后测得的声压级+响度级，3kHz附近实际听感响度更强

3.降噪深度的声压级与响度级，<1000 Hz时的降噪响度级更加显著

4.为了简化计算，尝试用平均SPL代替每个频点的SPL，替代后曲线如下，两者差距不明显。

Zwicker方法计算稳态声音：

原文解读

为了看懂上面这张图，给出新文档：Zwicker Loudness Measurement

首先是一个噪声频带：

在 Zwicker 的模型中，考虑到人耳频带的划分，对于300Hz以下的部分，我们通过如下方式进行加权：

合并原则如下：中心频率在 90Hz 以下的声压级相加为 L1；中心频率在 90-180Hz 的声压级相加为 L2；中心频率在 180-280Hz 的声压级相加为 L3。如下为例子：

给出如下的响度计算模板 (DIN, 1991-03) ，每个频段的声压级根据相应频段的标度用水平线蓝线输入。

每条蓝线的高度代表该频段的响度。上升的斜坡是垂直的。下降坡度与虚线"平行"输入。

蓝色曲线下方区域的大小代表总响度。可以通过找到一条平均蓝线的直线来估计区域的大小,即蓝线下方的面积约等于直线下方的面积。水平线与垂直刻度交叉处可以读取总响度。在本示例中，计算出的响度约为2.5 sone，响度为53.5phon。

然后是下面这张图

首先给出它纵坐标 $N'$ 的定义，specific loudness ，在crtitical band的宽带中心引起的响度,ISO标准后面给出了其计算公式：

$N_C'=0.0635\times 10^{0.025L_{TQ}}((1-s+s\times 10^{0.1(L_{CB}-L_{TQ})})^{0.25}-1)$

其中的 $L_{CB}$ 是在critical band内得到的声压级SPL， $L_{TQ}$ 则是critical band的中心频率在安静环境下的可听阈声压级SPL，其对应的值可以在ISO后面表中找到。这里给的图片介绍是：specific loudness versus critical band rate pattern for a single tone with the frequency f=1 kHz。但是这里实际上使用了三分之一倍频程滤波器，可以在相邻频段提供20dB的阻尼，因此在800Hz和1250Hz也有50dB的SPL输入。

Critical band level

$L_{CB}$

sound pressure level of sound contained within a critical band

specific loudness

$N'$

loudness (3.18) evoked over a frequency band with a bandwidth of a critical band (3.12) centred at the frequency of interest

Note 1 to entry:Specific loudnessis expressed in sones/Bark.

Note 2 to entry:The definition together with the stated unit are different from those in ISO 532-2.

各个cirtical band 中的安静听阈

接下来我们要验证上面提到的Specific loudness公式的准确性，这里首先给了一个50dB的全频白噪声，看到它在800Hz处的Specific loudness是0.72，对应第八个Bark。然后给一个70dB白噪声，看到它对应的Specific loudness是2.45，对应第九个Bark。两者都与ISO上给出的图片对应，而第十个Bark的Specific loudness也应处于0.7左右，因为差值过大，可以看到它都没有出现，被下降坡度所覆盖。这里Specific loudness的面积对应的就是这段信号的总响度，为8宋(soon),即70方(phon)。

Matlab实现

1.三分之一倍频程滤波器

首先，文中提到了滤波器的使用。滤波器针对的是时域信号，这里找到了一种巴特沃斯滤波器butter，能够实现ISO标准中提到了在相邻1/3倍频程内20dB的衰减，具体见如下代码：

function [p,f] = oct3bank(x);
% OCT3BANK Simple one-third-octave filter bank.
% OCT3BANK(X) plots one-third-octave power spectra of signal vector X.
% Implementation based on ANSI S1.11-1986 Order-3 filters.
% Sampling frequency Fs = 44100 Hz. Restricted one-third-octave-band
% range (from 100 Hz to 5000 Hz). RMS power is computed in each band
% and expressed in dB with 1 as reference level.
%
% [P,F] = OCT3BANK(X) returns two length-18 row-vectors with
% the RMS power (in dB) in P and the corresponding preferred labeling
% frequencies (ANSI S1.6-1984) in F.
% Author: Christophe Couvreur, Faculte Polytechnique de Mons (Belgium)
% couvreur@thor.fpms.ac.be
% Last modification: Aug. 23, 1997, 10:30pm.
Fs = 44100; % Sampling Frequency
N = 3; % Order of analysis filters.
F = [ 100 125 160, 200 250 315, 400 500 630, 800 1000 1250, ...
    1600 2000 2500, 3150 4000 5000 ]; % Preferred labeling freq.
ff = (1000).*((2^(1/3)).^[-10:7]); % Exact center freq.
P = zeros(1,18);
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 = 18:-1:13
    [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.
[Bu,Au] = oct3dsgn(ff(15),Fs,N); % Upper 1/3-oct. band in last octave.
[Bc,Ac] = oct3dsgn(ff(14),Fs,N); % Center 1/3-oct. band in last octave.
[Bl,Al] = oct3dsgn(ff(13),Fs,N); % Lower 1/3-oct. band in last octave.
for j = 3:-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 = 1; % Reference level for dB scale.
idx = (P>0);
P(idx) = 10*log10(P(idx)/Pref);
P(~idx) = NaN*ones(sum(~idx),1);
% Generate the plot
if (nargout == 0)
    bar(P);
    ax = axis;
    axis([0 19 ax(3) ax(4)])
    set(gca,'XTick',[2:3:18]); % Label frequency axis on octaves.
    % set(gca,'XTickLabels',F(2:3:length(F))); % Syntax for MATLAB 4.1c
    set(gca,'XTickLabel',F(2:3:length(F))); % Syntax for 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);
% OCT3DSGN Design of a one-third-octave filter.
% [B,A] = OCT3DSGN(Fc,Fs,N) designs a digital 1/3-octave filter with
% center frequency Fc for sampling frequency Fs.
% The filter is designed according to the Order-N specification
% of the ANSI S1.1-1986 standard. Default value for N is 3.
% Warning: for meaningful design results, center frequency used
% should preferably be in range Fs/200 < Fc < Fs/5.
% Usage of the filter: Y = FILTER(B,A,X).
% Note: Requires the Signal Processing Toolbox.
% Author: Christophe Couvreur, Faculte Polytechnique de Mons (Belgium)
% couvreur@thor.fpms.ac.be
% Last modification: Aug. 25, 1997, 2:00pm.

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
% the ANSI standard. BUTTER is part of the Signal Processing Toolbox.
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]);

clc;clear;
Fs = 44100;
time=0:1/Fs:1;

% x = 1*10^5.*rand(1, Fs); % 生成 N 个均匀分布的随机数
% EXAMPLE.M Example of frequency analysis
x=10^4.*sin(2*pi*1000*time);
n = length(x);

% Perform 1/3-octave band analysis.
% The plot is automatically generated by OCT3BANK.
figure(1)
oct3bank(x);
% Perform FFT-based analysis & generate separate plot.
X = fft(x)/n;
Pfft = 20*log10(abs(X(1:n/2)));
Ffft = (0:n/2-1)*Fs/n;
figure(2)
plot(Ffft,Pfft);
xlabel('Frequency [Hz]'); ylabel('Power Spectral Density [dB]');
title('DFT spectrum')

生成的是一个1000hz单频音，运行后得到如下结果，在相邻频带实现了接近20dB的衰减：

2.个人编写版本(不包含Zwicker模型中的下降沿)

clc;clear;
p_ref=2*10^(-5);
Fs = 44100; % Sampling Frequency
time=0:1/Fs:1;
% x = 1*ones(1, Fs); % 生成 N 个均匀分布的随机数
% x=0.02.*sin(2*pi*1000.*time);
x =20.* pinknoise(1*Fs);%一段粉噪声
% x =1.*rand(1,1*Fs);%一段白噪声

N = 3; % Order of analysis filters.
F = [ 100 125 160, 200 250 315, 400 500 630, 800 1000 1250, ...
    1600 2000 2500, 3150 4000 5000, 6300 8000 10000, 12500 ]; % Preferred labeling freq.
ff = (1000).*((2^(1/3)).^[-10:11]); % Exact center freq.
P = zeros(1,22);
m = length(x);

%%%%使用fft analysis
X = fft(x,m);
tol = 1e-6;
X(abs(X)<=tol)=0;
P=abs(X)/m*2;
SPL=20*log10(P./p_ref);
%%%Power_sum 
SPL_ave=10*log10(sum(x.^2)/m/p_ref.^2*2)
Ffft = (0:m/2-1)*Fs/m;
figure(1)
semilogx(Ffft,SPL(1:m/2));
xlim([80,20000])
xlabel('Frequency [Hz]'); ylabel('SPL [dB]');
title('FFT spectrum')

%%%%使用1/3 octave analysis
[B,A] = oct3dsgn(ff(11),Fs,N);
y = filter(B,A,x);
Y = fft(y);
tol = 1e-6;
Y(abs(Y)<=tol)=0;
P_2=abs(Y)/m*2;
SPL_2=20*log10(P_2./p_ref);
%%%Power sum
SPL_ave2=10*log10(sum(y.^2)/m/p_ref^2*2)
figure(2)
semilogx(Ffft,SPL_2(1:m/2))
xlim([80,20000])
xlabel('Frequency [Hz]'); ylabel('SPL [dB]');
title('FFT spectrum filtered')

%%bar of 1/3 ocatave power
Power=zeros(1,22);
for i = 1:1:22
    [B,A] = oct3dsgn(ff(i),Fs,N);
    y = filter(B,A,x);
    Power(i) = sum(y.^2)/m;
end
idx = find(Power>0);
Power(idx) = 10*log10(Power(idx)/p_ref^2*2);
Power(~idx) = NaN*ones(sum(~idx),1);
figure(3)
bar(Power)
set(gca,'XTick',[2:3:22]); % Label frequency axis on octaves.
set(gca,'XTickLabel',F(2:3:length(F))); % Syntax for MATLAB 5.1
xlabel('Frequency band [Hz]'); ylabel('Power [dB]');
title('One-third-octave spectrum')

%%%1/3 octave band to Bark band
%%threshold factor
s=0.25;
%%frequency of 1/3 octave,no.17 is 1000Hz
fr=[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];
%%Bark no.9 = 1000Hz
Bark=[0.9,1.8,2.8,3.5,4.4,5.4,6.6,7.9,9.2,...
    10.6,12.3,13.8,15.2,16.7,18.1,19.3,20.6,21.8,22.7,23.6];
%%Critical band level at the threshold in quiet 
% L_TQ=[30,30,30,30,30,30,18,18,18,12,12,8,7,6,5,4,3,3,3,3,3,3,3,3,3,3,3,3];
%%Bark band level at the threshold in quiet 
L_TQ=[30,18,12,8,7,6,5,4,3,3,3,3,3,3,3,3,3,3,3,3];
%%SPL level of critical band
L_CB=zeros(1,length(Bark));
L_CB(1)=30;
L_CB(2)=10.*log10(10.^(Power(1)/10)+10.^(Power(2)/10)+10.^(Power(3)/10));
L_CB(3)=10.*log10(10.^(Power(4)/10)+10.^(Power(5)/10));
L_CB(4:end)=Power(6:end);
N_c=zeros(1,length(Bark));
alpha=zeros(1,length(Bark));
for i=1:length(Bark)
    alpha(i)=1-s+s*10^(0.1*(L_CB(i)-L_TQ(i)));
    N_c(i)=0.0635*10^(0.025*L_TQ(i))*(alpha(i).^(0.25)-1);
end

tol = 1e-2;
N_c(abs(N_c)<=tol)=0;
%%%画出bar图
Bark2=zeros(1,2*length(Bark)+1);
N_c2=zeros(1,2*length(N_c)+1);

Bark2(1)=0;
N_c2(1)=N_c(1);

Result=Bark(1)*N_c(1);%total Specific loudness
for k=1:(length(Bark)-1)
    Bark2(2*k)=Bark(k);
    Bark2(2*k+1)=Bark(k);
    N_c2(2*k)=N_c(k);
    N_c2(2*k+1)=N_c(k+1);
    Result=Result+(Bark(k+1)-Bark(k)).*N_c(k+1);
end
Result=Result+(24-Bark(end))*N_c(20);
Bark2(40)=23.6;
N_c2(40)=N_c(20);
Bark2(41)=24;
N_c2(41)=N_c(20);
figure(4)
plot(Bark2,N_c2)
xlim([0,24])
ylim([0,5])
set(gca,'XTick',[0:2.5:24]); % Label frequency axis on octaves.
grid on
hold on
line([0,25],[0,0])
xlabel('Bark'); ylabel('Specific Loudness [Sone/Bark]');

%Loundness level (sone) to loudness (phon) 
if Result>=1
    L=40+33.22*log10(Result);
else
    L=40*(Result+0.0005)^0.35;
end

disp(sprintf('Loudness: %f (sone), \nLoudness: %f (phon)', Result, L));

代码中使用了一段粉红噪声进行测试，倍频谱能量以及Bark谱，响度如下：

3.MATLAB内置函数(重要！！！可惜前几天一直都没找到/(ㄒoㄒ)/)

做了很多天和Zwicker方法有关的调查才发现，matlab居然有内置函数：acousticLoudness

详情可见：Perceived loudness of acoustic signal - MATLAB acousticLoudness

这里只展示最简单的用法，给了一个粉红噪声：

fs = 48e3;
dur = 5;
pnoise = 2*pinknoise(dur*fs);
figure
acousticLoudness(pnoise,fs)

运行结果：

SoundCheck内置函数

声学测试软件SC中post-processing中内置有Zwicker响度函数，也可用来计算响度

Too_L

关注

18
点赞
踩
14

收藏

觉得还不错? 一键收藏
0
评论
Zwicker响度-Methods for calculating loudness

最近做了一些和Zwicker响度方面有关的调研，此文主要记录个人的一些日常总结，和分享相关的资料，如果只需要最终结论可以直接跳到Matlab实现中的第三节，有相关代码！ISO532系列标准的目的是指定基于声音物理特性的计算程序,以估计在特定聆听条件下耳科正常听力的人所感知的声音的响度和响度级别。ISO532-1和ISO532-2规定了两种不同的响度计算方法,对于给定的声音可能会产生不同的结果。由于目前无法说明对一种或另一种方法的一般偏好,因此由用户选择最适合给定情况的方法。
复制链接

扫一扫