PortAudio.jlのインストール(2019年11月3日時点)

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のレポジトリからインストールすることができます。

AES147 New York技術発表より

10月16日(水)〜10月19日(土)に行われた147th AES Conventionに行ってきました。AES Conventionは大きく製品展示と技術発表に分かれます。僕は技術発表を中心に参加して30件以上の発表を見聞きしましたが、その中で個人的に面白いと思った発表をいくつか書き残しておこうと思います。技術発表の中にもスタジオや放送の現場に近い内容を扱う発表もあれば、より基礎的な研究に焦点を絞った発表もあります。以下、どちらかというと研究方面に偏っています。

続きを読む

carとcdrとcons

学生のときに受けた講義に「プログラミング」がありました。「プログラミングI」「同II」「同III」と半期の科目が3つあり、それぞれ「基本のC言語」「色々経験するために、関数型言語LISP)・論理型言語(Prolog)・並列並行処理(AEGISという言語だったような?)」「(だいぶ昔のことなので)今後はオブジェクト指向プログラミングが主流になるから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を挿入するという動作をするので、図解すると以下のような印象でした。

f:id:amarui:20191006112409p:plain

この簡単なことが当時の僕にはよく分からず、

(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の見え方が変わっていました。値やリストを抽出・挿入という考え方から、括弧の位置が動くという考え方になりました。図解すると以下のような感じです。

f:id:amarui:20191006114038p:plain

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))が用意されています)

アダマール行列(N = 2^kのときのみ)

StautnerとPucketteが提案したFeedback Delay Networkという人工リバーブがあります。このFDNを実装するときに必要になるものに、N×Nの行列があります。これを行列Gとしたとき、各要素g_{i,j}はディレイラインiからディレイラインjへのフィードバック・ゲインが入っています。この行列はいろいろと選択肢があるのですが、J. Smithによると「アダマール行列も使えるよ」ということでした。

MaxやPureDataでFDNを実装するときにアダマール行列が欲しかったので、Juliaで計算したものを入力して使おうと思いました。ただ、MatlabOctaveにあるような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

ただ、アダマール行列を計算するプログラムはすでに色々とあります。上記のコードよりもっと高機能なものが以下のリンク先にあるので、それを使うのもよいですね。

nemocas.org

github.com

線形予測符号(LPC)でスペクトル包絡を計算する

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()で波形表示をしています。

f:id:amarui:20190814135831p:plain

音声の場合はフォルマントを把握しやすくするために高域を強調する「プリエンファシス・フィルタ」をかけるそうなのですが、今回は楽器音を使っているので、そのままにしておきます。

パワースペクトルの計算

パワースペクトルの計算は必須ではありませんが、スペクトル形状を確認したいので計算しています。

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を指定して振幅スペクトルの色を薄くしています。

f:id:amarui:20190814140424p:plain

LPC

ここからが本番です。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()を使えることをすっかり忘れていたので、上記ではフィルタにインパルスをつっこんで周波数応答を求めています。

f:id:amarui:20190814142702p:plain

おわりに

他にもケプストラムを使う方法もありますが、そちらについてはまた別の機会に試してみたいと思います。

上記の方法は厳密に正しくない部分があるかもしれません(たとえばFFT後に正負のエネルギを合算するところは本当に必要なのか分かりません)。「ここが間違ってるよ」「こうした方がいいよ」といった意見ありましたら、コメント欄からでも連絡ください。