convolution

Gaucheが半年ぶりに新しくなった。母語なんで追従しておく。なんか、簡単にインストール出来たり、ドッカーで手軽に試せるようになってたりで、頭が下がりますよ。

かの昔、SICPの読み合わせ会に出てた時、左畳み込みがいいか右畳み込みがいいか喧々諤々な議論になったのを思い出す。畳み込みなんで、foldって実に分かり易いネーミングだな。

[ham@cent tmp]$ gosh -V
Gauche scheme shell, version 0.9.8 [utf-8,pthreads], x86_64-pc-linux-gnu
[ham@cent tmp]$ gosh
gosh> (fold + 0 '(1 2 3 4 5 6 7 8 9 10))
55

引数のリストをプラス演算子を使って畳み込み(0は初期値)。結果は縮小されました。そう解釈すれば、平均なんてのも一種の畳み込みだわな。

畳み込みで検索してたら、懐かしいHPが出て来た。オイラーと同じ香りがする。

アマチュアのためのアマチュアによるアマチュアのTips

octaveに魅せられたのは、シミュレーションで学ぶディジタル信号処理の本かららしい。気が付かなかったわい。勉強される方は、ちゃんとこういう本で根本原理を習得されてるのね。

それに引き換え、今のpythonブーム。部品を糊で繋ぎ合わせただけで、分かった気分にさせてる。底が浅いように思うんだけど、どうよ？

help of octave

octaveのhelpの使い方を調べようとして、help helpしたら、help --listが有るよと。

*** operators:
 :
*** keywords:
 :
*** builtins:
 :
*** functions in /app/share/octave/5.1.0/m/time:

addtodate.m     clock.m         datenum.m       eomday.m        now.m
asctime.m       ctime.m         datestr.m       etime.m         weekday.m
calendar.m      date.m          datevec.m       is_leap_year.m

ってな具合に有益な情報が出てくる。けど、画面が流れてしまって、最初の方を見るには、tmuxのバックスクロールのお世話になる。ちょっと不便。パイプに流し込めたら検索も容易だろうに。何か解決策は無いかな。

こういう時は、helpを実現してるソースに当たれ。

/usr/share/octave/x.x.x/m/help/help.m (/app/share/octave/5.1.0/m/help/help.m)

    if (strcmp (name, "--list"))
      list = do_list_functions ();
      if (nargout == 0)
        printf ("%s", list);
      else
        retval = list;
      endif
      return;
    endif

ふむ、値を受け取る変数を指定すれば、画面に表示する代わりに、お持ち帰り出来るのね。そんな事、説明されていなかったぞ。

octave:13> text = help('--list');
octave:14> size(text)
ans =

       1   33019

octave:15> text
 :
addtodate.m     clock.m         datenum.m       eomday.m        now.m
asctime.m       ctime.m         datestr.m       etime.m         weekday.m
calendar.m      date.m          datevec.m       is_leap_year.m

後は、このtextをファイルに落とせばいいんだな。そうしておけばunix側でいかようにも出来る。どうやってファイルに落とす。writeぐらいで、検索(man -k key 相当)すればいいな。

octave:16> lookfor write
audiowrite          Write audio data from the matrix Y to FILENAME at sampling
                    rate FS with the file format determined by the file extensi
                    ion.
fputs               Write the string STRING to the file with file descriptor FI
                    D.
  :

fputsが良さそうだから、使ってみるか。

octave:17> fd = fopen('ocfunc.txt', 'w');
octave:18> fputs(fd, text);
octave:19> fclose(fd);
octave:21> ls -l ocfunc.txt
-rw-rw-r-- 1 ham ham 33019 Jun 24 16:22 ocfunc.txt

ここれ辺のオペレーションはunixと一緒だな。

octave:24> [~, out] = system('grep help ocfunc.txt');
octave:25> disp(out);
  :
__makeinfo__.m                 help.m

systemの使い方をhelpしてたら、more on/off なんてのがさりげなく使われていた。onにするとページャーが働き出すのね。これで、長い表示も安心。

which, type, who(s)

octave:34> which help
'help' is a function from the file /app/share/octave/5.1.0/m/help/help.m
octave:35> which diff
'diff' is a built-in function from the file libinterp/corefcn/data.cc

所在地の確認。

octave:36> who
Variables in the current scope:

ans   fd    out   text

octave:37> whos
Variables in the current scope:

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  =====
        ans         1x1                          8  double
        fd          1x1                          8  double
        out         1x296                      296  char
        text        1x33019                  33019  char

Total is 33317 elements using 33331 bytes

変数の一覧表。

octave:38> type fd
fd is a variable
 3
octave:39> type diff
diff is a built-in function
octave:40> type help
help is the user-defined function defined from: /app/share/octave/5.1.0/m/help/help.m

## Copyright (C) 2009-2019 Søren Hauberg
##
## This file is part of Octave.
 :

何と、ソースまで見えます。これ中々便利だな。

Octave Forge Documentation

flatpak 上の octave

なんか G線上のアリアを気取った積もりは無いんだけど、上記の構成で不都合な事がある。通常は、/tmpに移動して、そこでファイルを作ったり消したりしてる。が、それだとoctave側からファイルが見えないんだ。homeの中だと問題無いんだけどね。

これ、flatpakの弊害ではなかろうか？ systemとかunixって関数を使うと、何不自由なく unixのオペレーションが出来るので、少し調べてみた。

octave:1> ls /tmp
octave:2> ls -l /
合計 4
drwxr-xr-x   8 ham ham  115 Jun 15 06:41 app
lrwxrwxrwx   1 ham ham    7 Jun 25 14:04 bin -> usr/bin
drwxr-xr-x   5 ham ham  360 Jun 25 14:04 dev
drwxr-xr-x  15 ham ham  840 Jun 25 14:04 etc
drwxr-xr-x   4   0   0   30 Nov 20  2018 home
lrwxrwxrwx   1 ham ham    7 Jun 25 14:04 lib -> usr/lib
lrwxrwxrwx   1 ham ham    9 Jun 25 14:04 lib64 -> usr/lib64
drwx------   2   0   0    6 May  3  2018 lost+found
drwxr-xr-x   2   0   0    6 Apr 11  2018 media
drwxr-xr-x   2   0   0    6 Apr 11  2018 mnt
drwxr-xr-x   4   0   0   30 Apr 20 15:37 opt
dr-xr-xr-x 201   0   0    0 Jun 25 14:04 proc
drwxr-xr-x   6 ham ham  160 Jun 25 14:04 run
lrwxrwxrwx   1 ham ham    8 Jun 25 14:04 sbin -> usr/sbin
drwxr-xr-x   2   0   0    6 Apr 11  2018 srv
drwxr-xr-x   7 ham ham  140 Jun 25 14:04 sys
drwxr-xr-x   3 ham ham   60 Jun 25 14:04 tmp
drwxr-xr-x  16 ham ham 4096 Jun 15 06:41 usr
drwxr-xr-x   6 ham ham  140 Jun 25 14:04 var

普通、/tmpの下には忌まわしいsystemdの運用用のファイルが有るはずだけど、見えない。更にrootから下を見下ろすと、卒倒するようなパーミション関係になってる。

octave:4> ls /app
bin  include  jre  lib  libexec  manifest.json  share
octave:5> ls /app/bin
epstool              install-info     octave-5.1.0         remote_lat
fftw-wisdom          libpng-config    octave-cli           remote_thr
fftw-wisdom-to-conf  libpng16-config  octave-cli-5.1.0     texi2any
fftwf-wisdom         local_lat        octave-config        texi2dvi
gm                   local_thr        octave-config-5.1.0  texi2pdf
gnuplot              makeinfo         pdftexi2dvi          texindex
info                 mkoctfile        png-fix-itxt
inproc_lat           mkoctfile-5.1.0  pngfix
inproc_thr           octave           pod2texi

こういうの、flatpakが用意した砂場だよね。安全な代わりに制限を受けるとな。 Sandbox Permissionsこれだな。そんな訳なんで、ちょびっと、起動スクリプトを修正した。

[ham@cent ~]$ cat /usr/local/bin/octave
#! /bin/sh
flatpak run --share=network --filesystem=/tmp \
  --command=sh org.octave.Octave -c "DISPLAY=$DISPLAY octave $@"

物はついでなんで、大事なPATHも確認しておくかな。

octave:4> path = getenv('PATH');
octave:5> path
path = /app/bin:/usr/bin:/app/jre/bin:/app/libexec/octave/5.1.0/site/exec/x86_64-pc-linux-gnu:/app/libexec/octave/api-v53/site/exec/x86_64-pc-linux-gnu:/app/libexec/octave/site/exec/x86_64-pc-linux-gnu:/app/libexec/octave/5.1.0/exec/x86_64-pc-linux-gnu:/app/bin
octave:6> lookfor split
split_long_rows     Query or set the internal variable that controls whether ro
                    ws of a matrix may be split when displayed to a terminal wi
                    indow.
ostrsplit           Split the string S using one or more separators SEP and ret
                    urn a cell array of strings.
strsplit            Split the string STR using the delimiters specified by DEL
                    and return a cell string array of substrings.

ちょいと醜いので分割したいな。分割コマンドってあるのかな？

octave:7> strsplit
error: Invalid call to strsplit.  Correct usage is:

 -- [CSTR] = strsplit (STR)
 -- [CSTR] = strsplit (STR, DEL)
 -- [CSTR] = strsplit (..., NAME, VALUE)
 -- [CSTR, MATCHES] = strsplit (...)

octave:7> strsplit(path, ':')
ans =
{
  [1,1] = /app/bin
  [1,2] = /usr/bin
  [1,3] = /app/jre/bin
  [1,4] = /app/libexec/octave/5.1.0/site/exec/x86_64-pc-linux-gnu
  [1,5] = /app/libexec/octave/api-v53/site/exec/x86_64-pc-linux-gnu
  [1,6] = /app/libexec/octave/site/exec/x86_64-pc-linux-gnu
  [1,7] = /app/libexec/octave/5.1.0/exec/x86_64-pc-linux-gnu
  [1,8] = /app/bin
}

helpを引くのが面倒だと、ちょいと怒らせて、情報を仕入れます。そして実行しました。

octave:9> unix('java -version')
openjdk version "11.0.2" 2019-01-15
OpenJDK Runtime Environment (build 11.0.2+9)
OpenJDK 64-Bit Server VM (build 11.0.2+9, mixed mode)

こういう物が入っているから、膨れ上がるんだな。

compile octave

お決まりのソースからコンパイルするコースをやってみる。コンパイルには色々なツールが必要になるんで、一発で用意する方法が用意されてた。

apt-get build-dep octave
yum-builddep octave

Debianコースって事で、i386なマシンに入れる事にした。configureが警告を出してきた。あれが見つからなかったんで、これは出来ないとかね。取り合えずそんなのは無視して、gdbにかけられるようになれば、取り合えず満足。

コンパイルの時間を有効利用ってんで、ソースツリーをうろうろしてたら、リファレンスカードを発見。octaveするなら、これぐらいの事は知っておいて欲しいぞってのが列挙されてた。

そして、ネットには、信号処理のための関数ってな事で、特定分野のやつもあったぞ。

real    107m36.304s
user    100m53.824s
sys     3m18.948s
debian:octave-4.4.1$

待つ事2時間近く、3GのDiskを消費して(gdb用のシンボルテーブル付きの場合)、無事にコンパイルとインストールを完了した。

対話型のアプリなんで、立ち上げておいてから、gdbでアタッチするのが良い。下記は、convnを追った図。

gdb -p 1343 ./octave-cli
 :
Thread 1 "octave-cli-4.4." hit Breakpoint 1, convn (a=..., b=..., ct=convn_full) at liboctave/numeric/oct-convn.cc:198
198     CONV_DEFS ( , )
(gdb) bt
#0  convn (a=..., b=..., ct=convn_full) at liboctave/numeric/oct-convn.cc:198
#1  0xb6fc98ff in Fconv2 (args=...) at libinterp/corefcn/conv2.cc:211
#2  0xb6cf7ca2 in octave_builtin::call (this=0x7650e0, tw=..., nargout=<optimized out>, args=...) at libinterp/octave-value/ov-builtin.cc:65
#3  0xb6ee1afd in octave::tree_evaluator::visit_index_expression (this=0x741d10, idx_expr=...) at libinterp/parse-tree/pt-eval.cc:1381
#4  0xb6d4e1aa in octave::tree_evaluator::evaluate (this=0x741d10, expr=0x849c90, nargout=<optimized out>) at libinterp/parse-tree/pt-eval.h:300
#5  0xb6edd78e in octave::tree_evaluator::visit_statement (this=<optimized out>, stmt=...) at libinterp/parse-tree/pt-eval.cc:2337
#6  0xb6ed27c8 in octave::tree_statement::accept (tw=..., this=<optimized out>) at libinterp/parse-tree/pt-stmt.h:115
#7  octave::tree_evaluator::visit_statement_list (this=0x741d10, lst=...) at libinterp/parse-tree/pt-eval.cc:2406
#8  0xb72725ca in octave::tree_statement_list::accept (tw=..., this=<optimized out>) at libinterp/parse-tree/pt-stmt.h:190
#9  octave::interpreter::main_loop (this=<optimized out>) at libinterp/corefcn/interpreter.cc:995
#10 0xb7276271 in octave::interpreter::execute (this=<optimized out>) at libinterp/corefcn/interpreter.cc:717
#11 0xb68a0652 in octave::cli_application::execute (this=0xbff9312c) at libinterp/octave.cc:391
#12 0x00489e82 in main (argc=<optimized out>, argv=<optimized out>) at src/main-cli.cc:93

その後を追うと

(gdb) s
[Switching to thread 2 (Thread 0xaf6b7b40 (LWP 7476))](running)
0xb5800e95 in __x86.get_pc_thunk.di () at liboctave/array/Array.cc:550
550       ~rec_index_helper (void) { delete [] idx; delete [] dim; }
(gdb)
[Switching to thread 2 (Thread 0xaf6b7b40 (LWP 7476))](running)
convolve<double, double> (a=..., b=..., ct=convn_full) at liboctave/numeric/oct-convn.cc:132
132     convolve (const MArray<T>& a, const MArray<R>& b,

Cフラフラ語は、よう追えんわ。と、敗北宣言。

convolution

コンボリューション関数（conv）と相関関数(xcorr)と離散フーリエ変換関数（fft）の関係

conv by matlab

Matlab(Octave) による信号処理のシミュレーション基礎

select matrix (or vector)

   This shows an important general concept about using arrays for
indexing instead of looping over an index variable.  *Note Index
Expressions::.  Also use boolean indexing generously.  If a condition
needs to be tested, this condition can also be written as a boolean
index.  For instance, instead of

     for i = 1:n
       if (a(i) > 5)
         a(i) -= 20
       endif
     endfor

write

     a(a>5) -= 20;

which exploits the fact that ‘a > 5’ produces a boolean index.

これ、前回、目を付けておいた高速テクニック。forループを止めて、簡潔に書けるとな。簡潔==高速ってのは、コンピュータ界の真理だからな。

例では、要求を満たすエレメントを書き換えしてるけど、選び出すにはどうやれば良い？

octave:2> a = [3 6 1 2 9 4 8 2];
octave:3> out = [];
octave:4> for i = 1:columns (a)
> if ( a(i) > 5)
> out = [out a(i)]
> endif
> endfor
out =  6
out =

   6   9

out =

   6   9   8

Lispと言うかSchemeの知識を総動員して、アペンドもどきをやってみた。でも、もっと簡単になるはず。知恵を出せ。

octave:15> a
a =

   3   2
   5   1
   6   6
   3   3
   1   2
   7   6
   4   3

こんな行列の1列目が5未満のものだけを、抽出したい。

octave:16> a(a(:,1)<5,:)
ans =

   3   2
   3   3
   1   2
   4   3

1行で、指定できた。素晴しい表現力! これって、pythonやHaskellで自慢の種になる pythonの内包表記を少し詳しくもどきじゃなかろうか。

早速、応用してみよう。時々登場する、血圧記録データは、測定した年月日時でタグ付けされてる。タグを頼りに、起床時のデータだけを取り出してみる。

octave:1> bld = csvread('bld.csv');
octave:2> bld(1:4,:)
ans =

   12010104        118         73         57
   12010121        109         70         71
   12010204        128         82         60
   12010221        105         62         67

octave:3> bld( rem(bld(:,1),100)<12, : )
ans =

   12010104        118         73         57
   12010204        128         82         60
   12010304        122         81         59
   12010405        120         80         60
   12010505        136         84         55

remは余りを求める関数。この引数はbldの１列目だけをを100で割ってる。それが12時未満のものを対象に、全データを抽出。

csvreadって関数も、相手にするのは数値のみってのが潔くてなかなな良い選択だと思うよ。何たって、数値演算と言うか行列の扱いが得意なoctaveですからね。

octave:6> size(bld)
ans =

   10    4

こんな具合に、普通に行列になってるよ。

octave:12> bld(end,:)
ans =

   12010521        111         67         64

octave:13> bld(end-2:end,:)
ans =

   12010421        132         82         61
   12010505        136         84         55
   12010521        111         67         64

tailからのアクセスも慣れると簡単。