MATLAB Answers

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

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

댓글을 달려면 로그인하십시오.

채택된 답변

Kazuya
Kazuya 2019년 11월 20일
편집: Kazuya 2019년 11월 20일
% 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
Kazuya
Kazuya 2019년 11월 22일
いえいえ、ほぼ自己解決されたようなものですし。
こちらこそありがとうございました。

댓글을 달려면 로그인하십시오.

추가 답변(0개)

제품


릴리스

R2019a

Community Treasure Hunt

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

Start Hunting!