2022/03/29 追記：重大なバグを発見した．
[~, freq, time, ps] = spectrogram(z,Nw,Nol,NFFT,Fs,'yaxis');
とあるが，psは確率密度関数なのでspectrogramの結果を使用していない．
[s, freq, time, ~] = spectrogram(z,Nw,Nol,NFFT,Fs,'yaxis');
としてpsのかわりにsの値を使用すると良いかと思われる．
後ほど修正予定．
申し訳ないが取り急ぎ報告のみ．
スペクトログラムは横軸時間,縦軸周波数,レベルを色の濃淡で表したグラフである.
任意の周波数で横方向に切り取ったグラフおよび任意の時間で縦方向に切り取ったグラフを描画してみる.
図解するとこんな感じ.
今回はこれらを関数として用意し,mainで呼び出すことにした.
main
% 分析信号としてチャープ信号に白色ノイズを重畳したものを用いる Fs = 44100; % サンプリング周波数(Hz) T = 10; % 信号の時間長さ（秒） N = Fs * T; % 信号のサンプル数 n = randn(1, N); % N点の白色ノイズ % 線形チャープ信号の生成 f1 = 5000; % 開始周波数 f2 = 15000; % 終了周波数 t = ((1:N) - 1) / Fs; % 周期1/Fs, サンプル数Nの時間配列 s = chirp(t,5000,T,15000); % チャープ信号の生成 % ノイズとチャープ信号を足し合わせる SNRdB_original = snr(s,n); % 元のSN比を求める SNRdB = -10; % 設定するSN比 coef = 10^((SNRdB_original - SNRdB)/20); % SN比を調節するための係数 n = coef * n; % ノイズパワーを調整 z = s + n; % チャープ信号とノイズを重畳する % STFTする Nw = 512; % 時間窓の長さ Nol = round(Nw/2); % オーバーラップ NFFT = Nw; % FFT点数 [~, freq, time, ps] = spectrogram(z,Nw,Nol,NFFT,Fs,'yaxis'); % 周波数-レベルグラフを描く t = 2; figure; plot_freq_level_plane_of_spectrogram(ps, Fs, t); % 時間-レベルグラフを描く f = 10*10^3; figure; plot_time_level_plane_of_spectrogram(ps, Fs, f);
plot_freq_level_plane_of_spectrogram()
function [] = plot_freq_level_plane_of_spectrogram(ps,Fs,t) % スペクトログラムの時間t(秒)における周波数-レベル平面を描画する % 周波数配列の計算 N = size(ps,1); % 周波数方向のサイズ freq = ((1:N)-1)*Fs/N/2/1000; % 周波数配列(kHz単位) % 指定時間（秒）における時間-レベルの配列を取得 d_freq = Fs / N; % 周波数分解能 d_time = 1 / d_freq; % 時間分解能 time = round(t / d_time); % 時間軸インデックスを計算 P = ps(:,time); % 配列のコピー % グラフを描画する P = 10*log10(P); % 振幅のdBへの変換 plot(freq,P); % グラフの描画設定 xlabel('周波数（kHz）'); % X軸ラベル xlim([0, Fs/2/1000]); % X軸範囲(ナイキスト周波数(kHz)まで) ylabel('パワー/周波数(dB/Hz)') % Y軸ラベル ylim([-60, 0]); % Y軸範囲 str = sprintf('周波数-レベルグラフ(t=%0.1f秒)',t); title(str); % タイトル end
plot_time_level_plane_of_spectrogram()
function [] = plot_time_level_plane_of_spectrogram(ps,Fs,f) % スペクトログラムの周波数f(Hz)における時間-レベル平面を描画する % 時間軸配列を計算 N = size(ps,1); % 周波数軸方向のサイズ d_freq = Fs / N; % 周波数分解能 d_time = 1 / d_freq; % 時間分解能 N = size(ps,2); % 時間軸方向のサイズ time = ((1:N)-1) * d_time; % 時間軸配列 % 指定周波数（Hz）における時間-レベルの配列を取得 freq = round(f / d_freq * 2); % 周波数インデックスを計算 P = ps(freq,:); % 配列をコピーする % グラフを描画する P = 10*log10(P); % 振幅のdBへの変換 plot(time,P); % グラフの描画設定 xlabel('時間(secs)'); % X軸ラベル xlim([0, (N-1)*d_time]); % X軸範囲 ylabel('パワー/周波数(dB/Hz)') % Y軸ラベル ylim([-60, 0]); % Y軸範囲 str = sprintf('時間-レベルグラフ(f=%0.1fkHz)',f/1000); title(str); % タイトル end
