C++とgnuplotでリアルタイム3Dプロット

#リアルタイム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);
}

すると，こんな感じのプロットが見られます．

スペクトル波形に時間軸を追加して見た
電力上げてくと増幅器で非線形増幅するから高調波がリアルタイムで見られる
うねうねしてて可愛い pic.twitter.com/wUL0mRiOKL

— μsen (@DakoMusen) 2018年12月25日

これで高調波をリアルタイムで見られます（2Dの方が見やすかったかも）．
雑音を入れるとうねうねしてて可愛いですね．

fprintf(gp, "splot \"-\" with pm3d\n");としておくと，
カラーマップを使用できるので3Dプロットのときはオススメです．

次は外部装置からのデータを入力して，3Dのオシロスコープやスペクトルアナライザを作ってみたいですね．

#参考

C言語でGnuplotを動かしてアニメーションをする