Matlab实现计算声强与响度

Matlab实现计算声强与响度

实验环境：

matalb2019版本

一.声压计算

1.声压

声压是定量描述声波的最基本的物理量，它是由于声扰动产生的逾量压强，是空间位置和时间的函数。

2.有效声压

通常讲的声压指的是有效声压，即在一定时间间隔内将瞬时声压对时间求方均根值所 得。设语音长度为T, 离散点数为N, 则有效声压的计算公式为：
$$P_e=\sqrt{\frac{1}{T}\sum _{n=1}^N{x}^2\Delta t}=\sqrt{\frac{1}{N\Delta t}\sum _{n=1}^N{x}^2\Delta t}=\sqrt{\frac{1}{N}\sum _{n=1}^N{x}^2}$$
其中，x 表示语音信号的采样点。只要保证所取的点数N足够大，即可保证计算准确性。

3.声压级

声音的有效声压与基准声压之比，取以10 为底的对数，再乘以20, 即为声压级，通常以符号L, 表示，单位为dB。
$$L_p=20\lg \frac{P_e}{P_{\left(ref\right)}}\left(dB\right)$$
式中，pe为待测声压的有效值； Pref为参考声压，在空气中参考声压一般取2 x 10-5 Pa。

4.编写声压级计算函数

使用matlab编写声压级计算函数

%声压级计算函数
function spl = SPLCal( x,fs,flen )
% SPLCal函数用于计算一帧语音信号的声压级值
% 输入参数x为输入的语音信号
% 这里用x的值直接替代瞬时声压的值，理论上瞬时声压值跟x的值是成正比的
% 输入参数fs为采样率，单位为Hz
% flen为一帧信号的时间长度，单位为ms
% spl为函数输出的声压级值
% 输入的语音帧的长度
Length = length(x);
% 每一帧信号的离散点数
M = flen*fs/1000;
if Length~=M
   error('输入信号长度与所定义帧长不等！'); 
end
%%-------------------计算有效声压-----------------
% 根据定义计算有效声压，pa = sqrt((x(1)^2+x(2)^2+...+x(M)^2)/M)
% 单位为Pa
pp = 0;
for i = 1:M
   pp = pp + x(i)^2;
end
pa = sqrt(pp/M);
%%----------------计算声压级---------------------
% 声压级值spl=20*log10(pa/p0),单位为dB
% 基准声压p0，单位为Pa
p0 = 2*10^-5;
spl = 20*log10(pa/p0);
end

使用上面的函数计算有效声压

%实验要求一：计算有效声压
clear all;
%%-------------------读取文件--------------------
% 从wav文件读入语音数据，该语音采样率16k，故信号最高频率8k。
[x,fs]=audioread('E:\语音处理\配套代码\第02章 语音信号处理基础\2.3 声强与响度实验\C2_3_y.wav');
% 取x的一个通道
x=x(:,1);
% 将x从列向量转为行向量
x=x';
% max(x)
% 得到向量x的长度
Length=length(x);
%%-------------------语音分帧--------------------
% 每帧大小为M，当语音长度不是帧长的整数倍时：
% (1)若剩余长度大于等于帧长的二分之一，则补零至帧长
% (2)若剩余长度小于帧长的二分之一，则舍弃
% 用于计算声压级值的语音帧长分别为20ms、50ms、100ms、200ms以及500ms
framlen = 100;
% 每帧信号的离散点数
M=fs*framlen/1000;
% m为Length/M后得到的余数
m = mod(Length,M);
if m >= M/2 % 补零
    % 补零后的语音
    x = [x,zeros(1,M-m)];
    % 补零后的语音帧长
    Length = length(x);
else   % 即m < M/2，则将剩余的语音帧舍弃
    % l为Length/M后得到的商
    l = floor(Length/M);
    % 舍弃后的语音
    x = x(1,1:M*l);
    % 舍弃后的语音帧长
    Length = length(x);
end
% 最终的语音分帧总帧数
N = Length/M;
%%------------------------计算声压级-------------
s = zeros(1,N);
% N帧信号的声压级值存储在spl向量里
spl = zeros(1,N);
for k = 1:N
    % 读取第k帧信号
    s =x((k-1)*M + 1:k*M);
    % 计算第k帧信号的声压级值
    spl(1,k) = SPLCal(s,fs,framlen);
end
%%------------画图------
t = 1:Length;
SPL = zeros(1,Length);
for r = 1:N
    SPL(1,(r-1)*M+1:r*M) = spl(r);
end
figure(1)
subplot(211)
plot(t/fs,x);
grid on
xlabel('时间(s)');
title('输入语音波形');
subplot(212)
stairs(t/fs,SPL,'r');
grid on
xlabel('时间(s)');
ylabel('声压级(dB)');
title('语音信号的声压级(dB)/100ms');

运行结果

声压计算间隔为100ms

二.强度计算

1.强度

响度描述的是声音的响亮程度，表示人耳对声音的主观感受，其计量单位是宋。定义为声压级为40dB 的1 kHz纯音的响度为I Son (宋）。人耳对声音的感觉，不仅和声压有关，还和频率有关。声压级相同，频率不同的声音，听起来响亮程度也不同。

2.等响度曲线

人耳对不同频率的纯音有不同的敏感度。1kHz的纯音在0dBHL时就可被察觉；而20Hz的纯音则在70dBHL时才可被察觉。

根据1S0226 —2003标准H等响度曲线的定义如下。假设频率为J的纯音的响度级为Lv, 则其声压级LP 为：
$$L_p=\left(\frac{10}{a_f}.\lg A_f-L_u+94\right)\left(dB\right)$$

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

$$

Tf为听力阈值；af为响度感知指数；Lv为以1000Hz为标准所计算的线性传输函数的幅值。

3.响度计算函数的编写

根据计算公式，使用matlab进行响度计算函数编写

%响度计算函数
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=20:10:100;
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

运行结果：



由图可知， 当声压级在80dB以上时，各个频率的声压级与响度级的数值就比较接近了，这表明当声压级较高时，人耳对各个 频率的声音的感觉基本是一样的。