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

67 次查看(过去 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 个评论
Yoshio
Yoshio 2019-11-19
编辑:Yoshio 2019-11-19
1s = spectrogram()として、sの値とご自身の結果を突き合わせてみてはどうでしょうか。チェックポイントとしては、時系列データの二乗和とlogを取らない周波数成分の二乗和が一致するかも確認してみましょう。またsurfではなく, imageで両者をプロットしてみてはどうでしょうか。
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 个评论
Kei Manabe
Kei Manabe 2019-11-22
新たに気付いた事があるのですが、ビルトイン関数のプロットにおいて、3D回転のアイコンをクリックすると、surfを使った場合の自作コードと同じ精細な見た目に変化するようです。それはまるで、ビルトイン関数において、最初はimageでプロットされているのに、3D回転させようとすると、surfプロットに切り替わっているかの様です。
いずれにせよ、自作コードでもほぼ同じスペクトログラムを得る事が出来たので、この質問は閉じます。また改めて、surf関数、image関数でなぜ見た目が異なるのかについてと、eps整数倍の揺らぎについて、質問を投稿したいと考えております。
数日間、ご協力頂いて、大変助かりました。ありがとうございます。
Kazuya
Kazuya 2019-11-22
いえいえ、ほぼ自己解決されたようなものですし。
こちらこそありがとうございました。

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 データ分布プロット 的更多信息

产品


版本

R2019a

Community Treasure Hunt

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

Start Hunting!