power_spectrum の使い方

power_spectrum は ADS3KP + OCTAVIA のデータのパワースペクトルを出すソフトです。バンキャラ出す程度なら計算は一瞬で終わるのですが、MATLAB ランタイムの呼び出しに 10 秒近くかかります。

ファイルリスト

アプリ本体 (CentOS7, Mac)
インターフェーススクリプト
テストデータ (580 MHz 正弦波を 2 bit, 1024 MS/s で AD 変換したデータ)

インストール

以降 Linux ユーザー名は適宜読み替えてください。

適当なディレクトリにソフトをダウンロードし、実行権限を与え、実行してください。

$ wget http://astro.sci.yamaguchi-u.ac.jp/t-aoki/tools/pulse_searcher/power_spectrum_cent7gnome.install
$ chmod 555 power_spectrum_cent7gnome.install
$ ./power_spectrum_cent7gnome.install

ウィンドウ表示に少し時間がかかりますが、指示に沿ってインストールしてください。このアプリのインストールと MATLAB Runtime のインストールの二段階があります。インストール先は pulse_searcher と同じにして良いと思いますが、以下ではこちらの環境で記述します。

アプリのインストール先: /home/oper/MATLAB/power_spectrum

MATLAB Runtime のインストール先: /home/oper/MATLAB/MATLAB_Runtime

これですべてのファイルは /home/oper/MATLAB 以下に入ったはずです。

実行

バンキャラ出しに使うパラメータは限られるので、インターフェースとなるスクリプトを使うと便利です。使い方はデータファイルを引数に与えるだけです。

$ ./spectrum.sh hoge.vdif

このスクリプトのエイリアスをシステムに登録しておくとバンキャラ出しが楽になると思います。

$ vi ~/.bashrc
alias spectrum="/home/oper/MATLAB/spectrum.sh"

$ . ~/.bashrc

こうしておけばバンキャラ出しは以下のコマンドだけで出せます。

$ spectrum u18001.vdif

上記スクリプトの内容を見ればわかる通り、アプリ本体は power_spectrum で、次のように実行します。

$ ./power_spectrum {file} {nfft} {integration_s} {stft_s} {format} {sideband} {gpu_on} {figure_on} {outdir}

各引数は次の通りです。

{file} = データファイル
{nfft} = 1024 ... FFT 窓幅 (sample 単位)
{integration_s} = 0.1 ... 積分時間 (second 単位) の目安。積分時間が実際に指定した時間になるかはわかりません。標準出力の結果を確認してください。
{stft_s} = 0.1 ... STFT 幅 (second単位) で、一度にメモリ (main/GPU memory) に読み込むデータ幅にも相当し、処理速度に影響します。積分時間より短い必要があります。
{format} = ADS1K と指定してください。意味は ADS3KP + OCTAVIA のビット並びで、実際のADS1Kには対応していません。
{sideband} = USB ... 実数 FFT 後のスペクトルデータ上で有効な側波帯を指定。他にLSB, DSBを指定可能。
{gpu_on} = ON or OFF ... バンキャラ出すだけなら OFF (GPU使わない) で十分です。
{figure_on} = ON or OFF ... 解析のたびにグラフウィンドウを出力します。複数ファイルを連続処理するときは OFF に設定しないと、処理がいちいち止まります。
{outdir} = ./ ... テキストファイルを出力するディレクトリ。書き込み権限があることを確認してください。

出力

解析終了後テキストファイルが {outdir} に出力され、さらに引数 {figure_on} が ON の場合ウィンドウが開きスペクトルが表示されます。例えば入力ファイルが YAMAGU32_2018130143000.vdif (ヘッダなし) の場合に、出力スペクトルファイルは YAMAGU32_2018130143000_powerspectrum_fft1024_avg100000_rbw1000.000kHz.txt というファイル名になります。文字列の意味は下記の通りです。

fft1024 ... FFT 窓幅が 1024 sample (512 点分光) のスペクトルであることを示します。
avg100000 ... 10 万個のスペクトルを積分したことを示し、スペアナでいう AVG 100000 sweep です。例えば積分時間は (avg 1E+5 window) * (fft 1024 sample/window) / (sampling freq 1024E+6 sample/s) = 0.1 s になります。
rbw1000.000kHz ... Resolution bandwidth が 1000 kHz であることを表します。

ファイルの中身は以下のようなタブ区切りテキストで、spectral power は線形値で中心周波数成分が 1 となるように規格化しています。

#Freq/MHz	LinearPower/a.u.
1.000000	4.935
2.000000	4.913
3.000000	4.931
4.000000	5.007
5.000000	5.059

引数 {figure_on} が ON の場合、パワースペクトルが線形スケールと対数スケールの両方で表示されます。

出力例: 周波数 f の正弦波信号を、標本化周波数 fs (> 2f) で普通にサンプリングした場合

10 MHz の正弦波信号を、標本化周波数 1024 MHz で標本化したデータに対し、8192 点 FFT 処理した結果のスペクトルを以下に示します。正しく 10 MHz の信号が検出されています。量子化ノイズ、入力信号、FFTノイズフロアの関係についてはアナログ・デバイセズの資料 MT-001 などを参照してください。

出力例: 周波数 f の正弦波信号を、標本化周波数 fs (< 2f) で高次サンプリングした場合

580 MHz の正弦波信号を、標本化周波数 1024 MHz で標本化したデータ (当該データ) に対し、8192 点 FFT した結果のスペクトルを以下に示します。標本化周波数 1024 MHz の高次サンプリングによって、512 MHz 以上の周波数成分は 512 MHz を軸に折り返されて AD 変換されるが、以下の画像は反転した帯域を補正したものになっています。正しく 580 MHz の信号が検出されています。

出力例: メーザー天体

G9.62+0.20E のメタノールメーザー (6668.33 MHz) について、本アプリで得たスペクトルとスペアナで得たスペクトルを以下に示します。両者の結果は一致していることがわかります。