手动实现MATLAB spectrogram函数：如何保证nfft小于窗长时结果与原函数一致

现存问题 #

• 初版函数在nfft≥窗长场景下，输出仅存在虚部符号相反的差异，其余数值与MATLAB原生函数匹配
• nfft<窗长场景下，输出与MATLAB原生函数完全不匹配

问题根因解析 #

1. 虚部符号相反的本质原因 #

MATLAB原生FFT的计算定义为：
X(k) = sum_{n=0}^{N-1} x(n) * exp(-1j*2*pi*k*n/N)
如果移植时使用的FFT实现采用了正相位的指数项exp(+1j*2*pi*k*n/N)，就会导致所有非DC、非奈奎斯特频点的虚部符号取反，实部保持一致，和测试中观察到的现象完全吻合。

2. nfft<窗长时的MATLAB底层计算逻辑 #

MATLAB原生spectrogram在处理nfft小于窗长的场景时，不会直接截断窗长信号做nfft长度的FFT，实际流程为：

1. 按固定步长对原信号分帧，每帧长度严格等于输入的窗长wlen
2. 对每帧信号施加指定窗函数
3. 计算加窗后完整帧信号的傅里叶变换，在0~2π频率范围内按nfft指定的频率点数量做等间隔采样

对应测试3的复现验证：
测试参数为：窗长wlen=6，nfft=3，分帧第一帧信号为[1,2,3,4,5,6]
① 加hamming窗后做6点FFT，得到的有效频点结果为：9.73 + 0i, -5.2936 + 0.9205i, 0.7279 - 0.3737i, -0.1186 + 0i，与第一次测试的结果完全一致
② 因nfft=3对应的有效频点数量为ceil((3+1)/2)=2，等间隔采样时抽取6点FFT结果的第1个、第3个频点值，最终得到s = [9.73+0i, 0.7279-0.3737i]，与MATLAB原生函数输出完全匹配

初版函数的错误逻辑为：直接取加窗后信号的前nfft个点做nfft长度FFT，本质是对信号时域截断，和MATLAB的频域采样逻辑完全不同，因此结果不匹配。

函数优化实现 #

针对上述根因，对初版函数做两处修正：

修正1：适配FFT符号规则 #

如果存在虚部符号相反的问题，对所有FFT输出做实部保留、虚部取反的操作，匹配MATLAB的相位定义。

修正2：适配nfft<窗长的频域采样逻辑 #

• 当nfft ≥ wlen时，保留原有逻辑：对加窗后的wlen长度信号补零到nfft长度，计算nfft点FFT，直接取前nUnique个有效频点
• 当nfft < wlen时，先计算完整wlen长度加窗信号的FFT，再按等间隔规则抽取对应nfft长度的有效频点，抽取步长为wlen/nfft

优化后代码 #

function out = manualSpectrogram(x, win, noverlap, nfft)
x = x(:);
n = length(x);
wlen = length(win);
nUnique = ceil((1+nfft)/2);
L = fix((n-noverlap)/(wlen-noverlap));
out = zeros(L, nUnique);
index = 1:wlen;

for i = 0:L-1
    xw = win .* x(index);
    if nfft >= wlen
        % nfft≥窗长：补零到nfft长度做FFT
        X = fft(xw, nfft);
        out(i+1, :) = X(1:nUnique);
    else
        % nfft<窗长：计算整窗长FFT，等间隔抽点
        X_wlen = fft(xw, wlen);
        pick_idx = round( (0:nUnique-1) * wlen / nfft ) + 1;
        out(i+1, :) = X_wlen(pick_idx);
    end
    index = index + (wlen - noverlap);
end

% 若存在虚部符号相反的问题，打开下行注释
% out = real(out) - 1j*imag(out);
end

匹配度验证 #

• 测试1（nfft=6=wlen=6）：输出与MATLAB完全一致
• 测试2（nfft=6>wlen=3）：补零逻辑与MATLAB规则一致，输出匹配
• 测试3（nfft=3<wlen=6）：抽点逻辑输出与MATLAB测试结果完全一致

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