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

Question

Kei Manabe 2019-11-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/491806-spectrogram

评论： Kazuya 2019-11-22

采纳的回答： Kazuya

在 MATLAB Online 中打开

①ビルトインされているspectrogram関数

②自作コードによるspectrogram

上記２パターンで、プロットしたスペクトログラムが一致しません。

ビルトインの方が、見た目が荒く、周波数分解能が高く見えます。

（両者の値そのものが異なっている（ビルトインの方は-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]）のエネルギーを、
    % 取り出した正の周波数に組み入れる為に、２倍しています。
    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 个评论
显示 2更早的评论隐藏 2更早的评论

Yoshio 2019-11-19

编辑：Yoshio 2019-11-19

1s = spectrogram()として、sの値とご自身の結果を突き合わせてみてはどうでしょうか。チェックポイントとしては、時系列データの二乗和とlogを取らない周波数成分の二乗和が一致するかも確認してみましょう。またsurfではなく, imageで両者をプロットしてみてはどうでしょうか。

Kei Manabe 2019-11-20

ありがとうございます。両者で全く同じ値を得る事には成功しました。しかし、ビルトイン関数を使った場合、グラフの解像度が荒く、同様のグラフをsurf関数でプロット出来ませんでした。別の言い方をすると、ビルトイン関数でもっと解像度の高いスペクトログラムをプロットしたい場合、どうすればいいのでしょうか？

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Kazuya 2019-11-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/491806-spectrogram#answer_402331

编辑：Kazuya 2019-11-20

在 MATLAB Online 中打开

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

10 个评论
显示 8更早的评论隐藏 8更早的评论

Kei Manabe 2019-11-21

在 MATLAB Online 中打开

下記が違いが出るコードです。データは私のデータをアップして良いか機密的に微妙なので、handel.wavで試してみました所、こちらでも同様の解像度の違いが再現されました。

format compact
format long
% load data
filename = 'handel.wav';
[~, name, ~] = fileparts(filename);
[x, fs] = audioread(filename);
sample_length = length(x) / fs;
nfft = 128;
window = rectwin(nfft);
noverlap = nfft * 0;
spectrogram(x, window, noverlap, nfft, fs, 'yaxis');
% To compare with my own code
[~,~,~,spe] = spectrogram(x, window, noverlap, nfft, fs, 'yaxis');
colormap hot
J = fix((length(x)-noverlap) / (nfft-noverlap));
figure;
S = zeros(nfft/2+1, J);
for jj = 1 : J
    y = x(1 + (jj-1) * nfft : jj * nfft)';
    
    ydft = fft(y) .* rectwin(nfft)';
    
    % Calc of PSD 
    YDFT = abs(ydft).^2/(nfft*fs);
    % Extraction of positive frequency region for a lower frequency data
    YDFT = YDFT(1:length(ydft)/2+1);
    
    % not needed
    YDFT(2:end-1) = 2 * YDFT(2:end-1);
    
    YDFT = 10*log10(YDFT);
    S(:, jj) = YDFT;
end
F = 0 : (fs/2)/(nfft/2) : fs/2;
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)');

ここで、気付いた事があるのですが、両者の数字は厳密には全く同じではありませんでした。MATLAB上で配列変数をダブルクリックして中身を見ると、表示されている数字程度の精度では完全一致しているのですが、e-14オーダーの差がある事が新たに判明しました。

A = 10*log10(spe) - S;

として、Aの中身を見ますと、-7.1054e-15や、-1.4211e-14と言った数字がまばらに分布しております。見たところ、全て-7.1054e-15の倍数であるようなのですが、この数字に何か意味がありそうです。取り急ぎ、現状をご報告致しました。

Kazuya 2019-11-21

在 MATLAB Online 中打开

カラーマップの表示域を合わせて、surf 関数の 'EdgeColor' を 'none' にして表示してみました。

spectrogram 関数から

Surf 関数から

Surf 関数の表示の方が、若干表示に内挿がかかっているようにも見えますが、「解像度」とはこの違いのことをおっしゃっていますか？ 'FaceColor' を 'interp' にした状態だと、まさに内挿した表示になりますので、'none' に変更したほうがいいですね。それでも多少「内挿した感」が残っている気がするのは謎ですね・・。許容範囲内？

以下コード。

format compact
format long
% load data
filename = 'handel.wav';
[~, name, ~] = fileparts(filename);
[x, fs] = audioread(filename);
sample_length = length(x) / fs;
nfft = 128;
window = rectwin(nfft);
noverlap = nfft * 0;
figure(1)
spectrogram(x, window, noverlap, nfft, fs, 'yaxis');
ha1 = gca;
clim1 = ha1.CLim; % カラーマップ範囲を取得（後ほど利用）
% To compare with my own code
[~,~,~,spe] = spectrogram(x, window, noverlap, nfft, fs, 'yaxis');
colormap hot
J = fix((length(x)-noverlap) / (nfft-noverlap));
figure(2);
S = zeros(nfft/2+1, J);
for jj = 1 : J
    y = x(1 + (jj-1) * nfft : jj * nfft)';
    
    ydft = fft(y) .* rectwin(nfft)';
    
    % Calc of PSD 
    YDFT = abs(ydft).^2/(nfft*fs);
    % Extraction of positive frequency region for a lower frequency data
    YDFT = YDFT(1:length(ydft)/2+1);
    
    % not needed
    YDFT(2:end-1) = 2 * YDFT(2:end-1);
    
    YDFT = 10*log10(YDFT);
    S(:, jj) = YDFT;
end
F = 0 : (fs/2)/(nfft/2) : fs/2;
T = sample_length/J : sample_length/J : sample_length;
surf(T, F, S,'EdgeColor','none'); 
axis xy; axis tight; colormap(hot); view(2); colorbar;
c.Label.String = 'Power / frequency [dB / Hz]';
ha2 = gca;
ha2.CLim = clim1; % カラーマップ表示範囲を同じに設定
xlabel('Time');
ylabel('Frequency (Hz)');

Kei Manabe 2019-11-21

在 MATLAB Online 中打开

ありがとうございます。記載頂いたコードで試してみたところ、私のデータではやはり差がありました。そこで、自作コードにおいて、surfではなく、imageでプロットしてみたところ、差が無くなりました！

format compact
format long
% load data
filename = 'handel.wav';
[~, name, ~] = fileparts(filename);
[x, fs] = audioread(filename);
% 自作コードでは、128の窓長で11812回FFTした後に端数は無視していたので、
% ビルトイン関数でも端数を同様に処理しない様に修正しました。
x = x(1:1511936);　
sample_length = length(x) / fs;
nfft = 128;
window = rectwin(nfft);
noverlap = nfft * 0;
spectrogram(x, window, noverlap, nfft, fs, 'yaxis');
% To compare with my own code
[~,~,~,spe] = spectrogram(x, window, noverlap, nfft, fs, 'yaxis');
colormap hot
J = fix((length(x)-noverlap) / (nfft-noverlap));
figure;
S = zeros(nfft/2+1, J);
for jj = 1 : J
    y = x(1 + (jj-1) * nfft : jj * nfft)';
    
    ydft = fft(y) .* rectwin(nfft)';
    
    % Calc of PSD 
    YDFT = abs(ydft).^2/(nfft*fs);
    % Extraction of positive frequency region for a lower frequency data
    YDFT = YDFT(1:length(ydft)/2+1);
    
    % not needed
    YDFT(2:end-1) = 2 * YDFT(2:end-1);
    
    YDFT = 10*log10(YDFT);
    S(:, jj) = YDFT;
end
F = 0 : (fs/2)/(nfft/2) : fs/2;
T = sample_length/J : sample_length/J : sample_length;
image(T, F, S, 'CDataMapping','scaled')
axis xy; axis tight; colormap(hot); view(2); colorbar;
% caxis([-150, -50])
c.Label.String = 'Power / frequency [dB / Hz]';
xlabel('Time');
ylabel('Frequency (Hz)');

ビルトイン関数

自作コード

つまり、結論としては、MATLABのspectrogram関数は、surf関数ではなく、image関数を使っているっぽいという事です。これで当初の疑問は解決しました。ただ、eps（浮動小数点相対精度）の整数倍、両者の間で値が揺らぐ事については未だ謎ですが…。

Kei Manabe 2019-11-22

新たに気付いた事があるのですが、ビルトイン関数のプロットにおいて、3D回転のアイコンをクリックすると、surfを使った場合の自作コードと同じ精細な見た目に変化するようです。それはまるで、ビルトイン関数において、最初はimageでプロットされているのに、3D回転させようとすると、surfプロットに切り替わっているかの様です。

いずれにせよ、自作コードでもほぼ同じスペクトログラムを得る事が出来たので、この質問は閉じます。また改めて、surf関数、image関数でなぜ見た目が異なるのかについてと、eps整数倍の揺らぎについて、質問を投稿したいと考えております。

数日間、ご協力頂いて、大変助かりました。ありがとうございます。

Kazuya 2019-11-22

いえいえ、ほぼ自己解決されたようなものですし。

こちらこそありがとうございました。

请先登录，再进行评论。

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

采纳的回答

10 个评论
显示 8更早的评论隐藏 8更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

4 个评论 显示 2更早的评论隐藏 2更早的评论

采纳的回答

10 个评论 显示 8更早的评论隐藏 8更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

10 个评论
显示 8更早的评论隐藏 8更早的评论