リアルタイム3Dプロットの目的
- 楽しい
- 可愛い
概要
- テスト用コード
- 応用（周波数スペクトルの時間変化をプロット）
テスト用コード
まずはテスト用のコードを見ながら楽しさを実感しましょう．
事前にgnuplotをインストールしておく必要があります．
私は今回ターミナルにwxtを使用しました．
test.cpp
#include <vector>
#include <cmath> //for sin, cos
int main(){
int loop = 1000;
int N = 256;
int K = 16;
FILE *gp;
gp = popen("gnuplot", "w");
fprintf(gp, "unset key\n");
fprintf(gp, "set xrange[0:%d]\n", N);
fprintf(gp, "set yrange[0:%d]\n", K);
fprintf(gp, "set zrange[-2:2]\n");
fprintf(gp, "set xlabel \"Time\"\n");
fprintf(gp, "set ylabel \"Time\"\n");
fprintf(gp, "set zlabel \"Amplitude\"\n");
std::vector<std::vector<double> > x(K, std::vector<double>(N,0));
int k = 0;
for(int i = 0; i < loop; i++){
for(int n = 0; n < N; n++){
x[k][n] = std::sin(2.0 * std::acos(-1.0) * (n + i * 10) / N);
}
fprintf(gp, "splot \"-\" with points \n");
for(int j = 0; j < K; j++){
for(int n = 0; n < N; n++){
fprintf(gp, "%d, %d, %lf\n", n, j, x[(K + k - j)%K][n]);
}
}
fprintf(gp, "e\n");
fflush(gp);
k = (k + 1)%K;
}
}
そうすると，こんな感じのプロットができます．
リアルタイム3Dプロットのテストコード pic.twitter.com/kheSD3xvZt
— μsen (@DakoMusen) 2018年12月27日
美しいですね．
ちなみに，愚直にループ毎に二次元配列を生成すると重くて動きが鈍かったです．
gnuplot部分↓
FILE *gp;
gp = popen("gnuplot", "w");
fprintf(gp, "unset key\n");
fprintf(gp, "set xrange[0:%d]\n", N);
fprintf(gp, "set yrange[0:%d]\n", K);
fprintf(gp, "set zrange[-2:2]\n");
fprintf(gp, "set xlabel \"Time\"\n");
fprintf(gp, "set ylabel \"Time\"\n");
fprintf(gp, "set zlabel \"Amplitude\"\n");
fprintf(gp, "splot \"-\" with points \n");
for(int j = 0; j < K; j++){
for(int n = 0; n < N; n++){
fprintf(gp, "%d, %d, %lf\n", n, j, x[(K + k - j)%K][n]);
}
}
fprintf(gp, "e\n");
fflush(gp);
ターミナルによってはgp = popen("gnuplot -persist", "w");とすれば，
プログラムが終了しても画面が残ります．
また，メインが終了する前にsystem("pause")などでプログラムを一時停止することもできます．
描画しているデータ↓
for(int n = 0; n < N; n++){
x[k][n] = std::sin(2.0 * std::acos(-1.0) * (n + i * 10) / N);
}
しがないsin波ですが，ループ毎にちょっとズラしておくと綺麗です．
応用(周波数スペクトルの時間変化をプロット)
以下のように書いてみました．
信号生成やFFT，アンプ，雑音はここでは冗長なので省略しました．
vectorは長いのでusingしておくと便利です．
main.cpp
#include <iostream>
#include <vector>
#include <complex>
using COM = std::complex<double>;
using V = std::vector<double>;
using VV = std::vector<V>;
using VCOM = std::vector<COM>;
using VVCOM = std::vector<VCOM>;
int main(){
int N = 256; //FFTサイズ
int K = 16; //
double sps = 50; //MHz
VVCOM wave(K, VCOM(N,0));
VV power_dB(K, V(N));
//値ファイル
int Lfile = 2; //ファイルの長さ
V Dfile(Lfile,0); //0 -> スイッチ，1 -> 信号電力
FILE *fp;
fp = fopen("value.csv", "r");
//アンプ
Amplifier amp;
/*
略
*/
// gnuplot
FILE *gp;
gp = popen("gnuplot", "w");
fprintf(gp, "unset key\n");
fprintf(gp, "set xrange[%lf:%lf]\n", -sps/2.0, sps/2.0);
fprintf(gp, "set yrange[0:%d]\n", K);
fprintf(gp, "set zrange[-150:50]\n");
fprintf(gp, "set zlabel rotate by 90\n");
fprintf(gp, "set xlabel \"Frequency [MHz]\"\n");
fprintf(gp, "set ylabel \"Sample time\"\n");
fprintf(gp, "set zlabel \"Power (dBm/Hz)\"\n");
fprintf(gp, "set isosamples 50\n");
int k = 0;
while(Dfile[0] == 0){
for(int i = 0; i < Lfile; i++){
fscanf(fp, "%lf", &(Dfile[i]) ); //外部入力
}
rewind(fp);
/*
信号生成
IFFT
アンプ入力
雑音入力
FFT
*/
//Power
for(int i = 0; i < N/2; i++){
power_dB[k][i + (N / 2)] = 10.0 * log10(pow(abs(wave[k][i]),2));
power_dB[k][i] = 10.0 * log10(pow(abs(wave[k][i + (N / 2)]),2));
}
fprintf(gp, "splot \"-\" with pm3d\n"); //カラーマップ使用
for(int i = 0; i < N; i++){
for(int j = 0; j < K; j++){
fprintf(gp, "%lf, %d, %lf\n", -sps/2.0 + (sps/(double)N) * (double)i, j ,power_dB[(K + k - j)%K][i]);
}
fprintf(gp, "\n");
}
fprintf(gp, "e\n");
fflush(gp);
fprintf(gp, "set label 5 center at screen 0.5,0.9 \"Power = %d dBm\" \n", (int)Dfile[1]);
k = (k + 1)%K;
}
fclose(fp);
pclose(gp);
}
すると，こんな感じのプロットが見られます．
スペクトル波形に時間軸を追加して見た
— μsen (@DakoMusen) 2018年12月25日
電力上げてくと増幅器で非線形増幅するから高調波がリアルタイムで見られる
うねうねしてて可愛い pic.twitter.com/wUL0mRiOKL
これで高調波をリアルタイムで見られます（2Dの方が見やすかったかも）．
雑音を入れるとうねうねしてて可愛いですね．
fprintf(gp, "splot \"-\" with pm3d\n");としておくと，
カラーマップを使用できるので3Dプロットのときはオススメです．
次は外部装置からのデータを入力して，3Dのオシロスコープやスペクトルアナライザを作ってみたいですね．