企业官网建设流程全解析

简介：这份MATLAB脚本面向从事声学测量、环境噪声评估与工业设备噪声检测的工程师及研究人员，用于对时域噪声信号执行标准A计权处理，并完成1/3倍频程频谱划分与声压级计算。脚本内置ISO 266与IEC 61672推荐的A计权滤波器系数，支持自定义采样率与频带范围，输出各1/3倍频程中心频率对应的A加权声压级，可直接用于合规性评估。资源包共8个文件，约250KB，包含核心脚本、Python对照实现、A计权曲线与NR曲线示意图、依赖说明及配置文件，结构紧凑便于快速上手。已有106人学习下载。读者可获得一套无需额外工具箱、兼容MATLAB R2015b及以上版本的完整分析流程，输入列向量时域信号即可得到结构化频带级数据，方便后续绘图或导出，同时可借助曲线图与对照脚本理解A计权与1/3倍频程的实现细节，适合声学测试数据处理与算法验证。

1. 从一段噪声录音说起：这套 A 计权 1/3 倍频程脚本到底能干什么

手上有一段现场录回来的噪声 wav，领导要你出一张「各频段声压级」的柱状图，还要标出 A 计权总级。你打开 MATLAB，发现手头只有pwelch和fft，算出来的谱是窄带谱，跟声学报告里那种一根一根的 1/3 倍频程柱子完全对不上。这时候大多数人会去翻 IEC 61672 和 GB/T 3785，然后被 A 计权传递函数、频带归一化、参考声压 20 μPa 这一堆东西绕进去。这套 MATLAB 脚本解决的正是这个断层：把一段时域声压信号，直接算成符合声学惯例的 A 计权 1/3 倍频程声压级，输出中心频率、带宽、各带级和总级。

它适合三类人：做 NVH 或环境噪声测试、需要快速出频带报告的工程师；带声学实验课、要给学生一套可复现流程的教师；以及手上有实测数据、想自己搭一套分析链路的开发者。核心链路是「时域采样 → 窄带谱 → 频带能量归并 → A 计权修正 → 声压级换算」，每一步都有明确的物理量纲和参数约定，脚本把这四步封成了可调用的函数。下面按「先立住原理、再动手复现、最后讲坑」的顺序拆开讲，参数怎么设、哪里容易翻车，都会落到具体代码上。

2. A 计权与 1/3 倍频程的物理约定：先把量纲和频带定死

2.1 A 计权到底在修正什么

人耳对低频不敏感、对 2～5 kHz 附近最敏感，A 计权就是把这个主观响度特性做成一条频率响应曲线，叠加到物理声压级上。它的传递函数在 IEC 61672-1 里有解析式，脚本里一般用双线性变换或直接按标准公式算幅度平方：

$$R_A(f) = \frac{12194^2 f^4}{(f^2+20.6^2)\sqrt{(f^2+107.7^2)(f^2+737.9^2)}(f^2+12194^2)}$$

A 计权级就是 $L_A = 10\log_{10}(R_A(f))$，单位 dB。注意这是「修正量」，加在窄带谱的每个频点上，不是加在总级上。很多人图省事，先算总声压级再减一个固定 A 值，这是错的——A 计权是频率相关的，必须逐频点乘进去再积分。脚本里对应的是A_weighting(f)这个函数，输入频率向量，返回同长度的 dB 修正量。

2.2 1/3 倍频程的中心频率与带宽

1/3 倍频程不是随便分的，中心频率按 $f_c = 1000 \cdot 10^{n/10}$ 生成，n 取整数，标准序列从 12.5 Hz 到 20 kHz。每个带的上下限是 $f_u = f_c \cdot 2^{1/6}$、$f_l = f_c / 2^{1/6}$，带宽约为中心频率的 23%。脚本里通常预置一张标准中心频率表，避免自己算错边界。归并时把落在 $[f_l, f_u]$ 内的窄带谱线能量求和，得到该带的均方声压，再转声压级：

$$L_p = 10\log_{10}\left(\frac{p_{rms}^2}{p_0^2}\right), \quad p_0 = 20,\mu Pa$$

这里有个容易忽略的点：窄带谱如果是功率谱密度（PSD），归并时要乘频率分辨率 $\Delta f$ 才是能量；如果是已经积分的功率谱，就直接求和。脚本里用pwelch默认输出 PSD，所以归并前必须乘df，这一步漏了，结果会随 FFT 点数漂移，是血泪经验里最常见的翻车点。

2.3 采样率、窗长与频率分辨率的三角约束

频率分辨率 $\Delta f = f_s / N$，要覆盖 20 kHz 的 1/3 倍频程上限，采样率至少 48 kHz；要让最低带 12.5 Hz 有足够谱线，窗长 N 得够大。常见做法是 $f_s = 48k$、$N = 32768$，此时 $\Delta f \approx 1.46$ Hz，12.5 Hz 带内有约 2 条谱线，勉强够用。如果只关心 100 Hz 以上，$N = 16384$ 也行。脚本把fs和nfft做成入参，就是让你按被测信号的频率范围去权衡，而不是无脑套默认值。

3. 脚本落地：从 wav 到 1/3 倍频程柱状图的完整复现

3.1 读取信号与预加重、去直流

第一步是把 wav 读进来，转成声压。录音设备给的通常是电压或归一化浮点，要乘灵敏度系数（比如 50 mV/Pa 对应的倒数）才回到 Pa。脚本里我一般留一个sensitivity参数，默认 1，让用户自己填。

% 读取噪声录音并转成声压 (Pa) [raw, fs] = audioread('noise_sample.wav'); % raw 为归一化电压 raw = raw(:,1); % 多声道取第一路 sensitivity = 50e-3; % 传声器灵敏度 50 mV/Pa p = raw / sensitivity; % 转成声压 Pa p = p - mean(p); % 去直流，避免 0 Hz 能量污染低频带

逻辑说明：audioread返回的是 -1 到 1 的浮点，代表电压归一化值；除以灵敏度得到 Pa。去直流用减均值，比高通滤波更干净，因为 12.5 Hz 带对直流泄漏很敏感。参数sensitivity必须按实际传声器标定值改，填错的话总级会整体平移，这是量纲层面的错，后面怎么调都救不回来。

3.2 窄带功率谱与 A 计权逐点修正

用pwelch算 PSD，再逐频点乘 A 计权修正量。注意pwelch返回的是单边谱，频率向量从 0 到 $f_s/2$。

% 计算窄带 PSD 并施加 A 计权 nfft = 32768; window = hann(nfft); overlap = nfft/2; [pxx, f] = pwelch(p, window, overlap, nfft, fs); % pxx 单位 Pa^2/Hz df = fs / nfft; % 频率分辨率 A_dB = A_weighting(f); % 逐频点 A 计权修正量 pxx_A = pxx .* 10.^(A_dB/10); % 修正到 A 计权 PSD

逻辑说明：10.^(A_dB/10)是把 dB 修正量转成线性功率因子，乘到 PSD 上。A_weighting内部按 2.1 的公式实现，输入频率向量返回同长度 dB。参数nfft决定df，进而影响 3.3 的归并精度；overlap取 50% 是常规做法，能平滑谱估计但不改变能量总量。

3.3 频带能量归并与声压级换算

把 A 计权 PSD 按 1/3 倍频程边界归并，乘df得带内均方声压，再转 dB。

% 1/3 倍频程归并 fc = [12.5 16 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 16000 20000]; % 标准中心频率 L_band = zeros(size(fc)); for k = 1:length(fc) fl = fc(k) / 2^(1/6); fu = fc(k) * 2^(1/6); idx = (f >= fl) & (f < fu); ms_p = sum(pxx_A(idx)) * df; % 带内均方声压 L_band(k) = 10*log10(ms_p / (20e-6)^2); % 转声压级 dB end L_total = 10*log10(sum(10.^(L_band/10))); % A 计权总级

逻辑说明：sum(pxx_A(idx)) * df是带内能量积分，量纲回到 Pa²；除以参考声压平方再取对数得 dB。L_total用能量求和而非算术平均，因为各带是不相干声源。参数fc表可按需截断，比如只关心 100 Hz～10 kHz，就把两端删掉，但总级要注明覆盖范围，否则和报告对不上。

3.4 出图与结果核对

最后画柱状图，横轴用对数刻度，纵轴声压级，总级标在标题或图例里。

% 绘制 1/3 倍频程柱状图 figure; bar(log10(fc), L_band); set(gca, 'XTick', log10(fc(1:3:end)), ... 'XTickLabel', num2str(fc(1:3:end)', '%g')); xlabel('中心频率 (Hz)'); ylabel('A 计权声压级 (dB)'); title(sprintf('A 计权 1/3 倍频程分析，总级 %.1f dB(A)', L_total)); grid on;

逻辑说明：横轴取log10(fc)让倍频程等距，XTick每隔三个标一次避免拥挤。出图后先看总级是否落在合理范围（正常环境噪声 40～80 dB(A)），再看低频带有没有异常高值——如果有，多半是去直流没做干净或灵敏度填错。

4. 避坑与排查：五个让结果偏掉 10 dB 的常见问题

4.1 现象：总级比参考仪器低 20 dB 左右

原因：灵敏度系数填成了 mV/Pa 的数值但没转成 V/Pa，或者干脆没除。解决：确认sensitivity单位是 V/Pa，50 mV/Pa 要写50e-3，填完用一段已知声压级的校准信号验一遍。

4.2 现象：低频带（12.5～31.5 Hz）数值忽高忽低

原因：nfft太小，df大于 1 Hz，12.5 Hz 带内谱线数不足，归并时能量被截断。解决：把nfft提到 32768 以上，或对信号做零填充；如果信号本身长度不够，就降低最低分析频率，别硬算 12.5 Hz。

4.3 现象：A 计权总级和线性总级差得离谱

原因：把 A 计权修正量加到了总级上，而不是逐频点乘。解决：检查pxx_A是不是在归并前就乘了10.^(A_dB/10)，顺序错了结果必错。

4.4 现象：换一段录音，同一频带结果差 3 dB 以上

原因：pwelch的窗和重叠没固定，或者信号有截断没加窗。解决：把window、overlap、nfft写成固定参数，别用默认值；对非整周期信号，加 Hann 窗并保证重叠 50%。

4.5 现象：总级和商用声级计对不上，差 1～2 dB

原因：商用仪器可能用了不同的频带边界约定或计权实现细节，也可能你的传声器频响没修正。解决：先确认双方都用 IEC 61672 的 A 计权公式，再检查传声器校准证书上的频响修正曲线，必要时在pxx上再乘一个频响修正因子。

5. 进阶技巧：把脚本改成可批量、可验证的分析工具

5.1 批量处理与参数外置

实际项目里不会只分析一段录音。把fs、nfft、sensitivity、fc范围抽到一个config.m或结构体里，主脚本循环读文件夹下的 wav，结果写进表格。这样换设备、换工况只改配置，不动算法。

% 批量分析配置 cfg.fs = 48000; cfg.nfft = 32768; cfg.sensitivity = 50e-3; cfg.fc = [100 125 160 200 250 315 400 500 630 800 1000 ... 1250 1600 2000 2500 3150 4000 5000 6300 8000 10000]; files = dir('data/*.wav'); results = table(); for i = 1:length(files) [p, fs] = read_and_calibrate(fullfile('data', files(i).name), cfg); [L_band, L_total] = third_octave_A(p, fs, cfg); results.File{i} = files(i).name; results.L_total(i) = L_total; results.Bands{i} = L_band; end writetable(results, 'result_summary.csv');

逻辑说明：read_and_calibrate封装读取和灵敏度换算，third_octave_A封装 3.2～3.3 的算法，主循环只负责调度和汇总。参数全走cfg，避免函数里散落魔法数字。这样一套脚本能直接接进测试流程，而不是每次手动改。

5.2 用白噪声和纯音验证算法

算法写完必须验证。生成一段已知 A 计权总级的白噪声，或一个 1 kHz、94 dB 的纯音，跑一遍看输出对不对。1 kHz 纯音在 A 计权下约 -0.0 dB 修正，如果输入是 94 dB 线性，输出应该接近 94 dB(A)；白噪声则要看它的谱密度是否平坦，总级应等于谱密度积分加 A 计权修正。

5.3 一个我常犯的错

早些年我图快，直接把pwelch的 PSD 当功率谱用，归并时没乘df，结果同一段信号换个nfft总级就变，排查了一下午才发现是量纲问题。从那以后我每次改完nfft或fs，都强制用 1 kHz 纯音跑一遍验证，确认总级不随分辨率漂移才继续。这个习惯帮我省掉了无数次返工。希望帮到你。

本文还有配套的精品资源，点击获取