スペクトログラム表示におけるよくある誤解と修正点

MATLABのspectrogram関数を用いて短時間フーリエ変換（STFT）を実行する際、出力されるグラフの時間軸・周波数軸が物理単位と合わず、FFTの結果との対応が取りにくい場合がある。これは主に、サンプリング周波数の未指定と、window引数の誤解に起因する。

以下の例では、周波数が時間とともに変化する信号を生成し、誤った方法と正しい方法でスペクトログラムを求める。

問題のある実装例

samplingRate = 50; % サンプリング周波数
timeVec = 0:(1/samplingRate):5; % 時間ベクトル
sig = [sin(2*pi*8*timeVec(timeVec <= 2.5)), sin(2*pi*20*timeVec(timeVec > 2.5))]; % 8Hzから20Hzへ遷移する信号

numSamples = length(sig);
fftSize = 2^nextpow2(numSamples);
sigFFT = fft(sig, fftSize) / numSamples;
freqAxis = samplingRate/2 * linspace(0, 1, fftSize/2 + 1);

figure;
subplot(3,1,1);
plot(timeVec, sig, 'k');
xlabel('Time (s)'); ylabel('Amplitude');

subplot(3,1,2);
plot(freqAxis, 2*abs(sigFFT(1:fftSize/2+1)));
xlabel('Frequency (Hz)'); ylabel('Amplitude');

subplot(3,1,3);
% 誤ったウィンドウ指定：時間ベクトルの部分集合を渡している
win = timeVec(1:20:end); 
[S, F, T] = spectrogram(sig, win);
pcolor(T, F, abs(S)); shading interp;
xlabel('Time (s)'); ylabel('Frequency (Hz)');

上記のコードでは以下の2つの問題が発生する。

1. spectrogramの呼び出し時にサンプリング周波数を指定していないため、出力されるTとFが正規化された無次元の値となり、物理スケール（秒・ヘルツ）と一致しない。
2. window引数に時間ベクトルの間引きデータを渡しているが、この引数はウィンドウを適用する時間的位置を指定するものではない。ここにはウィンドウ関数の係数ベクトル、あるいはウィンドウ長（スカラ値）を指定する必要がある。

修正後の実装例

% ウィンドウ長をスカラ値で指定（この場合、既定のハミング窓が適用される）
winLen = 80;

% 第5引数にsamplingRateを指定することで、TとFのスケールが正確になる
[S, F, T, P] = spectrogram(sig, winLen, [], [], samplingRate);

% 代替案：矩形窓を明示的に定義する場合
% winCoeff = rectwin(winLen); 
% [S, F, T, P] = spectrogram(sig, winCoeff, [], [], samplingRate);

figure;
subplot(3,1,1);
plot(timeVec, sig, 'k');
xlabel('Time (s)'); ylabel('Amplitude');

subplot(3,1,2);
plot(freqAxis, 2*abs(sigFFT(1:fftSize/2+1)));
xlabel('Frequency (Hz)'); ylabel('Amplitude');

subplot(3,1,3);
% パワースペクトル密度をデシベルスケールで表示
imagesc(T, F, 10*log10(abs(P)));
axis xy; axis tight;
colormap(jet);
xlabel('Time (s)');
ylabel('Frequency (Hz)');

修正のポイントとして、スペクトログラムの縦軸は周波数を表すのに対し、通常のFFT振幅スペクトルの縦軸は振幅を表すことに注意する。両者のY軸を直接比較することはできないが、スペクトログラムのカラー強度（Z軸）がFFTの縦軸の振幅情報に相当する。