DSP of julia
いささか旧聞になるけど、
5 Programming Languages That Are Probably Doomed
オワンコに挙げられている言語らしい。
PerlやObject-Cが完全にオワンコってのは激しく同意。Haskellは頭が固くて扱いづらい。Rは特定部門に特化してて伸びしろ無し。昔LLでPHPで色々な事が出来るか試されたけど、哀れなくらい駄目だった。それと同じ構図だな。
じゃRubyはどうよ？ みんなレールに乘って楽しようとしたけど、よそ様も同じ事をやるにつれて、優位性が失われた。やはりPython一択か。
まあ、世の中ML（機械学習）とWebが出来ればそれで良しと言う風潮だものな。
昔、コンピューター、ソフト無ければただの箱って標語が有ったけど、それを今風に置き換えれば、
スマートフォン通信出来なきゃ、ただの板
かな。通信はWebでしょ。その裏でMLが出来ればそれで満足。パイソン万歳の風潮に逆らって、今日もjulia道に邁進するのさ。
ギリシャ文字
REPL内から、ギリシャ文字に変換出来るとな。
入力文字列 タブキーで変換できる記号 \Alpha Α \alpha α \beta β \gamma γ \delta δ \epsilon ϵ \zeta ζ \eta η \theta θ \iota ι \kappa κ \lambda λ \mu μ \nu ν \pi π \rho ρ \sigma σ \tau τ \upsilon υ \phi ϕ \chi χ \psi ψ \omega ω \sqrt √
平方根の記号にも変換出来て、数学屋さんにはたまらんでしょう。よく使う大文字のシグマ Σ は、残念ながら変換定義されていないっぽい。
でも、電気で使うオメガが有るから許したる。 sin( ω * t ) と書けるからかっこいい。本当は、sin( ωt ) がいいんだけど、こういう風に変数を続けては書けない。
julia> ω = 2pi 6.283185307179586 julia> t=3 3 julia> ωt ERROR: UndefVarError: ωt not defined
どうせなら、ギリシャ文字は特別扱いしてくれればいいのに。数学者は泣いて喜ぶぞ。
上で例に使っちゃったけど、2 * x + 1 みたいな、数字と変数を乗算する場合、2x + 1 と、数学風に書いて良い事になっている。
それから、大きな数値、10000 みたいのは、10_000 と、書けるように、年寄りモードが用意されてる。若い人よ、年寄りにも優しくしておくれ。そのうちに若者も年寄りになるんだから！
femtoLisp
debian:AA$ julia --lisp ; _ ; |_ _ _ |_ _ | . _ _ ; | (-||||_(_)|__|_)|_) ;-------------------|---------------------------------------------------------- > (cons 1 '(a b c)) (1 a b c) > (exit)
femtoLispが裏方になって、サービスされてるのだろうね。
julia> aa(x) = 2x + 1234 aa (generic function with 1 method) julia> code_warntype(aa, (Int,)) Body::Int64 1 ─ %1 = (Base.mul_int)(2, x)::Int64 │ %2 = (Base.add_int)(%1, 1234)::Int64 └── return %2
llvm風だな。ちょいと面倒。
julia> @code_warntype aa(3.14) Body::Float64 1 ─ %1 = (Base.mul_float)(2.0, x)::Float64 │ %2 = (Base.add_float)(%1, 1234.0)::Float64 └── return %2 julia> @code_warntype aa(π) Body::Float64 1 ─ return 1240.2831853071796
そんな時はマクロの出番で、勝手に解析してくれる。マクロはユーザーフレンドリー。ああ、畳み込みが行われている。
julia> @code_lowered aa(3) CodeInfo( 1 ─ %1 = 2 * x │ %2 = %1 + 1234 └── return %2 )
これは、ちょい見によさそう。
julia> Meta.parse("bb(x) = 4x")
:(bb(x) = begin
#= none:1 =#
4x
end)
julia> Meta.parse("x,y = 12, 34")
:((x, y) = (12, 34))
これが、femtoLispに一番近い出力だろうか。関数型の関数定義って、begin - end で包むと 複数行でも良いとな。こりゃもう、もろに scheme の構文だな。
DSP of julia
octaveからjuliaにシフトして、ふと信号処理ってサポートしてるのか気になった。夜も眠れない悩みだ（単に暑くて眠れないのかもしれないけど）。調べてみた所、
Speech Signal Processing in Julia
こういうのが出てきた。DSPってパッケージが鉄板みたいだ。早速オイラーもやってみる。
音源は、今年の5月末頃にsoxで試した、からおけCDを再び使う。誰が謳ってたか調べたら、 千賀かほるさんの、真夜中のギターって曲だった。
[ham@cent SND]$ soxi S.wav Input File : 'S.wav' Channels : 1 Sample Rate : 8000 Precision : 16-bit Duration : 00:00:10.00 = 80000 samples ~ 750 CDDA sectors File Size : 160k Bit Rate : 128k Sample Encoding: 16-bit Signed Integer PCM
例に倣ってコード（ギターのやつじゃなくてスマソ）を書いたら、REPLみたいになったぞ。 WAVモジュールを使って音源を取り込み、DSPで評価、結果はGnuplotで目に見える形にする。
using WAV
using DSP
using Gnuplot
x, fs = wavread("S.wav", format="native");
x, fs = convert(Vector{Float64}, vec(x)), convert(Int, fs);
winsize = 512;
hopsize = 256;
S = spectrogram(x, winsize, winsize-hopsize, fs=fs, window=hanning);
p = power(S);
q = log10.(p');
t = time(S);
f = freq(S);
@gp( "set nokey", "set grid", xlabel="Time [S]", ylabel="Freq [Hz]",
# "unset xtics", "unset ytics",
# "set xtics('3' 100, '6' 200, '9' 300)",
q, "with image")
# save(term="png", output="DSPspect.png")
wavreadな石とDSPな石とのマッチングを取る為、間にレベルコンバータ（じゃなくて型コンバータ）を挟んでいます。これが無いと、石から火もしくは発煙しちゃいます。（いかん、いかん、ついついハード寄りな表現になっちゃうな）
またspectrogramからの映像出力は位相がπ/2だけ回っているので、これを補正します。ああ、またハード的な表現だなあ。p' して、行列を転置しておきます。出力はパワー表現したいので、本当は、10 * log10 なんだけど、ごまかして、dB じゃなく B ベルさんにしています。グラハム・ベルさんに敬意を表してね。
出力結果を貼っておく。DSPspect.png
X軸、Y軸の数値が出鱈目（のように見える）。それには、下記のような理由がある。
julia> size(q)
(311, 257)
julia> t
0.032:0.032:9.952
julia> f
257-element Frequencies:
0.0
15.625
3984.375
4000.0
qと言うマトリックスの内容を画像とみなして表示されてるんだ。それじゃ理解に苦しむ（今みたいにね）。そんな事で、補足情報が得られるようになってる。X軸への補足情報はｔというレンジ型。Y軸は、周波数情報。これらを使って、上手に表示してねってのがDSP作者の願い。
でも、それを受け取ってgnuplotに反映するのは至難の技になる。ちと、その辺を調べてみたら 目盛見出し が出てきた。
set xticsとかを使って、地道にX軸の目盛り数字を再マップするのが正攻法。それとも、unsetxticsとかして、潔く表示を消しちゃうか。
そうよのう、soxiの情報を見れば、X軸、Y軸のスケールが分かるんだから、後は心眼でデータを読めよ。それでこそ達人さ。
修行を積んだ坊主は、心頭滅却すれば火もまた涼し と言うじゃないですか。最近の呆けた老人は、こう暑くても暑さを感じず、クーラーを入れないものだから、熱中症でバタバタと逝っている。自分では大丈夫と思っていても体がついていかないのだな。
range
上で、ｔはレンジ型とかやったけどこれを使って、xticsとかを書き換える事を考える。
julia> t
0.032:0.032:9.952
julia> typeof(t)
StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
julia> size(t) (311,) julia> t[1] 0.032 julia> t[311] 9.952 julia> t[end] 9.952 julia> t[100] 3.2 julia> t[150] 4.8
レンジ型と言っても、ベクターと同様に扱えるとな。
あれ、spectrogramからの出力は、Sで受けていたけど、octaveみたいに分解した形にならないのかな？
julia> p,t,f = spectrogram(x, winsize, winsize-hopsize, fs=fs, window=hanning);
ERROR: MethodError: no method matching iterate(::DSP.Periodograms.Spectrogram{Float64,Frequencies})
Closest candidates are:
iterate(::Core.SimpleVector) at essentials.jl:568
iterate(::Core.SimpleVector, ::Any) at essentials.jl:568
iterate(::ExponentialBackOff) at error.jl:199
怒られてしまった。I/Fが微妙に違うのね。
julia> S = spectrogram(x, winsize, winsize-hopsize, fs=fs, window=hanning); julia> S.time 0.032:0.032:9.952
こういうのは許されている。Sはコンポジット信号（色々な信号が重なった物を言う）、もとえ構造体なんだな。話が逸れた。
レンジ型の値を取り出すのは簡単だった。gnuplotはレンジ型の言わばインデックス番号をX軸なりの表示使ってる。ならば、レンジ型を参照して、その値を利用すれば、容易に見易い軸値になるだろう。
上のスクリプト中で、コメントになってた部分を、下記のように書き換える。
"set xtics('$(t[100])' 100, '$(t[200])' 200, '$(t[300])' 300)",
すると、X軸は、3.2, 6.4, 9.6 と表示が変わった。ちょっと端数が出ててかっこわるい。ああ、ここに出てきた、$( .. ) は、変数値のその場展開、インターポーレーションの書式ね。
julia> t[95] 3.04 julia> t[94] 3.008 julia> t[93] 2.976
ぴったりの数字はおぼつかないので、丸めは必要だろう。更に、幾つに分解したらいいかな？ ざっと、4分割ぐらいが、くどくなくていいのかな。
julia> 311/4 77.75 julia> t[78] 2.496
小数点以下2桁で丸めて、2.5にする方法を考えねば。なんかこんな事、現役時代にgnuplotと付き合ってて、真剣に実装した覚えがあるぞ。
conv
DSPのパッケージを入れると、もれなく畳み込みが出来るようになる。普通はこの畳み込み関数 conv を信号処理に使うけど、数式の畳み込みにも使える。
例として、(x + 1)(x - 1) を畳んでみる。式の係数だけをベクターにする。(x + 1) なら、(1 * x + 1) とみなして、ベクター [1,1]を得る。そんな具合で、畳んでみる。
julia> conv([1,1], [1,-1])
3-element Array{Int64,1}:
1
0
-1
係数は、[1,0,1] なので、(1 * x^2 + 0 * x -1)となる。簡約すると (x^2 - 1) こうなるね。どうやって実現してる？ ソース嫁。
"""
conv(u,v)
Convolution of two vectors. Uses FFT algorithm.
"""
function conv(u::StridedVector{T}, v::StridedVector{T}) where T<:BLAS.BlasFloat
nu = length(u)
nv = length(v)
n = nu + nv - 1
np2 = n > 1024 ? nextprod([2,3,5], n) : nextpow(2, n)
upad = [u; zeros(T, np2 - nu)]
vpad = [v; zeros(T, np2 - nv)]
if T <: Real
p = plan_rfft(upad)
y = irfft((p*upad).*(p*vpad), np2)
else
p = plan_fft!(upad)
y = ifft!((p*upad).*(p*vpad))
end
return y[1:n]
end
なんとFFTが顔を出してきたぞ。
octaveのconvはどうかと思って調べたら、主要部分はbuilt-inのconv2にお任せしてた。
octave:19> help conv2 'conv2' is a built-in function from the file libinterp/corefcn/conv2.cc
暑くて、やる気 Nothing.
plotsの市場調査
表示器と言うか可視化ツールをJulia Homeから見繕ってみた。
juliaだけで出来ているという表示器。表示はブラウザーに依頼。むき出しのブラウザー画面が出て来て、ちょっと原始的。
流行りのjavascriptで描画するって事は、これもブラウザーに依頼の口だな。変に自前の表示器を用意するより、ブラウザーに任せてしまうってのは妙案だと思うぞ。
julia> @time using PlotlyJS 10.277952 seconds (12.24 M allocations: 629.339 MiB, 3.34% gc time) julia> @time plot(1:10) ┌ Warning: Accessing `scope.id` is deprecated, use `scopeid(scope)` instead. │ caller = ip:0x0 └ @ Core :-1 7.846489 seconds (13.56 M allocations: 674.209 MiB, 4.47% gc time) julia> @time plot(1:20) 0.147175 seconds (6.71 k allocations: 423.750 KiB)
ぎりぎり待てるかな。表示画面が結構モダンな感じ。ひょっとして自前で用意した表示器かしら？ 疑問に思ったら即確認。psしてみる。
[ham@cent ~]$ ps -xocmd CMD : emacs snd.jl julia /home/ham/.julia/packages/Blink/AO8uN/deps/atom/electron /home/ham/.julia/packages/Blink/AO8uN/src/AtomShell/main.js port 5729 /home/ham/.julia/packages/Blink/AO8uN/deps/atom/electron --type=zygote --no-sandbox /home/ham/.julia/packages/Blink/AO8uN/deps/atom/electron --type=gpu-process --enable-features=SharedArrayBuffer --no-sandbox --gpu-prefere /home/ham/.julia/packages/Blink/AO8uN/deps/atom/electron --type=renderer --no-sandbox --enable-features=SharedArrayBuffer --service-pipe-t ps -xocmd
どうやら、自前で用意したっぽい。そのせいか重いブラウザーの起動がなく、心持ち速いのかな。
なんでも、元はd3.jsっていうjavascrit製の可視化ツールらしい。それをpythonの連中が何でもpythonでパッケージングしたとな。その上前にjuliaの連中が乗っかったと言う構図。
従って、
Plotly Python Open Source Graphing Library
matplotlib使いづらくない？plotlyで可視化しようよ
とかを見ておくと、幸せになれるでしょう。マウスでぐりぐりするにはうってつけです。
３次元構造のグラフを回転させて眺めるってのがストレスなく出来る。マウスをグラフ上に置けば、その点の数値情報がポップアップされて（ちとウザイけど）きて、未来型の表示器でしたよ。 素人受け、間違い無し。
記念に、上記のスペクトラムをPlotlyJSで、絵にしてみた。PlotlyJSでは、グラフを書くのは 関数にまとめるってのが流儀なのね。こういう風にまとめてしまうのは、後で見た時、すっきりしてて良いかも。
using PlotlyJS
function heatmap2(q)
trace = heatmap(
x=1:size(q)[1],
y=1:size(q)[2],
z=q
)
plot(trace)
end
heatmap2(q)
gnuplotの図と内容的にはほぼ同じはずなんだけど、PlotlyJSのカメラで撮影した図は、３倍も大きいなあ。何か余分についているのだろうか。
想像するに、マウスでポイントした時にポップ表示される生データが入ってるんかな。（掲載したものは、データがポップ表示されません。あくまでも建前はPNGなファイルですから）
etc
Linux development with WSL and Visual Studio Code
なかなか面白い。