PortAudio.jlが公式レジストリから消えてしまったのか、これまでのようにインストールしようとすると「見つからない」と言われてしまいます。
(v1.2) pkg> add PortAudio#julia1 ERROR: The following package names could not be resolved: * PortAudio (not found in project, manifest or registry) Please specify by known `name=uuid`.
GitHubにはこれまでどおりコードが置いてありますので、
add https://github.com/JuliaAudio/PortAudio.jl#julia1
とURLを指定することで直接GitHubのレポジトリからインストールすることができます。
10月16日(水)〜10月19日(土)に行われた147th AES Conventionに行ってきました。AES Conventionは大きく製品展示と技術発表に分かれます。僕は技術発表を中心に参加して30件以上の発表を見聞きしましたが、その中で個人的に面白いと思った発表をいくつか書き残しておこうと思います。技術発表の中にもスタジオや放送の現場に近い内容を扱う発表もあれば、より基礎的な研究に焦点を絞った発表もあります。以下、どちらかというと研究方面に偏っています。
続きを読む学生のときに受けた講義に「プログラミング」がありました。「プログラミングI」「同II」「同III」と半期の科目が3つあり、それぞれ「プログラミングI」は「基本のC言語」、「プログラミングII」は「色々経験するために、関数型言語（LISP）・論理型言語（Prolog）・並列並行処理（AEGISという言語だったような？）」、「プログラミングIII」は「（だいぶ昔のことなのでうろ覚えですが）今後はオブジェクト指向プログラミングが主流になるからC++とJava」というようになっていた気がします。
「プログラミングII」で関数型言語に触れる部分ではLISPを習いました。いま思えばMLとかでも良かったかと思いますが、当時はLISPは実用性のあるメジャーな言語という雰囲気だったのかもしれません。歴史的にも勉強しておいた方がいい言語の一つですしね。
LISPで最初につまづいたのがリスト操作です。いまから思い返すとなぜ理解できないのか分からない、なんて気がしますが、当時はcar、cdr、consでのリスト操作に違和感があり、何が何だか分かりませんでした。おそらくcar+cdrとconsが対称な操作であることが分かっていなかったのではないかと思います。
たとえば
(car '(1 2 3)) -> 1 (cdr '(1 2 3)) -> (2 3) (cons 1 '(2 3)) -> (1 2 3)
について、carはリストの最初の要素をとってくる、cdrはリストの残りをリストとして得る、consは引数2の先頭に引数1を挿入するという動作をするので、当時の自分の理解を図解すると以下のような感じです。
この簡単なことが当時の僕にはよく分からず、
(cons '(1 2) '(3 4)) -> (1 2 3 4)
を期待しているのに
(cons '(1 2) '(3 4)) -> ((1 2) 3 4)
となるので「先頭に値を入れてくれるのではないのか？」などと混乱していました。同じ頃にMatlabも使い始めていたので、
>> A = [1 2];
>> B = [3 4];
>> [A B]
ans =
1 2 3 4
と同じような動作を期待していたのだと思います。
LISPに触れるのは半期の授業の1/3だけなので4回程度。話題はあっという間に他の言語に移ってしまったので「論理型言語おもしろいけど、家系図つくって実用性あるのか？」だの「並列・並行は意味が分からない、円卓の哲学者は馬鹿か？」などと言い始め、期末試験が終わるとそれら全てを忘れてしまい、卒業までずっとC言語を使っていました。
その後、何年も経ってからLISPに再入門すると、car+cdrとconsの見え方が変わっていました。値やリストを抽出・挿入という考え方から、括弧の位置が動くという考え方になりました。図解すると以下のような感じです。
carやcdrは括弧が1つ右に動く感覚です。たとえば(car '(1 2 3))は(1 2 3)の先頭の括弧が1 (2 3)と1つぶん右に移動し、そのときに括弧外に出たものがcarで、括弧内に残っているものがcdrです。一方でconsは括弧が1つ左に動く感覚です。リストは1つのまとまりなので、(cons '(1 2) '(3 4))をすると、(3 4)の先頭の括弧が(1 2)のまとまりを超えて左側に動き((1 2) 3 4)となるのでした。
結局、(1 2)と(3 4)から僕が期待したとおりの(1 2 3 4)を作るためには、(1 2)の各要素を取り出しつつconsを使い
(cons (car '(1 2)) (cons (car (cdr '(1 2))) '(3 4))) -> (1 2 3 4)
とすることになります。（これをするために(append '(1 2) '(3 4))が用意されています）
StautnerとPucketteが提案したFeedback Delay Networkという人工リバーブがあります。このFDNを実装するときに必要になるものに、N×Nの行列があります。これを行列Gとしたとき、各要素g_{i,j}はディレイラインiからディレイラインjへのフィードバック・ゲインが入っています。この行列はいろいろと選択肢があるのですが、J. Smithによると「アダマール行列も使えるよ」ということでした。
MaxやPureDataでFDNを実装するときにアダマール行列が欲しかったので、Juliaで計算したものを入力して使おうと思いました。ただ、MatlabやOctaveにあるようなhadamard関数がJuliaにはなかったので、作りました。参考にしたのはWikipediaに載っていた“他の生成法”。MatlabだとN=2k * pでp=1, 12, 20, 28のときの計算ができるけど、とりあえず今回はp=1のときだけ。
""" hadamard(n) Compute the Hadamard matrix of size `n`-by-`n`. `n` must be a positive integer of power of 2. """ function hadamard(n) function hadamardF(n) if n <= 1 return [0 1]; else return [zeros(1, 2^(n-1)) ones(1, 2^(n-1)) hadamardF(n-1) hadamardF(n-1)]; end end F = hadamardF(floor(Int, log2(n))); FF = F' * F; FF[(FF .% 2) .== 1] .= -1; FF[(FF .% 2) .== 0] .= 1; return map(Int, FF); end
以下のように使えます。
julia> hadamard(4)
4×4 Array{Int64,2}:
1 1 1 1
1 -1 1 -1
1 1 -1 -1
1 -1 -1 1
ただ、アダマール行列を計算するプログラムはすでに色々とあります。上記のコードよりもっと高機能なものが以下のリンク先にあるので、それを使うのもよいですね。
Juliaを使った音の処理については、これまでにいくつかのエントリを書きましたが、今日は線形予測符号（LPC; Linear Predictive Coding）を用いてスペクトル包絡を求めてみます。ここでは線形予測符号の説明は省きますが、音声信号処理の分野では昔からよく使われている方法で、Juliaコードを書くにあたっては以下のページを参考にしています。
なお、JuliaでLPCを自作することはせず、DSPパッケージに入っているlpc()をありがたく利用します。
まず、使用するパッケージ群の読み込みを行います。
using FileIO: load import LibSndFile using DSP using FFTW using Plots
LibSndFileを使って音ファイルを読み込んできます。今回は音声ではなく、とある管楽器の音を使ってみます。
snd = load(expanduser("~/instrument.wav")) x = snd.data[1:60000]; fs = snd.samplerate; t = [0:length(x)-1] / fs; plot(t, x, ylim=[-1, 1], xlabel="Time (s)", ylabel="Amplitude", legend=false)
楽器音を読み込んで、そこから音信号xと標本化周波数fsを取り出します。また波形のプロットに使いたいので、xの長さに合わせた各サンプルの時刻tを計算します。plot()で波形表示をしています。
音声の場合はフォルマントを把握しやすくするために高域を強調する「プリエンファシス・フィルタ」をかけるそうなのですが、今回は楽器音を使っているので、そのままにしておきます。
パワースペクトルの計算は必須ではありませんが、スペクトル形状を確認したいので計算しています。
y = fft(x) / length(x); yy = abs.(y); yy = yy[1 : div(length(y), 2)+1]; yy[2:length(yy)-1] .*= 2; f = range(0, fs/2, length=length(yy)); plot(f, 20*log10.(yy), alpha=0.5, label="Power Spectrum", xlabel="Frequency (Hz)", ylabel="Power (dB)")
フーリエ変換（fft）をして、正の周波数領域の振幅スペクトルだけを抜き出す処理をしています。
各行をなぞっていきます。まずabs()で振幅スペクトルを得て（2行目）、DCからナイキストまでの成分を抜き出します（3行目）。すると負の周波数側の成分を捨てたことになるので2倍して補正します（4行目）。yyは各周波数ビンの振幅が入っているので、各ビンの中心周波数を求めます（5行目）。
プロットでは、あとでLPCの結果を重ねたいので、alphaを指定して振幅スペクトルの色を薄くしています。
ここからが本番です。DSPパッケージにlpc()が入っているので、それを使います。LPCの次数は適当に決めて、20次でやってみます。lpc()を使うとフィルタ係数aが求められます。
p = 20 a, e = lpc(float(x), p); env_x = filt([1], [1; a], [1; zeros(length(x)-1)]); env_y = fft(env_x) / length(env_x); env_yy = abs.(env_y); env_yy = env_yy[1 : div(length(env_y), 2)+1]; env_yy .*= 2; plot!(f, 20*log10.(env_yy), width=1.5, label="Spectral Envelope (LPC p=$(p))")
フィルタ係数からエンベロープの形状を計算するのにDSPパッケージのfreqz()を使えることをすっかり忘れていたので、上記ではフィルタにインパルスをつっこんで周波数応答を求めています。
他にもケプストラムを使う方法もありますが、そちらについてはまた別の機会に試してみたいと思います。
上記の方法は厳密に正しくない部分があるかもしれません（たとえばFFT後に正負のエネルギを合算するところは本当に必要なのか分かりません）。「ここが間違ってるよ」「こうした方がいいよ」といった意見ありましたら、コメント欄からでも連絡ください。
（2020年8月15日：Julia 1.5で、REPLではlocalスコープ内からglobal変数を参照できるようになりました。Julia 0.6以前の挙動に戻ったことになります）
Juliaの変数スコープはMatlabやC言語と少し違うので混乱した話。（2019年8月25日：global/localスコープの説明が公式ドキュメントと違ったので修正しました）
たとえば、1から10までの整数の和を求めたいとき、C言語では以下のようにやるかと思います。
#include <stdio.h> int main(int argc, char *argv[]) { int total = 0; int x; for (x=1; x<=10; x++) { total += x; } printf("Total: %d\n", total); return(0); }
totalをfor文の内側から更新することができます。
あるいはPython 3のREPLでも以下のようになります。forの内側から、外側にあるtotalにアクセスできていることが分かります。
>>> total = 0 >>> for x in range(1, 11): ... total += x ... >>> total 55
同じことをJulia 1.1でやってみます。ここでは（立ち上げたばかりのREPLなど）トップレベルでの実行を想定しています。
julia> total = 0; julia> for x = 1:10 total += x; end ERROR: UndefVarError: total not defined Stacktrace: [1] top-level scope at ./REPL[2]:2 [inlined] [2] top-level scope at ./none:0
公式ドキュメントによると、for文は新しいlocalスコープを作るため、内部（localスコープ）から外側（globalスコープ）の変数を見ることができないとのことです。回避するにはglobalを使うらしく、
total = 0; for x = 1:10 global total += x; end println("Total: $(total)")
とすれば期待通りの動作になります。globalというキーワードを付けてあげれば上位にあるglobalスコープの変数にもアクセスできるのだな、ふむふむ、という感じです。
ただ、問題はここからで、二重のループにするとそういったことが起こりません。for yの中から、その外側にあるはずのsubtotalを更新できます。
for x = 1:10 subtotal = 0; for y = 1:10 subtotal += y; end end
これは（Local Scopesに書いてありますが）、localスコープ内に作られたlocalスコープからは上位の変数にアクセスできるからなのだそうです。外側のforがlocalスコープを作ったので、その内側のfor文内部から見ると上位がlocalスコープなので、for yからfor x内の変数にアクセスできるのです。
同じものを関数にした場合はどうなるでしょうか。
function testsum() total = 0; for x = 1:10 total += x; end return(total); end
この場合は、functionが新しいlocalスコープを作るので、その内側に作られたlocalスコープのforから見て上位にあるtotalにはアクセスできます。
globalはキモい（し、不具合の温床になる）のでできるだけ使わない方がよさそう、と考えると、Julia言語としては「できるだけ関数にしよう」という方向にユーザーを誘導していますのでるのではないかと感じます。そういえば、関数化した方が実行も速くなることが多いです。（追記：公式スタイルガイドに「スクリプトではなく関数にせよ」と書いてありました）
同じような操作であっても言語ごとに少しずつ違うので、細かいところで「どうなっていたっけ？」とドキュメントを参照する必要がありますが、こういうところから言語設計者の考え方に思いをはせるのもプログラミングの楽しいところです。
| 給油日 | オドメーター (km) | 給油量 (L) | 単価 (円/L) | 燃費 (km/L) | 距離単価 (円/km) |
|---|---|---|---|---|---|
| 2017-09-15 | 10004.7 | 3.09 | 123.96 | 53.03 | 2.34 |
| 2017-10-04 | 10152.2 | 2.08 | 125.96 | 70.91 | 1.78 |
| 2017-10-18 | 10314.9 | 3.42 | 125.73 | 47.57 | 2.64 |
| 2018-05-05 | 10499.5 | 3.18 | 137.74 | 58.05 | 2.37 |
| 2019-06-08 | 10607.9 | 2.40 | 148.75 | 45.17 | 3.29 |
| 2019-07-13 | 10753.3 | 3.20 | 137.81 | 45.44 | 3.03 |
日野まで行った帰りに給油。街乗りが多かったので燃費が悪い……。いちばん良いときは70 km/Lを超えたこともあったんだけど。
