Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
Help us understand the problem. What is going on with this article?

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

More than 1 year has passed since last update.

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

そうすると,こんな感じのプロットができます.

美しいですね.
ちなみに,愚直にループ毎に二次元配列を生成すると重くて動きが鈍かったです.
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);
}

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

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

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

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

参考

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

DakoMusen
趣味で無線通信と信号処理してます
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away