文章目录
声压级计算一个文件夹下多个文件的SPL执行文件exert国标计权A_curve一个文件夹下多组文件的SPL差执行文件exert一个总文件夹下多个子文件夹，每个子文件夹有多组文件，计算spl差。执行文件exert调用函数importdata声压级计算
一个文件夹下多个文件的SPL
执行文件exert
% 功能：遍历data文件夹下的所有文件，计算带A计权的声压级曲线，并在命令行显示SPL（dB）% 编辑者：Heart_sea% 日期：,4,24clear;clc;close all;% ======================= preferences set =================================fs = 12000;% sampling frequency (unit: Hz)NFFT = fs; % the number of samples used in FFTpl = 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; % 灵敏度级dBmic_sen = (10^(MIC_SEN./20)).*OP_AMP;%麦克风的灵敏度% ======================== import data ====================================maindir = './Data';File_L1 =dir(maindir); % 获取路径下的所有文件信息PathNumber=numel(File_L1);% 获取元素的数量m = 1;% 存放数据列数初始化switch timecomcase 1for i = 3:PathNumberPath=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);endtitle(['时域曲线 '])xlabel('t (s)')ylabel('pressure (Pa)')case 2endif timecom == 1figure;end% ======================== 频域 ====================================m = 1; %列放数据初始化for i = 3:PathNumberPath=fullfile(maindir,File_L1(i).name);% 文件夹下文件路径Path2 = fullfile(File_L1(i).name); % 文件夹下的文件的名字Path2(end-3:end) = []; %去掉文件的txtnote(:,m) = {[File_L1(i,1).name]}; % 用作legendm = 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-ASPL = 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. % /help/dsp/ref/fdesign.audioweighting.htmlFilter_Obj = fdesign.audioweighting('WT,Class',Weighting,1,fs);H_Filter = design(Filter_Obj);[Hf,~] = freqz(H_Filter,N_f,fs);for i = 1:N_fSxx(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 == 0Sxx_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 == 1semilogx(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;endif A_weighting_flag == 0 title(['SPL '])elsetitle(['SPLA '])end xlabel('f (Hz)')ylabel('(dB)')endlegend(note)if save_image_flag == 1set(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% 日期：,4,25clear;clc;close all;% ======================= preferences set =================================fs = 12000;% sampling frequency (unit: Hz)NFFT = fs; % the number of samples used in FFTfl = 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; % 灵敏度级dBmic_sen = (10^(MIC_SEN./20)).*OP_AMP;%麦克风的灵敏度% ======================== import data ====================================figure(1)for micNumber = micbegin:micEndfilename = ['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-ASPL_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_fS_off(i) = S_off(i)*((abs(Hf(i)))^2);S_on(i) = S_on(i)*((abs(Hf(i)))^2);endSyy_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 == 0Syy_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 onplot(f,SPL_on)grid onelseif 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'); %降噪前offhold onS_on = csaps(f,SPL_on,z);fnplt(S_on,'r'); % 降噪后oN% set(gca,'YTick',-40:10:50);% set(gca,'XTick',fl:500:fh);grid onendelseif coordinate == 1 % 1代表x对数坐标% ==========是否需要平滑滤波if smooth == 0 subplot (micEnd,1,micNumber)semilogx(f,SPL_off); hold onsemilogx(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 onelseif 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 onS_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 onendendk = SPL_total_off - SPL_total_on;%s是off，c是ondisp(['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)])elsetitle(['SPLA ' num2str(micNumber)])end xlabel('f (Hz)')ylabel('(dB)')end% =============保存图片if save_image_flag == 1set(gca,'FontSize',12)saveas(gcf,'ANCEffectFig.fig') %可以保存成jpg，fig等不同格式的图片end
一个总文件夹下多个子文件夹，每个子文件夹有多组文件，计算spl差。
执行文件exert
## 一个总文件夹下多个子文件夹，每个子文件夹有多组文件### 执行文件exert```java% 功能：遍历data文件夹下的所有子文件夹下的文件（off-on），计算带A计权的每组的声压级曲线，并在命令行显示SPL（dB）% 编辑者：lily% 日期：,4,25clear;clc;close all;% ======================= preferences set =================================maindir = './Data';fs = 12000;% sampling frequency (unit: Hz)NFFT = fs; % the number of samples used in FFTfl = 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; % 灵敏度级dBmic_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 == 1set(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);%前两阶文件夹的名称，字符串格式% ==============设置lengendPath2 = 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) = [];% 将字符串的后四个字符串置成空的，去掉.txtNote_legend(row,i-2) = {str}; row = row+1;end% ======================== 循环dat数据（子文件夹里的文件） ================colum =1; % 存放数据到列初始化v = 2;% subplot第几个fig初始化for m = 3:2:length( dat )% ==============设置dispstr1 = 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-ASPL_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_fS_off(i) = S_off(i)*((abs(Hf(i)))^2);S_on(i) = S_on(i)*((abs(Hf(i)))^2);endSyy_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 == 0Syy_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 onplot(f,SPL_on)grid onelseif 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,''); %降噪前offhold onS_on = csaps(f,SPL_on,z);fnplt(S_on,''); % 降噪后oN% set(gca,'YTick',-40:10:50);% set(gca,'XTick',fl:500:fh);grid onendelseif coordinate == 1 % 1代表x对数坐标if smooth == 0 subplot (max,1,fig)semilogx(f,SPL_off); hold onsemilogx(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 onelseif 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 onS_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 onendendlegend(char(Note_legend(m-2:m-2+1,:)))v = v+1;k = SPL_total_off - SPL_total_on;%s是off，c是ondisp([num2str(Note_disp{fig}) ':' num2str(k) 'dB']) % ============= 标题if A_weighting_flag == 0 title(['SPL ' num2str(fig)])elsetitle(['SPLA ' num2str(fig)])end xlabel('f (Hz)')ylabel('(dB)')endendend
