MATLAB之声压级计算

￡神魔★判官ぃ 2022-02-14 07:23 378阅读 0赞

### 文章目录 ###

*  声压级计算
 *   *  一个文件夹下多个文件的SPL
     *   *  执行文件exert
         *  国标计权A\_curve
     *  一个文件夹下多组文件的SPL差
     *   *  执行文件exert
     *  一个总文件夹下多个子文件夹，每个子文件夹有多组文件，计算spl差。
     *   *  执行文件exert
         *  调用函数importdata

# 声压级计算 #

## 一个文件夹下多个文件的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