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

①ビルトインされている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)');

[Kazuya]

自作コードの処理内容をもう少しコメント付けて頂けると読み解きやすいと思いまして・・例えば

ydft = fft(y) .* rectwin(nfft)';

の部分、箱型ウィンドウなので影響無いですが、fft 後にウィンドウかけています？

[Kei Manabe]

コメントありがとうございます。少しですが、取り急ぎコメントを付けました。箱型ウィンドウの件ですが、ビルトインspectrogram関数との対比の為に、明示的に書いただけです。ご指摘の通り、この部分は何も仕事をしておりません。

[Yoshio]

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

[Kei Manabe]

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

Answer 1

[Kazuya]

% 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

[Kazuya]

ちなみに spectrogram 関数からパワー スペクトル推定値を返すと数値の直接の比較ができるので、是非どうぞ。

% サンプルデータ
load handel.mat
filename = 'handel.wav';
audiowrite(filename,y,Fs);
%%

[x, fs] = audioread(filename);
nfft = 128;
window = rectwin(nfft);
noverlap = nfft * 0;

% spectrogram でプロット
figure(1)
spectrogram(x, window, noverlap, nfft, fs, 'yaxis');
ha1 = gca;
clim1 = ha1.CLim; % カラーマップ範囲を取得（後ほど利用）

% spectrogram で計算した値をプロット
figure(2)
[~,f,t,ps] = spectrogram(x, window, noverlap, nfft, fs, 'yaxis');
s = 10*log10(ps);
surf(t,f,10*log10(ps),'EdgeColor','none'); 
ha2 = gca;
ha2.CLim = clim1; % カラーマップ表示範囲を同じに設定
axis tight; 
view(0,90); colorbar;

同じ結果になるはず。

[Kei Manabe]

ありがとうございます。遅くなりましたが、いろいろ試してみて、自作コードもビルトイン関数も全く同じ値を得る事が出来ましたので、両者が合わなかった理由を記載します。

①自作コードにおいて、下記右部分にある、２乗とPSD計算部分が抜けていた。

YDFT = abs(ydft).^2/(nfft*fs);

②自作コードにおいて、負の周波数分を正の周波数に組み入れようとして、logを取ってから２倍していた。

しかしそれでも、ビルトイン関数を使った場合、グラフの解像度が荒く、同様のグラフをsurf関数でプロット出来ませんでした。別の言い方をすると、ビルトイン関数でもっと解像度の高いスペクトログラムをプロットしたい場合、どうすればいいのでしょうか？

[Kazuya]

自分が上に書いたコードとデータでは同じように描けているのですが、なんででしょうね。

その違いが出るデータも合わせて、コードを実行すればご指摘の現象(解像度の違い)が再現できる形で記載してもらえますでしょうか？ こちらでも確認して原因を考えてみます。

[Kei Manabe]

下記が違いが出るコードです。データは私のデータをアップして良いか機密的に微妙なので、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の倍数であるようなのですが、この数字に何か意味がありそうです。取り急ぎ、現状をご報告致しました。

[Kei Manabe]

-7.1054e-15は、eps（浮動小数点相対精度）を32倍した値のようですね。

[Kazuya]

カラーマップの表示域を合わせて、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]

ありがとうございます。記載頂いたコードで試してみたところ、私のデータではやはり差がありました。そこで、自作コードにおいて、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（浮動小数点相対精度）の整数倍、両者の間で値が揺らぐ事については未だ謎ですが…。

[Kazuya]

なるほどー、確かに値がそのまま色に反映されそうな image 関数の方が適しているかもしれませんね。勉強になりました。

[Kei Manabe]

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

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

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

[Kazuya]

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

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