matlabにはspectrogramという関数があり,これを使うとスペクトログラムが描ける.
今回はmesh関数で同じように描けるか試してみる.
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の値を使用すると良いかと思われる．
文中で「色合いが若干違うのはなんだろ・・・計算等や使っている変数が違うかもしれないので,要調査.」と書いているが，この辺のバグのせいかと思われる．
後ほど修正予定．
申し訳ないが取り急ぎ報告のみ．
準備
% 分析する信号zを生成する Fs = 44100; % サンプリング周波数(Hz) T = 10; % 信号の時間長さ（秒） N = Fs * T; % 信号のサンプル数 n = randn(1, N); % N点の白色ノイズ % 線形チャープ信号の生成 f1 = 5000; % 開始周波数(Hz) f2 = 15000; % 終了周波数(Hz) t = ((1:N) - 1) / Fs; % 周期1/Fs, サンプル数Nの時間配列 s = chirp(t,f1,T,f2); % チャープ信号の生成 % ノイズとチャープ信号を足し合わせる SNRdB_original = snr(s,n); % 元のSN比を求める SNRdB = -10; % 設定するSN比 coef = 10^((SNRdB_original - SNRdB)/20); % SN比を調節するための係数 n = coef * n; % ノイズパワーを調整 z = s + n; % チャープ信号とノイズを重畳する
分析する信号はチャープ信号に白色ノイズを重畳したものにした.
spectrogramで描画
% STFTする Nw = 512; % 時間窓の長さ Nol = round(Nw/2); % オーバーラップ NFFT = Nw; % FFT点数 spectrogram(z,Nw,Nol,NFFT,Fs,'yaxis'); % スペクトログラムを描画 caxis([-60 0]); % レベルの範囲を指定 colormap jet; % カラーマップを設定
fig.spectrogramで描画したグラフ
これをmeshで再現してみる.
meshを使って描画(修正前)
% meshを使って描画する [~, freq, time, ps] = spectrogram(z,Nw,Nol,NFFT,Fs,'yaxis'); psdB = 10*log10(ps); % dBに変換 figure; % figureウィンドウの作成 mesh(time, freq, psdB); % パワースペクトルを描画 caxis([-60 0]); % レベルの範囲を指定 colormap jet; % カラーマップを設定
fig.meshをそのまま使ったグラフ
視点や軸ラベルが設定されていないので,これを設定する.
meshにいろいろ設定する
% meshを使って描画する（修正後） [~, freq, time, ps] = spectrogram(z,Nw,Nol,NFFT,Fs,'yaxis'); psdB = 10*log10(ps); % dBに変換 figure; % figureウィンドウの作成 freq = freq / 1000; % HzをkHzに変換 mesh(time, freq, psdB); % パワースペクトルを描画 caxis([-60 0]); % レベルの範囲を指定 colormap jet; % カラーマップを設定 view(0,90); % 見る角度を指定 xlabel('時間(secs)'); % X軸ラベル ylim([0,Fs/2/1000]); % Y軸範囲を指定(ナイキスト周波数(kHz)まで) ylabel('周波数(kHz)'); % Y軸ラベル h = colorbar; % カラーバーを表示 ylabel(h, 'パワー/周波数(dB/Hz)');% カラーバーのラベル
fig.meshで描画したグラフ(いろいろ設定後)
結果
上で貼ったグラフを再度貼り付ける.
fig.spectrogramで描画したグラフ
fig.meshで描画したグラフ(いろいろ設定後)
ほぼ同一のグラフを描くことができた.
色合いが若干違うのはなんだろ・・・
計算等や使っている変数が違うかもしれないので,要調査.
