MATLAB Answers

0

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

Asked by Kei Manabe on 19 Nov 2019 at 1:01
Latest activity Commented on by Kazuya
on 22 Nov 2019 at 3:41
①ビルトインされているspectrogram関数
②自作コードによるspectrogram
上記2パターンで、プロットしたスペクトログラムが一致しません。
ビルトインの方が、見た目が荒く、周波数分解能が高く見えます。
(両者の値そのものが異なっている(ビルトインの方は-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])のエネルギーを、
% 取り出した正の周波数に組み入れる為に、2倍しています。
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)');

  4 Comments

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

Sign in to comment.

Products


Release

R2019a

1 Answer

Kazuya
Answer by Kazuya
on 20 Nov 2019 at 1:19
Edited by Kazuya
on 20 Nov 2019 at 1:19
 Accepted Answer

% 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])のエネルギーを、
% 取り出した正の周波数に組み入れる為に、2倍しています。
YDFT(2:end-1) = 2 * YDFT(2:end-1);
YDFT = 10*log10(YDFT); % modified, Kazuya
S(:, jj) = YDFT;
end

  10 Comments

Kazuya
on 21 Nov 2019 at 19:17
なるほどー、確かに値がそのまま色に反映されそうな image 関数の方が適しているかもしれませんね。勉強になりました。
Kei Manabe on 22 Nov 2019 at 3:05
新たに気付いた事があるのですが、ビルトイン関数のプロットにおいて、3D回転のアイコンをクリックすると、surfを使った場合の自作コードと同じ精細な見た目に変化するようです。それはまるで、ビルトイン関数において、最初はimageでプロットされているのに、3D回転させようとすると、surfプロットに切り替わっているかの様です。
いずれにせよ、自作コードでもほぼ同じスペクトログラムを得る事が出来たので、この質問は閉じます。また改めて、surf関数、image関数でなぜ見た目が異なるのかについてと、eps整数倍の揺らぎについて、質問を投稿したいと考えております。
数日間、ご協力頂いて、大変助かりました。ありがとうございます。
Kazuya
on 22 Nov 2019 at 3:41
いえいえ、ほぼ自己解決されたようなものですし。
こちらこそありがとうございました。

Sign in to comment.