spectrogram関数を用いたスペクトログラムと、自作のコードでプロットしたスペクトログラムが完全一致しない。
31 次查看(过去 30 天)
显示 更早的评论
①ビルトインされている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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/248939/image.png)
②自作コードによる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)');
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/248940/image.png)
4 个评论
采纳的回答
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 个评论
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!