こんにちは世界。

内容
VcXsrvを入れる
gnuplot を使えるようにする
Raku の Chart::Gnuplot を使う
Raku で計算してプロット
png を動画にする
アニメ化
余弦波の和がデルタ関数になる理由
References

内容

WSL の Ubuntu20.04 で gnuplot を使おうとしたときのメモです。VcXsrv を使います。あと Raku の Chart::Gnuplot を使ってみました。

VcXsrvを入れる

https://sourceforge.net/projects/vcxsrv/ からVcXsrvをダウンロードし、

を参考にVcXsrvを入れます。

ただし私の環境では .profile の末尾には export DISPLAY=:0 を追加しました（nameserver からとってくるIPアドレスが windows のIPアドレスと異なるIPアドレスなのでうまくいかなかった）。ここの設定は環境によって変わりそうです。

https://qiita.com/nishemon/items/bb3aca972404f68bfcd6 に書いてあるようにwindowsのフォントを使えるようにシンボリックリンクを作る

sudo apt install fontconfig
sudo ln -s /mnt/c/Windows/Fonts /usr/share/fonts/windows
sudo fc-cache -fv
fc-list

windowsのフォントがいっぱい表示されればOKです。

gnuplot を使えるようにする

インストールしたものをまとめておくと

sudo apt install libsdl2-dev libcairo2-dev gnuplot gnuplot-x11 eog yaru-theme-icon yaru-theme-gtk

gnuplot したときに Error loading theme icon 'document-save-as' for stock: Icon 'document-save-as' not present in theme Yaru が出たのでテーマを入れています。これを入れないとファイルを保存のアイコン画像が表示されなかったりします。

Raku の Chart::Gnuplot を使う

Raku の Chart::Gnuplot を使うには

zef install SDL2::Raw
zef install Cairo
zef install Chart::Gnuplot

particle や snake などのデモができます。

忘れていたので aliasにraku='rakudo'を登録しました。 ~/.bashrc の末尾に alias raku='rakudo' を追加し、~/.profile に source ~/.bashrc を追加して設定が反映されるようにします。

/home/ユーザ名/.p6chart-gnuplot/bin にある gnuplot が自分で入れたやつじゃない（verは一緒だけどpatchlevelが違った）のでシンボリックリンクを作りました。

which gnuplot

で gnuplot の場所を探したら /usr/bin/gnuplot に置かれていたので

cd ~/.p6chart-gnuplot/bin
ln -s /usr/bin/gnuplot gnuplot

Raku で計算してプロット

今回新たに

zef install Math::Libgsl::Function
zef install Math::Libgsl::Constant

をインストールしました。

前回インストールした Math::Libgsl::Vector のバージョンが古いままでした。 Vector::View の subvector が使いたいので

zef update
zef upgrade

を実行しました。

Chart::Gnuplot は動いたんですけど multiplot があまりうまく動きませんでした。番号の割り振りに問題があるのかな？

use Math::Libgsl::Vector;
use Math::Libgsl::Function :trig;
use Chart::Gnuplot;

my constant f0 = 1.0;
my constant fmax = 100.0;
my constant fs = fmax * 2; # to plot smooth line
my constant xmax = 2.5 / f0;
my constant Ns = ceiling(2 * fs * xmax); # number of sample points

my Math::Libgsl::Vector $t .= new(:size(Ns));
my Math::Libgsl::Vector $v1 .= new(:size(Ns));
my Math::Libgsl::Vector $v2 .= new(:size(Ns));
my Math::Libgsl::Vector $average .= new(:size(Ns));


$t[$_] = $_ / fs for 0..^Ns;
$t.add-constant(-1 * ceiling(Ns / 2) / fs);
$v1[$_] = gsl-cos(2 * pi * f0 * $t[$_]) for 0..^Ns;
$v2[$_] = gsl-cos(2 * pi * 2 * f0 * $t[$_]) for 0..^Ns;
$average.setall(0); # initialize

#my @x-data = 0 xx Ns;
#my @y-data = 0 xx Ns;
#my @x-data = 0..20;
#{@x-data[$_] = $t[$_]} for 0..^Ns;
#{@y-data[$_] = $v1[$_]} for 0..^Ns;
my @vertices1[Ns];
my @vertices2[Ns];
#{@vertices[$_] = [$t[$_], $v1[$_]]} for 0..^Ns;
{@vertices1[$_] = [$t[$_], $v1[$_]]} for 0..^Ns;
{@vertices2[$_] = [$t[$_], $v2[$_]]} for 0..^Ns;

my $gnu = Chart::Gnuplot.new(:terminal("png"), :filename("cos.png"));
$gnu.command("set multiplot layout 2,1");
#$gnu.command("plot sin(x)");
#$gnu.command("plot cos(x/2)");
# $gnu.title(:text("test"));
$gnu.plot(:vertices(@vertices1), :style("lines"), :title("cos(2pi ft)"));
sleep 3.0;
#$gnu.plot(:function("sin(x)"), :style("lines"), :title("cos(2pi ft)"));
$gnu.command("plot cos(x/2)");
$gnu.command("unset multiplot");
#$gnu.command("set terminal x11");
#$gnu.command("set output");
#$gnu.command("exit");
$gnu.dispose;
# $gnu.multiplot(:layout([2,1]));

#my $size = 3;

#$v1.setall(1); # $v1 = (1 1 1)
#$v2.setall(2); # $v2 = (2 2 2)
#$v1.add($v2); # $v1 += $v2

#print "$v1[$_] " for 0..^3;
#say "";

なのでいったん dat として出力して gnuplot しました。何一つ整理していないコードですがメモですのでご容赦ください。Raku 内で gnuplot 用のスクリプトを出力して実行させてます。

#!/usr/bin/rakudo

use Math::Libgsl::Vector;
use Math::Libgsl::Function :trig; # to use gsl-cos()
use Chart::Gnuplot;

my constant f0 = 1;
my constant kmax = 50;
my constant fmax =kmax * f0; 
my constant fs = 301; #fmax * 2; # fmax + 2.01; #101; #fmax * 10; #2.01; # actually Nyquist rate is not so important in this program
my constant tmax = 2.5; # 3.0; #2.5 / f0;
my constant Ns = ceiling(2 * fs * tmax); # number of sample points

my Math::Libgsl::Vector $t-vec .= new(:size(Ns));
#my Math::Libgsl::Vector $v1 .= new(:size(Ns));
my Math::Libgsl::Vector $v1 .= new(:size(kmax + 1));
$v1.zero();

sub sum-vector(Math::Libgsl::Vector $v) {
    my $parity = $v.size % 2;
    # my $length-ceil = ceiling($v.size / 2);
    my $half-length = floor($v.size / 2);
    if ($half-length == 0) {
        return $v[0];
    } else {
        my Math::Libgsl::Vector::View $vl .= new;
        my Math::Libgsl::Vector::View $vr .= new;
        my Math::Libgsl::Vector::View $v-accum .= new;
        my Math::Libgsl::Vector $v-half = $vl.subvector($v, 0, $half-length);
        $v-half.add($vr.subvector($v, $half-length + $parity, $half-length).reverse());
        sum-vector($v-accum.subvector($v, 0, $half-length + $parity));
    }
}

my $cos-file-name = 'cos.dat';
my $sum-file-name = 'sum.dat';

sub pnggen($f, $f0) {
    shell "rm $cos-file-name" if "$cos-file-name".IO.e;
    shell "rm $sum-file-name" if "$sum-file-name".IO.e;


    $t-vec[$_] = $_  / fs for 0..^Ns;
    $t-vec.add-constant(-1 * ceiling(Ns / 2) / fs);
    for (0..^Ns) -> $i {
        my $t = $t-vec[$i];
        $v1.zero();
        spurt($cos-file-name, "$t ", :append);
        for (1..$f) -> $k { #for (1..kmax) -> $k {
            $v1[$k] = $f0 * gsl-cos(2 * pi * $k * $f0 * $t);
            # $v1[$k] = gsl-cos(2 * pi * $k * $f0 * $t);
            spurt($cos-file-name, "$v1[$k] ", :append);
        }
        spurt($cos-file-name, "\n", :append);
        my $sum = sum-vector($v1);
        # say "$t $sum";
        spurt($sum-file-name, "$t $sum\n", :append);
    }


    my $file-gp = 'script-gp.txt';
    #my $serial = "$f".fmt('%03d');
    my $serial = "{10/$f0}".fmt('%03d');
    my $plt-str = 'plot "cos.dat" using 1:2 title \'\' w l, ';
    for (1..^$f) -> $r {
        my $r-num = $r + 2;
        $plt-str = $plt-str ~ " '' using 1:$r-num title '' w l, ";
    }
    $plt-str = $plt-str ~ "\n";
    shell "rm $file-gp" if "$file-gp".IO.e;
    spurt($file-gp, "set terminal pngcairo\n", :append);
    spurt($file-gp, "set output 'test-$serial.png'\n", :append);
    ## spurt($file-gp, "set terminal gif\n", :append);
    ## spurt($file-gp, "set output 'test-$serial.gif'\n", :append);
    ## spurt($file-gp, "set size 2.0, 2.0\n", :append);
    spurt($file-gp, "set xrange [-2.5:2.5]\n", :append);
    spurt($file-gp, "set grid\n", :append);
    spurt($file-gp, "set multiplot layout 2,1\n", :append);
    spurt($file-gp, 'set title "f_0 cos(2{/Symbol p}n{f_0}t)' ~ " for n=1:$f, f0 = $f0" ~ '"' ~ "\n", :append);
    spurt($file-gp, "set yrange [-1.1:1.1]\n", :append);
    ## spurt($file-gp, 'plot "cos.dat" using 1:2 w l' ~ "\n", :append);
    spurt($file-gp, $plt-str, :append);
    spurt($file-gp, 'set title "{/Symbol S}^{' ~ "$f" ~ '}_{n = 1} f_0 cos(2{/Symbol p}n{f_0}t)' ~ "\n", :append);
    # spurt($file-gp, 'set title "{/Symbol S}^{/Symbol \245}_{n = 1} cos(2{/Symbol p}n{f_0}t)' ~ "\n", :append);
    # spurt($file-gp, "set yrange [-2:{$f+1}]\n", :append);
    spurt($file-gp, "set yrange [-2:11]\n", :append);
    # spurt($file-gp, "set yrange [-5:30]\n", :append);
    # spurt($file-gp, "set yrange [-5:15]\n", :append);
    spurt($file-gp, 'plot "sum.dat" using 1:2 title "sum" w l' ~ "\n", :append);
    spurt($file-gp, "unset multiplot\n", :append);

    shell "gnuplot-x11 $file-gp";
    # shell "eog test.png";
}
#pnggen(50);

my $f0 = Int(100/30)/10.0;
my $x = ceiling(10/$f0);
for (1..20) -> $ff {
    say "$x, $f0";
    pnggen($x, $f0 + $ff);
}

vector の要素の総和を計算する関数が見つけられなかったのでベクトルを折り返しながら足していく関数を作りました。ただ上手くやらないと積み残しの誤差とかが出ちゃうので気を付けたいです。

書き出されて実行される gnuplot のスクリプトはこんな感じです。

set terminal pngcairo
set output 'test-000.png'
set xrange [-2.5:2.5]
set grid
set multiplot layout 2,1
set title "f_0 cos(2{/Symbol p}n{f_0}t) for n=1:34, f0 = 20.3"
set yrange [-1.1:1.1]
plot "cos.dat" using 1:2 title '' w l,  '' using 1:3 title '' w l,  '' using 1:4 title '' w l,  '' using 1:5 title '' w l,  '' using 1:6 title '' w l,  '' using 1:7 title '' w l,  '' using 1:8 title '' w l,  '' using 1:9 title '' w l,  '' using 1:10 title '' w l,  '' using 1:11 title '' w l,  '' using 1:12 title '' w l,  '' using 1:13 title '' w l,  '' using 1:14 title '' w l,  '' using 1:15 title '' w l,  '' using 1:16 title '' w l,  '' using 1:17 title '' w l,  '' using 1:18 title '' w l,  '' using 1:19 title '' w l,  '' using 1:20 title '' w l,  '' using 1:21 title '' w l,  '' using 1:22 title '' w l,  '' using 1:23 title '' w l,  '' using 1:24 title '' w l,  '' using 1:25 title '' w l,  '' using 1:26 title '' w l,  '' using 1:27 title '' w l,  '' using 1:28 title '' w l,  '' using 1:29 title '' w l,  '' using 1:30 title '' w l,  '' using 1:31 title '' w l,  '' using 1:32 title '' w l,  '' using 1:33 title '' w l,  '' using 1:34 title '' w l,  '' using 1:35 title '' w l, 
set title "{/Symbol S}^{34}_{n = 1} f_0 cos(2{/Symbol p}n{f_0}t)
set yrange [-2:11]
plot "sum.dat" using 1:2 title "sum" w l
unset multiplot

png を動画にする

convert -delay 50 -loop 0 test-*.png matome.gif

を使ってgifアニメにして

eog matome.gif

で確認できます。あとは sudo apt install ffmpeg して

ffmpeg -r 2 -i 'out(%d).png' -pix_fmt yuv420p output.mp4

で mp4 にできます。元画像のサイズは

mogrify -resize 640x480 test-*.png

でリサイズできます。

png ファイルのナンバリングは rename 使った shellスクリプトを書くか、Windows だとファイルをドラッグアンドドロップで選択して番号を付けると一括でナンバリングできるのでそれを使うのが楽です。

アニメ化

こんな感じです。

余弦波の和がデルタ関数になる理由

動画をご覧ください。

https://youtu.be/hQVbLKiYeVo

どうもありがとうございました。

References

文中記載に加えて

Chidori Adachi, dashing on each thing

足立千鳥があちこちぶちあたりながら学んでいく記録です。

WSL で gnuplot を使う

内容