MATLAB之声压级计算

最新推荐文章于 2023-03-06 18:51:44 发布

Heart_Sea

最新推荐文章于 2023-03-06 18:51:44 发布

阅读量1.1w

点赞数 9

分类专栏：信号处理

本文链接：https://blog.csdn.net/heart_sea/article/details/89529517

版权

信号处理专栏收录该内容

10 篇文章 6 订阅

订阅专栏

文章目录

声压级计算

声压级计算

一个文件夹下多个文件的SPL

执行文件exert

% 功能：遍历data文件夹下的所有文件，计算带A计权的声压级曲线，并在命令行显示SPL（dB）

% 编辑者：Heart_sea
% 日期：2019,4,24

clear;
clc;
close all;

% ======================= preferences set =================================
fs = 12000;          % sampling frequency (unit: Hz)
NFFT = fs;           % the number of samples used in FFT
pl = 20;             % 下限频率
ph = 6000;           % 上限频率
N_win = NFFT;        % the length of the window，% 规定窗长度是NFFT,暂时这样理解，后面学窗再详细探讨
Window =  hamming(N_win);      % 或者 WINDOW =  NFFT;默认是汉明窗，只给出长度即可
N_overlap = floor(N_win/2);    % 重叠是窗长度的一半
A_weighting_flag = 2;          % 是否要A计权,0表示不计权，1表示国标的A计权，2表示调用库的A计权
save_image_flag =0;            % 是否保存图片,0不存
timecom = 1;                   % 1 需显示时域曲线，2 无需显示时域曲线
coordinate = 1;                % 坐标选择，0表示直角坐标，1表示x轴对数坐标
% ===========计算麦克风的灵敏度
OP_AMP = 5.1;                  % 麦克风板的放大倍数，采集的是经过放大倍数的数据
MIC_SEN = -38;                 % 灵敏度级dB
mic_sen = (10^(MIC_SEN./20)).*OP_AMP;%麦克风的灵敏度
% ======================== import data ====================================
maindir = './Data';
File_L1 =dir(maindir);    % 获取路径下的所有文件信息
PathNumber=numel(File_L1);% 获取元素的数量
m = 1;                    % 存放数据列数初始化
switch timecom
    case 1
        for i = 3:PathNumber
            Path=fullfile(maindir,File_L1(i).name);
            Path2 = fullfile(File_L1(i).name);
            Path2(end-3:end) = [];

            note(:,m) = {[File_L1(i,1).name]};
            m = m+1;

            x = load(Path);   % importdata(Path)
            x = x./mic_sen;
            t = 1/fs:1/fs:length(x)/fs;
            plot(t,x)
            hold on;
            legend(note);
        end
        title(['时域曲线 '])
        xlabel('t (s)')
        ylabel('pressure (Pa)')
      case 2
end

if timecom == 1
    figure;
end
% ======================== 频域 ====================================
m = 1; %列放数据初始化
for i = 3:PathNumber
    Path=fullfile(maindir,File_L1(i).name);% 文件夹下文件路径
    Path2 = fullfile(File_L1(i).name);  % 文件夹下的文件的名字
    Path2(end-3:end) = [];              %去掉文件的txt
    note(:,m) = {[File_L1(i,1).name]}; % 用作legend
    m = m+1;
    
    x = load(Path);
    x = x./mic_sen;
 
    [Sxx,f] = cpsd(x,x, Window, N_overlap, NFFT, fs);
    f = f(pl:ph);
    Sxx = Sxx(pl:ph);
% ======================== A-weighting-国标 ===============================
    if A_weighting_flag == 1 
        WeightA = A_curve(f); % get the vaule of weight-A
        SPL = 10*log10(Sxx/4e-10)+ WeightA;%各频率上的A计权声压级
        SPLA_S = sum(10.^(SPL./10));
        SPLA_total = 10*log10(SPLA_S);
        
% ======================== A-weighting-调用MATLAB库 ===============================
    elseif A_weighting_flag == 2 
        N_f = sum(f<=fs/2);
        Weighting = 'A';%Weighting支持‘A’计权或者‘C’计权
    % The valid weighting types are: A weighting, C weighting, C-message, ITU-T 0.41, and ITU-R 468–4 weighting. 
    % Filter class is only applicable for A weighting and C weighting filters. 
    % https://ww2.mathworks.cn/help/dsp/ref/fdesign.audioweighting.html
        Filter_Obj = fdesign.audioweighting('WT,Class',Weighting,1,fs);
        H_Filter = design(Filter_Obj);
        [Hf,~] = freqz(H_Filter,N_f,fs);
        for i = 1:N_f
            Sxx(i) = Sxx(i)*((abs(Hf(i)))^2);
        end 
        Sxx_total = sum(Sxx);
        SPLA_total = 10*log10(Sxx_total/4e-10);
        SPL = 10*log10(Sxx/4e-10);
        
% ======================== 不计权
    elseif A_weighting_flag == 0
        Sxx_total = sum(Sxx);
        SPLA_total = 10*log10(Sxx_total/4e-10);
        SPL = 10*log10(Sxx/4e-10);
    end   
    
       disp([num2str(Path2) ':   ' num2str(SPLA_total) 'dB'])
% ===================== 坐标选择 
     if coordinate == 0   % 0代表直角坐标
         plot(f,SPL)
         hold on;
    elseif coordinate == 1
         semilogx(f,SPL); 
         set(gca, 'XTick',[20,50,100,200,500,1e3,2e3,6e3] );
%        set(gca, 'XTick',[20,50,100,200,500,1e3,2e3,5e3,10e3,20e3] ); % X轴的记号点
         set(gca,'xticklabel',get(gca,'xtick'),'yticklabel',get(gca,'ytick'));
         hold on;
     end
    
    if  A_weighting_flag == 0 
        title(['SPL '])
    else
        title(['SPLA '])
    end  
    xlabel('f (Hz)')
    ylabel('(dB)')
end
legend(note)
    if save_image_flag == 1
        set(gca,'FontSize',12)
        saveas(gcf,'SPL.fig')  %可以保存成jpg，fig等不同格式的图片
    end

国标计权A_curve

% 功能：根据国标模拟滤波估计计算A计权曲线，单位为dB
% 输入  f：频率
% 输出  y：A计权曲线对应的值
function [y] = A_curve(f) 
    n1 = 20.6;
    n2 = 107.7;
    n3 = 737.9;
    n4 = 12194;

    p  = n4^2 .* f.^4;
    p1 = n1^2 + f.^2;
    p2 = sqrt(n2^2 + f.^2);
    p3 = sqrt(n3^2 + f.^2);
    p4 = n4^2 + f.^2;

    y  = 20.*log10(p./(p1 .* p2 .* p3 .* p4))+2;
 
end

一个文件夹下多组文件的SPL差

执行文件exert

% 功能：遍历data文件夹下的所有文件（off-on），计算带A计权的每组的声压级曲线，并在命令行显示SPL（dB）
% 编辑者：lily
% 日期：2019,4,25
clear;
clc;
close all;
% ======================= preferences set =================================
fs = 12000;          % sampling frequency (unit: Hz)
NFFT = fs;           % the number of samples used in FFT
fl = 20;             % 下限频率
fh = 6000;           % 上限频率
point_number = 540000;   % the number of sample points 
jump_line = 0;           % 跳过第jump_line行开始取，0表示从第一行取
micbegin = 1;            % 文件序号的第几个开始
micEnd = 4;              % 文件序号的第几个结束
N_win = NFFT;            % the length of the window，% 规定窗长度是NFFT,暂时这样理解，后面学窗再详细探讨
Window =  hamming(N_win);      % 或者 WINDOW =  NFFT;默认是汉明窗，只给出长度即可
N_overlap = floor(N_win/2);    % 重叠是窗长度的一半
A_weighting_flag = 1;          % 是否要A计权,0表示不计权，1表示国标的A计权，2表示调用库的A计权
save_image_flag =0;            % 是否保存图片,0不存
coordinate = 0;                % 坐标选择，0表示直角坐标，1表示x轴对数坐标
smooth = 1;                    % the switch of SPL smooth, ON:1 OFF:0
% ===========计算麦克风的灵敏度
% 灵敏度与降噪前后的dB值有关，但与降噪前后的dB差值无关
OP_AMP = 5.1;                  % 麦克风板的放大倍数，采集的是经过放大倍数的数据
MIC_SEN = -38;                 % 灵敏度级dB
mic_sen = (10^(MIC_SEN./20)).*OP_AMP;%麦克风的灵敏度

% ======================== import data ====================================
figure(1)
for micNumber = micbegin:micEnd
    filename = ['err mic' num2str(micNumber) ' - off.txt'];
    off = ReadData_Raw(filename,point_number,jump_line);
    off = off /mic_sen;%电压/灵敏度=声压
    [S_off,f] = cpsd(off,off,Window,N_overlap,NFFT,fs);

    filename = ['err mic' num2str(micNumber) ' - on.txt'];
    on = ReadData_Raw(filename,point_number,jump_line);
    on = on/mic_sen;
    [S_on,~] = cpsd(on,on,Window,N_overlap,NFFT,fs);
    
    f = f(fl+1:fh+1);
    S_off = S_off(fl+1:fh+1); 
    S_on = S_on(fl+1:fh+1); 
    
% ======================== A-weighting-国标 ===============================
%     此处关键点：dB相加
    if A_weighting_flag == 1 
        WeightA = A_curve(f); % get the vaule of weight-A
        SPL_off = 10*log10(S_off/4e-10)+ WeightA;%各频率上的A计权声压级
        SPL_on = 10*log10(S_on/4e-10)+ WeightA;%各频率上的A计权声压级
        
        SPL_off_sum = sum(10.^(SPL_off./10));
        SPL_on_sum = sum(10.^(SPL_on./10));
        
        SPL_total_off = 10*log10(SPL_off_sum);
        SPL_total_on = 10*log10(SPL_on_sum);
        
  % ======================== A-weighting-调用MATLAB库 ===============================
    elseif A_weighting_flag == 2 
        N_f = sum(f<=fs/2);%返回满足条件的数据个数
            Weighting = 'A';%Weighting支持‘A’计权或者‘C’计权
            Filter_Obj = fdesign.audioweighting('WT,Class',Weighting,1,fs);
            H_Filter = design(Filter_Obj);
            [Hf,~] = freqz(H_Filter,N_f,fs);

            for i = 1:N_f
                S_off(i) = S_off(i)*((abs(Hf(i)))^2);
                S_on(i) = S_on(i)*((abs(Hf(i)))^2);
            end

        Syy_off_sum = sum(S_off);
        Syy_on_sum = sum(S_on);

        SPL_total_off = 10*log10(Syy_off_sum/4e-10);
        SPL_total_on = 10*log10(Syy_on_sum/4e-10);

        SPL_off = 10*log10(S_off/4e-10);
        SPL_on = 10*log10(S_on/4e-10);
  % ======================== 不计权
    elseif A_weighting_flag == 0
         Syy_off_sum = sum(S_off);
         Syy_on_sum = sum(S_on);
        
         SPL_total_off = 10*log10(Syy_off_sum/4e-10);
         SPL_total_on = 10*log10(Syy_on_sum/4e-10);
         
         SPL_off = 10*log10(S_off/4e-10);
         SPL_on = 10*log10(S_on/4e-10);
    end
 % ============================= 坐标选择 
     if coordinate == 0   % 0代表直角坐标
 % ======是否需要平滑滤波
        if smooth == 0 %不需滤波
            subplot (micEnd,1,micNumber)
            plot(f,SPL_off);
            hold on
            plot(f,SPL_on)
            grid on
        elseif  smooth == 1  %平滑滤波  
            subplot (micEnd,1,micNumber)
            % csaps函数可以给定各种边界条件的三次样条插值函数
            % csaps（x,y,z）实现光滑拟合，z为权因子，0<z<1,z值越大，越接近真实值，z=0实现线性拟合，z=1实现自然线条
            z = 1E-4;
            S_off = csaps(f,SPL_off,z);
            fnplt(S_off,'b');  %降噪前off
            hold on
            S_on = csaps(f,SPL_on,z);
            fnplt(S_on,'r');  % 降噪后oN
%             set(gca,'YTick',-40:10:50);
%             set(gca,'XTick',fl:500:fh);
            grid on
        end

    elseif coordinate == 1   % 1代表x对数坐标
% ==========是否需要平滑滤波
         if smooth == 0  
             subplot (micEnd,1,micNumber)
             semilogx(f,SPL_off); 
             hold on
             semilogx(f,SPL_on); 
             set(gca, 'XTick',[20,50,100,200,500,1e3,2e3,6e3] );
    %        set(gca, 'XTick',[20,50,100,200,500,1e3,2e3,5e3,10e3,20e3] ); % X轴的记号点
             set(gca,'xticklabel',get(gca,'xtick'),'yticklabel',get(gca,'ytick'));
             grid on
         elseif smooth == 1  %平滑滤波
                subplot (micEnd,1,micNumber)
                % csaps函数可以给定各种边界条件的三次样条插值函数
                % csaps（x,y,z）实现光滑拟合，z为权因子，0<z<1,z值越大，越接近真实值，z=0实现线性拟合，z=1实现自然线条
                z = 1E-4;
                S_off = csaps(f,SPL_off,z);
                semilogx(f,([S_off.coefs(:,end);S_off.coefs(end,end)]),'b');  % 绘制样条函数图形
                hold on
                S_on = csaps(f,SPL_on,z);
                semilogx(f,([S_on.coefs(:,end);S_on.coefs(end,end)]),'r');
                set(gca, 'XTick',[20,50,100,200,500,1e3,2e3,6e3] );
                set(gca,'xticklabel',get(gca,'xtick'),'yticklabel',get(gca,'ytick'));
                grid on
         end
     end
        k = SPL_total_off - SPL_total_on;%s是off，c是on
        disp(['err' num2str(micNumber) ' = ' num2str(k)  'dB'] );
        disp(['第' num2str(micNumber) '组SPL：降噪前=' num2str(SPL_total_off) 'dB' '   降噪后=' num2str(SPL_total_on) 'dB']);
     
% ============= 标题     
        if  A_weighting_flag == 0 
            title(['SPL ' num2str(micNumber)])
        else
            title(['SPLA ' num2str(micNumber)])
        end  
        xlabel('f (Hz)')
        ylabel('(dB)')
end
% =============保存图片
if save_image_flag == 1
    set(gca,'FontSize',12)
    saveas(gcf,'ANCEffectFig.fig')  %可以保存成jpg，fig等不同格式的图片
end

一个总文件夹下多个子文件夹，每个子文件夹有多组文件，计算spl差。

执行文件exert

## 一个总文件夹下多个子文件夹，每个子文件夹有多组文件
### 执行文件exert
```java
% 功能：遍历data文件夹下的所有子文件夹下的文件（off-on），计算带A计权的每组的声压级曲线，并在命令行显示SPL（dB）
% 编辑者：lily
% 日期：2019,4,25
clear;
clc;
close all;
% ======================= preferences set =================================
maindir = './Data';
fs = 12000;          % sampling frequency (unit: Hz)
NFFT = fs;           % the number of samples used in FFT
fl = 20;             % 下限频率
fh = 6000;           % 上限频率
max = 4;             % 子文件夹最多的文件个数
point_number = 80000;   % the number of sample points 
jump_line = 0;           % 跳过第jump_line行开始取，0表示从第一行取
N_win = NFFT;            % the length of the window，% 规定窗长度是NFFT,暂时这样理解，后面学窗再详细探讨
Window =  hamming(N_win);      % 或者 WINDOW =  NFFT;默认是汉明窗，只给出长度即可
N_overlap = floor(N_win/2);    % 重叠是窗长度的一半
A_weighting_flag = 1;          % 是否要A计权,0表示不计权，1表示国标的A计权，2表示调用库的A计权
save_image_flag =0;            % 是否保存图片,0不存
coordinate = 1;                % 坐标选择，0表示直角坐标，1表示x轴对数坐标
smooth = 1;                    % the switch of SPL smooth, ON:1 OFF:0
% ========计算麦克风的灵敏度
% 灵敏度与降噪前后的dB值有关，但与降噪前后的dB差值无关
OP_AMP = 5.1;                  % 麦克风板的放大倍数，采集的是经过放大倍数的数据
MIC_SEN = -38;                 % 灵敏度级dB
mic_sen = (10^(MIC_SEN./20)).*OP_AMP;%麦克风的灵敏度

% ======================== 循环 import data ===============================
[Note_legend,~] = importdata(maindir,fs,NFFT,fl,fh,max,point_number,jump_line,...
    N_win,Window,N_overlap,A_weighting_flag,coordinate,smooth,mic_sen);
% =============保存图片
if save_image_flag == 1
    set(gca,'FontSize',12)
    saveas(gcf,'ANCEffectFig.fig')  %可以保存成jpg，fig等不同格式的图片
end

调用函数importdata


function [Note_legend,Note_disp] = importdata(maindir,fs,NFFT,fl,fh,max,...
    point_number,jump_line,~,Window,N_overlap,A_weighting_flag,coordinate,smooth,mic_sen)

File_L1 =dir(maindir);    % 获取data里的内容
PathNumber=numel(File_L1);% 获取结构体的行数

for i = 3:PathNumber      %读取数据从第三行到最后一行
    figini = 2;%fig初始化,3-2=1,每一组都要初始化
    Path = fullfile(maindir,File_L1(i).name);%前两阶文件夹的名称，字符串格式
% ==============设置lengend
    Path2 = fullfile(File_L1(i).name);%获取第二层文件夹的名字，字符串格式
    dat = dir(Path);%获取第二层文件夹的所有txt文件，struct结构体格式
    row = 1; 
    for  j = 3 : length( dat )   %循环dat名称
        str = strcat(Path2,fullfile(dat( j ).name)); %  str获取文件夹及其文件名字，并连接
        str(end-3:end) = [];                         % 将字符串的后四个字符串置成空的，去掉.txt
        Note_legend(row,i-2) = {str}; 
        row = row+1;
    end
   
% ======================== 循环dat数据（子文件夹里的文件） ================
    colum =1;  % 存放数据到列初始化
    v = 2;     % subplot第几个fig初始化
    for m = 3:2:length( dat )
% ==============设置disp
        str1 = strcat(Path2,fullfile(dat( m ).name)); 
        str1(end-8:end) = [];
        Note_disp (:,colum) = {str1 };
        colum = colum+1; 
% ==============开始导入数据
%     off = load(fullfile(Path,dat(m).name));
        filename = fullfile(Path,dat(m).name);
        off = ReadData_Raw(filename,point_number,jump_line);
        off = off /mic_sen;%电压/灵敏度=声压
        [S_off,f] = cpsd(off,off,Window,N_overlap,NFFT,fs);
  
        filename = fullfile(Path,dat(m+1).name);
        on = ReadData_Raw(filename,point_number,jump_line);
        on = on /mic_sen;
        [S_on,~] = cpsd(on,on,Window,N_overlap,NFFT,fs);

        f = f(fl+1:fh+1);
        S_off = S_off(fl+1:fh+1); 
        S_on = S_on(fl+1:fh+1);
% ======================== A-weighting-国标 ===============================  
%     此处关键点：dB相加
    if A_weighting_flag == 1 
        WeightA = A_curve(f); % get the vaule of weight-A
        SPL_off = 10*log10(S_off/4e-10)+ WeightA;%各频率上的A计权声压级
        SPL_on = 10*log10(S_on/4e-10)+ WeightA;%各频率上的A计权声压级
        
        SPL_off_sum = sum(10.^(SPL_off./10));
        SPL_on_sum = sum(10.^(SPL_on./10));
        
        SPL_total_off = 10*log10(SPL_off_sum);
        SPL_total_on = 10*log10(SPL_on_sum);
 % ======================== A-weighting-调用MATLAB库 ===============================
    elseif A_weighting_flag == 2 
        N_f = sum(f<=fs/2);%返回满足条件的数据个数
            Weighting = 'A';%Weighting支持‘A’计权或者‘C’计权
            Filter_Obj = fdesign.audioweighting('WT,Class',Weighting,1,fs);
            H_Filter = design(Filter_Obj);
            [Hf,~] = freqz(H_Filter,N_f,fs);

            for i = 1:N_f
                S_off(i) = S_off(i)*((abs(Hf(i)))^2);
                S_on(i) = S_on(i)*((abs(Hf(i)))^2);
            end

        Syy_off_sum = sum(S_off);
        Syy_on_sum = sum(S_on);

        SPL_total_off = 10*log10(Syy_off_sum/4e-10);
        SPL_total_on = 10*log10(Syy_on_sum/4e-10);

        SPL_off = 10*log10(S_off/4e-10);
        SPL_on = 10*log10(S_on/4e-10);
        
  % ======================== 不计权
    elseif A_weighting_flag == 0
         Syy_off_sum = sum(S_off);
         Syy_on_sum = sum(S_on);
        
         SPL_total_off = 10*log10(Syy_off_sum/4e-10);
         SPL_total_on = 10*log10(Syy_on_sum/4e-10);
         
         SPL_off = 10*log10(S_off/4e-10);
         SPL_on = 10*log10(S_on/4e-10);
    end
    
     % ============================= 坐标选择 
     fig = m-v;% fig的第几个。m=3时，fig=1，m=5时，fig=2,依次类推
     if coordinate == 0   % 0代表直角坐标
        if smooth == 0 %不需滤波
            subplot (max,1,fig)
            plot(f,SPL_off);
            hold on
            plot(f,SPL_on)
            grid on
        elseif  smooth == 1  %平滑滤波  
            subplot (max,1,fig)
            % csaps函数可以给定各种边界条件的三次样条插值函数
            % csaps（x,y,z）实现光滑拟合，z为权因子，0<z<1,z值越大，越接近真实值，z=0实现线性拟合，z=1实现自然线条
            z = 1E-4;
            S_off = csaps(f,SPL_off,z);
            fnplt(S_off,'');  %降噪前off
            hold on
            S_on = csaps(f,SPL_on,z);
            fnplt(S_on,'');  % 降噪后oN
%             set(gca,'YTick',-40:10:50);
%             set(gca,'XTick',fl:500:fh);
            grid on
        end

    elseif coordinate == 1   % 1代表x对数坐标
         if smooth == 0  
             subplot (max,1,fig)
             semilogx(f,SPL_off); 
             hold on
             semilogx(f,SPL_on); 
             set(gca, 'XTick',[20,50,100,200,500,1e3,2e3,6e3] );
    %        set(gca, 'XTick',[20,50,100,200,500,1e3,2e3,5e3,10e3,20e3] ); % X轴的记号点
             set(gca,'xticklabel',get(gca,'xtick'),'yticklabel',get(gca,'ytick'));
             axis([fl,fh,-inf,inf]);
             grid on
         elseif smooth == 1  %平滑滤波
                subplot (max,1,fig)
                % csaps函数可以给定各种边界条件的三次样条插值函数
                % csaps（x,y,z）实现光滑拟合，z为权因子，0<z<1,z值越大，越接近真实值，z=0实现线性拟合，z=1实现自然线条
                z = 1E-4;
                S_off = csaps(f,SPL_off,z);
                semilogx(f,([S_off.coefs(:,end);S_off.coefs(end,end)]),'');  % 绘制样条函数图形
                hold on
                S_on = csaps(f,SPL_on,z);
                semilogx(f,([S_on.coefs(:,end);S_on.coefs(end,end)]),'');
                set(gca, 'XTick',[20,50,100,200,500,1e3,2e3,6e3] );
                set(gca,'xticklabel',get(gca,'xtick'),'yticklabel',get(gca,'ytick'));
                axis([fl,fh,-inf,inf]);
                grid on
         end
     end
        legend(char(Note_legend(m-2:m-2+1,:)))
        v = v+1;
        k = SPL_total_off - SPL_total_on;%s是off，c是on
        disp([num2str(Note_disp{fig}) ':' num2str(k) 'dB'])    
        
   % ============= 标题     
        if  A_weighting_flag == 0 
            title(['SPL ' num2str(fig)])
        else
            title(['SPLA ' num2str(fig)])
        end  
        xlabel('f (Hz)')
        ylabel('(dB)')
    end
end
end

Heart_Sea

关注

9
点赞
踩
90

收藏

觉得还不错? 一键收藏
21
评论
MATLAB之声压级计算

文章目录声压级计算一个文件夹下多个文件的SPL执行文件exert国标计权A_curve一个文件夹下多组文件的SPL差执行文件exert一个总文件夹下多个子文件夹，每个子文件夹有多组文件，计算spl差。执行文件exert调用函数importdata声压级计算一个文件夹下多个文件的SPL执行文件exert% 功能：遍历data文件夹下的所有文件，计算带A计权的声压级曲线，并在命令行显示SPL（...
复制链接

扫一扫