MATLAB Answers

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

66 views (last 30 days)
Kei Manabe
Kei Manabe on 19 Nov 2019
Commented: Kazuya on 22 Nov 2019
①ビルトインされている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
Kei Manabe
Kei Manabe on 20 Nov 2019
ありがとうございます。両者で全く同じ値を得る事には成功しました。しかし、ビルトイン関数を使った場合、グラフの解像度が荒く、同様のグラフをsurf関数でプロット出来ませんでした。別の言い方をすると、ビルトイン関数でもっと解像度の高いスペクトログラムをプロットしたい場合、どうすればいいのでしょうか?

Sign in to comment.

Accepted Answer

Kazuya
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])のエネルギーを、
% 取り出した正の周波数に組み入れる為に、2倍しています。
YDFT(2:end-1) = 2 * YDFT(2:end-1);
YDFT = 10*log10(YDFT); % modified, Kazuya
S(:, jj) = YDFT;
end
  10 Comments
Kazuya
Kazuya on 22 Nov 2019
いえいえ、ほぼ自己解決されたようなものですし。
こちらこそありがとうございました。

Sign in to comment.

More Answers (0)

Products


Release

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!