maMATLABにおいて、spectrogram関数を用いたスペクトログラムと、自作のコードでプロットしたスペクトログラムが完全一致しない。

Question

MATLABでの話です。

①ビルトインされているspectrogram関数
②自作コードによるspectrogram

上記２パターンで、プロットしたスペクトログラムが一致しません。
ビルトインの方が、見た目が荒く、周波数分解能が高く見えます。
（両者の値そのものが異なっている（ビルトインの方は-150dB ~ -50dB、自作の方は-80dB ~ +10dB）事も気にはなっているのですが、直接は関係ないと認識しております。）

どなたか、原因を見つけられる方はおりませんでしょうか？
よろしくお願いいたします。

①ビルトインされているspectrogram関数

[x, fs] = audioread(filename);
nfft = 128;
window = rectwin(nfft);
noverlap = nfft * 0;
spectrogram(x, window, noverlap, nfft, fs, 'yaxis');
colormap hot

②自作コードによるspectrogram

[x, fs] = audioread(filename);

sample_length = length(x)/fs;
nfft = 128;
window = rectwin(nfft);
noverlap = nfft * 0;
J = fix((length(x)-noverlap) / (nfft-noverlap));

S = zeros(nfft/2+1, J);

for jj = 1 : J
    y = x(1 + (jj-1) * nfft : jj * nfft)';

　　% 矩形窓をかけている部分は形式的なものなので特に何もしていません。
    ydft = fft(y) .* rectwin(nfft)';
    
    % 絶対値のlogを取っています。
    YDFT = abs(ydft);
    YDFT = 10*log10(YDFT);
    
    % FFTの出力が、正負の周波数に分布しているので、
    % 正の周波数と0[Hz]成分（0[Hz]~64[Hz]）のデータだけ取り出しています。
    YDFT = YDFT(1:length(ydft)/2+1);
    
    % 取り残された負の周波数（-64[Hz]~-1[Hz]）のエネルギーを、
    % 取り出した正の周波数に組み入れる為に、２倍しています。
    YDFT(2:end-1) = 2 * YDFT(2:end-1);
    
    S(:, jj) = YDFT;
end

% スペクトログラムをsurf関数でプロットする為に、周波数軸のグリッドを作っています。
F = 0 : (fs/2)/(nfft/2) : fs/2;

% スペクトログラムをsurf関数でプロットする為に、時間軸のグリッドを作っています。
T = sample_length/J : sample_length/J : sample_length;

surf(T, F, S, 'EdgeColor', 'flat', 'FaceColor', 'interp', 'FaceLighting', 'none')
axis xy; axis tight; colormap(hot); view(2); colorbar;

c.Label.String = 'Power / frequency [dB / Hz]';
xlabel('Time');
ylabel('Frequency (Hz)');

Answer 1

とりあえず原因が分かりましたので記載いたします。

①自作コードにおいて、下記右部分にある、２乗とPSD計算部分が抜けていた。
YDFT = abs(ydft).^2/(nfft*fs);

つまり、ビルトイン関数では、PSDを計算していたのに、自作コードでは、PSDになっていなかったという事です。

②自作コードにおいて、負の周波数分を正の周波数に組み入れようとして、logを取ってから２倍していた。

両者で数字上は、全く同じになりましたが、得られたプロットの解像度は、ビルトイン関数の方が悪く見えます。この件に関しては、MATHWORKS社に聞いてみたいと思います。

ありがとうございます。

Answer 2

for jj = 1 : J
    y = x(1 + (jj-1) * nfft : jj * nfft)';
    ydft = fft(y(:).*rectwin(nfft));
    S(:, jj) = ydft(1:length(ydft)/2+1);
end

と

P=spectrogram(x, window, noverlap, nfft, fs, 'yaxis');

を

 max(abs(P(:)-S(:)))

により比較してみると、同じものを計算していることがわかると思います。
違いはプロットの方法です。その点についてはご自身で調べてください。