WAV文件谱密度分析dB值异常偏低（低约80dB）问题求助

我改编了一段Matlab代码用于WAV文件的谱密度分析，代码运行正常且绘图符合预期，但dB SPL值比预期低约80dB。我推测问题可能与未加入1μPa参考值有关，但无法定位问题根源。作为一名正在学习编程的生物学者，恳请各位提供解决思路。

代码如下：

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FFT parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
NFFT = 96000;    % 
NOVERLAP = round(0.75*NFFT);
w = hanning(NFFT);
% spectrogram dB scale
spectrogram_dB_scale = 100;  % dB range scale (means , the lowest displayed level is XX dB below the max level)
gain = 20
sensitivity = - 165
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% load signal
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[data,Fs]=audioread('sample.wav'); %(newer matlab)
samples = length(data);
dt = 1/Fs;
t = (0:dt:(samples-1)*dt);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% display 1 : averaged FFT spectrum
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[sensor_spectrum, freq] = pwelch(data,w,NOVERLAP,NFFT,Fs);
% convert to dB scale (ref = 1)
sensor_spectrum_dB = 20*log10(sensor_spectrum) - sensitivity - gain ;


figure(1),semilogx(freq,sensor_spectrum_dB);grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(freq(2)-freq(1)) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB (L))');

你的推测方向正确，核心问题是没有正确将WAV归一化数据转换为实际声压，也未结合1μPa参考值计算dB SPL，同时对pwelch输出的功率谱密度（PSD）的转换逻辑有误，具体修正点如下：

1. 关键概念梳理 #

• WAV文件的data是归一化到[-1,1]的数值，对应声卡采集的电压范围（16位WAV通常对应±1V或±3V，若不确定可按通用逻辑处理）。
• 麦克风灵敏度sensitivity=-165dB re 1V/μPa：表示1μPa的声压会产生10^(-165/20) V的输出电压，电压值除以该系数即可得到对应声压（单位Pa）。
• pwelch输出的是功率谱密度（PSD），单位为(幅值²/Hz)，需通过sqrt(PSD * 频率分辨率)转换为RMS幅值。

2. 具体修正步骤 #

步骤1：计算频率分辨率 #

df = freq(2) - freq(1); % 单个频率点的带宽

步骤2：将PSD转换为RMS电压幅值 #

先把pwelch的PSD结果转成RMS电压，同时还原20dB增益的影响：

gain_factor = 10^(gain/20); % 20dB增益对应的放大倍数
sensor_rms_volts = sqrt(sensor_spectrum * df) / gain_factor;

步骤3：转换为实际声压（Pa） #

利用麦克风灵敏度，将电压值转换为声压：

% 1μPa对应的输出电压是10^(sensitivity/20) V，电压除以该值得到声压（Pa）
sensor_rms_pa = sensor_rms_volts / 10^(sensitivity/20);

步骤4：计算dB SPL（参考值为1μPa=1e-6 Pa） #

sensor_spl_dB = 20*log10(sensor_rms_pa / 1e-6);

3. 修正后的核心代码片段 #

替换原代码中sensor_spectrum_dB的计算及绘图部分：

[sensor_spectrum, freq] = pwelch(data,w,NOVERLAP,NFFT,Fs);

% 新增声压转换步骤
df = freq(2) - freq(1);
gain_factor = 10^(gain/20);
sensor_rms_volts = sqrt(sensor_spectrum * df) / gain_factor;
sensor_rms_pa = sensor_rms_volts / 10^(sensitivity/20);
sensor_spl_dB = 20*log10(sensor_rms_pa / 1e-6);

% 调整绘图标签
figure(1),semilogx(freq,sensor_spl_dB);grid
title(['Averaged FFT Spectrum  / Fs = ' num2str(Fs) ' Hz / Delta f = ' num2str(df) ' Hz ']);
xlabel('Frequency (Hz)');ylabel('Amplitude (dB SPL)');

额外说明 #

• 若WAV为多通道文件，需先提取单通道数据（如data = data(:,1);）再处理。
• 若不确定声卡满量程电压，上述修正已覆盖核心偏差来源，足以解决80dB左右的差值问题。

内容的提问来源于stack exchange，提问作者Benjamin Colbert