# spectrogra​m関数を用いたスペク​トログラムと、自作の​コードでプロットした​スペクトログラムが完​全一致しない。

119 views (last 30 days)
Kei Manabe on 19 Nov 2019
Commented: Kazuya on 22 Nov 2019
①ビルトインされているspectrogram関数
②自作コードによるspectrogram

ビルトインの方が、見た目が荒く、周波数分解能が高く見えます。
（両者の値そのものが異なっている（ビルトインの方は-150dB ~ -50dB、自作の方は-80dB ~ +10dB）事も気にはなっているのですが、直接は関係ないと認識しております。）

よろしくお願いいたします。
①ビルトインされているspectrogram関数
nfft = 128;
window = rectwin(nfft);
noverlap = nfft * 0;
spectrogram(x, window, noverlap, nfft, fs, 'yaxis');
colormap hot ②自作コードによるspectrogram
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)'); Show 1 older comment
Kei Manabe on 19 Nov 2019
コメントありがとうございます。少しですが、取り急ぎコメントを付けました。箱型ウィンドウの件ですが、ビルトインspectrogram関数との対比の為に、明示的に書いただけです。ご指摘の通り、この部分は何も仕事をしておりません。
Yoshio on 19 Nov 2019
1s = spectrogram()として、sの値とご自身の結果を突き合わせてみてはどうでしょうか。チェックポイントとしては、時系列データの二乗和とlogを取らない周波数成分の二乗和が一致するかも確認してみましょう。またsurfではなく, imageで両者をプロットしてみてはどうでしょうか。
Kei Manabe on 20 Nov 2019
ありがとうございます。両者で全く同じ値を得る事には成功しました。しかし、ビルトイン関数を使った場合、グラフの解像度が荒く、同様のグラフをsurf関数でプロット出来ませんでした。別の言い方をすると、ビルトイン関数でもっと解像度の高いスペクトログラムをプロットしたい場合、どうすればいいのでしょうか？

Kazuya on 20 Nov 2019
Edited: Kazuya on 20 Nov 2019
% modified の部分確認してみてください。log 取るタイミングも重要です。
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);
YDFT = abs(ydft).^2/(nfft*fs); % modified, Kazuya
% 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);
YDFT = 10*log10(YDFT); % modified, Kazuya
S(:, jj) = YDFT;
end

Kazuya on 21 Nov 2019
なるほどー、確かに値がそのまま色に反映されそうな image 関数の方が適しているかもしれませんね。勉強になりました。
Kei Manabe on 22 Nov 2019

いずれにせよ、自作コードでもほぼ同じスペクトログラムを得る事が出来たので、この質問は閉じます。また改めて、surf関数、image関数でなぜ見た目が異なるのかについてと、eps整数倍の揺らぎについて、質問を投稿したいと考えております。

Kazuya on 22 Nov 2019
いえいえ、ほぼ自己解決されたようなものですし。
こちらこそありがとうございました。