MATLAB之声压级计算

声压级计算

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