Matlab自定义pwelch实现与官方函数的增益及频点差异问询

1. 核心问题：分块逻辑错误 #

你的代码中第一个块的构造逻辑完全错误：

if i == 1
    block = [zeros(1, w * (1-overlap)) data(1:w * overlap)] .* hamming_win;

这相当于把第一个块做成了128个零 + 前128个数据，再乘以256点的汉明窗。但官方pwelch的分块逻辑是：第一个块直接取前256个完整数据，后续块每次滑动w*(1-overlap)=128点，取连续的256点数据。这种错误的零填充分块会导致第一个块能量异常，直接拉低整体平均后的PSD增益，同时引入频谱泄漏，造成首尾频点的结果偏差。

正确的分块逻辑应该去掉特殊判断，所有块统一取连续的w点数据：

for i = 1:num_w
    start_idx = 1 + (i-1)*w*(1-overlap);
    end_idx = start_idx + w - 1;
    block = data(start_idx:end_idx) .* hamming_win;
    % 后续FFT和累加逻辑不变
end

你的参数刚好满足len/(w*overlap)为整数（65536/(256*0.5)=512），所以所有块都能完整取到数据，不会越界。

2. 增益差异的根源 #

你手动添加的0.318是巧合的近似值，真正的问题来自两个方面：

（1）错误的第一个块拉低了整体平均 #

第一个块的零填充导致其能量远低于正常块（只有前128个数据有值），使得最终平均结果被大幅压低。手动乘以0.318只是临时抵消了这个错误带来的增益损失，并非正确解法。

（2）单边PSD的归一化缺失 #

对于实信号的PSD计算，官方pwelch默认返回单边PSD：除了直流（第一个点）和Nyquist（最后一个点）分量，其余频点的能量需要乘以2（因为双边谱的能量分布在正负频率，单边谱要合并这两部分能量）。你的代码中没有做这个处理，修正分块逻辑后需要补充该步骤：

Pxx = (block_fd_half .* conj(block_fd_half)) / (w * win_normal);
% 单边PSD修正：中间频点乘以2
Pxx(2:end-1) = Pxx(2:end-1) * 2;

注：官方pwelch默认采样频率Fs=1，所以你的代码无需额外除以Fs，结果单位与官方一致。

3. 首尾频点差异的解决 #

首尾频点的偏差完全来自第一个块的零填充错误。修正分块逻辑后，所有块都是完整的有效数据，频谱泄漏消失，首尾频点结果会和官方pwelch完全对齐。

修正后的完整代码 #

len = 65536;
data = randn(1, len);
overlap = 0.5;
w = 256;
num_w = len / (w * overlap); % 非整数场景需调整为floor((len - w)/(w*(1-overlap)))+1
hamming_win = hamming(w)';
win_normal = mean(hamming_win.^2);
Pxx_sum = zeros(1, w/2+1);

for i = 1:num_w
    start_idx = 1 + (i-1)*w*(1-overlap);
    end_idx = start_idx + w - 1;
    block = data(start_idx:end_idx) .* hamming_win;
    
    block_fd = fft(block);
    block_fd_half = block_fd(1:w/2+1);
    Pxx = (block_fd_half .* conj(block_fd_half)) / (w * win_normal);
    % 单边PSD修正
    Pxx(2:end-1) = Pxx(2:end-1) * 2;
    
    Pxx_sum = Pxx_sum + Pxx;
end

psd_avg = Pxx_sum / num_w;
figure, plot(psd_avg);

% 官方函数对比
[psd_mat, freq] = pwelch(data, w, 0.5*w);
figure, plot(psd_mat);

验证说明 #

修正后的代码运行后，psd_avg与psd_mat的结果会几乎完全重合，无需再手动添加增益系数，首尾频点的差异也会彻底消失。

内容的提问来源于stack exchange，提问作者Ben Li