前回はコンスタレーション・プロットで「星座が見えた」と喜びました。今回はパワースペクトル密度(PSD)プロットです。実はScilab “comm_tbx”を使い始めた第75回から既にお世話になっていたもの。要はFFTプロットにお化粧を施しただけのユーティリティなのですが便利なものです。スペアナみたいなグラフが出てきていい感じ？
※「手習ひデジタル信号処理」投稿順 Indexはこちら
※動作確認に使用させていただいているのは、Scilab 6.1.1(Windows版）です。
PSD(Power spectral density)
パワースペクトル密度は、単位周波数あたりのパワースペクトルであると。なんだかな～。御専門の小野測器殿の以下のホームページに解説と計算式などが紹介されてます。
パワースペクトル密度（Power Spectral Density）
ぶっちゃけ、対象の信号にFFTかけて適切な処理してプロットするのがScilab comm_tbxのplot_psd()関数のようです。上記のページにも書いてありますが、計算にあたっては「ハニング窓」使えということです。plot_psd()もそうしているようです。
comm_tbxのドキュメントを参照するとplot_psd()の大まかな説明は以下みたいっす。知らんけど。
DFTかけて（ハニング窓で）、ノルムをとって二乗して（パワー）、常用対数とって１０倍すればdB（パワーなので１０倍）だ、ということみたい。
前回例の8PSK例題データをPSDプロットしてみる
前回は以下のような処理をして最後plot_eye()つかってアイプロットしてました。最後のplot_eye()をplot_psd()に変更すれば、パワースペクトル密度プロットを得ることができます。
fs = 1e6;
fi = 0;
fsym = 100e3;
mod = mod_init('8psk', fs, fi, fsym);
testpat = prbs(1e6);
[mod, x] = mod_process(mod, testpat);
xn = awgn(x, 0.1);
clf();
plot_psd(xn);
プロット結果が以下に。
上記の信号列xnは複素信号です。その場合、plot_psd()は正規化周波数で-Fs/2～+Fs/2の区間に対してプロットを行ってくれるようです(正規化周波数で0.5Hzがナイキスト周波数。）
上記のプロットに使ったデータはawgn()関数（正規化ガウス分布）でσ=0.1と少し「盛った」ノイズを載せてます。以下のようにノイズをσ=0.01と控えめにして数列を再生成してみました。
xn2 = awgn(x, 0.01); clf(); plot_psd(xn2);
このときのプロットが以下に。上のプロットに比べると「すこし控えめ？」な感じがしないでもない。ホントか？
実数ベクトルへの適用
上記では、I/Qで表現される複素信号としてサンプル数列を作ってました。plot_psd()関数は、複素信号だけでなく実数信号も扱うこともできて、第７５回では既に実数列でお世話になってました。
FSK(設定は第７５回参照）のときの例が以下に。
実数列に対しては、横軸は正規化周波数で0Hz～+Fs/2(正規化周波数0.5Hz)の正の範囲になるようです。
ASK（デフォルト設定）のときの例が以下に。
スペアナみたいなグラフが描けていい感じ？