線形予測符号(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後に正負のエネルギを合算するところは本当に必要なのか分かりません)。「ここが間違ってるよ」「こうした方がいいよ」といった意見ありましたら、コメント欄からでも連絡ください。

Juliaの変数スコープについて

Juliaの変数スコープMatlabC言語と少し違うので混乱した話。(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);
}

totalfor文の内側から更新することができます。

あるいは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を超えたこともあったんだけど。

REAPERのJSFXでピンクノイズを生成してみた

REAPERのEffectsフォルダにはIvanov氏が作った「Pink Noise Generator」というJSFXが入っていたのですが、ソースを読んでみるとホワイトノイズにフィルタをかけてピンクにしているものでした。ホワイトノイズにフィルタをかけるのではなくピンクノイズを直接生成するプログラムを、以前Pythonで作ったことがあるので、それをリアルタイムで動作するようにしてみます。

marui.hatenablog.com

Voss-McCartney法そのままで、効率化もまったくしていません。しかも定数をそのまま書いているので汚いコードですが、低域から高域まできれいなピンクノイズが出力されます。

f:id:amarui:20190704141710p:plain
ピンクノイズJSFXからの出力のスペクトル

以下、JSFXコードです。

続きを読む

Gtk.jlの準備

JuliaでGUIを作りたいと思ってJulia Observerに行くと、いくつかの選択肢が掲載されています。掲載されているものの他にもWebIO.jlGtk.jlもありますし、以前はFLTKバインディングがあったような気もします。また、SDLを使うためのSimpleDirectMediaLayer.jlもあります。その中でも、ImageView.jlなどで使用されているGtk.jlを試してみたいと思います。

まずはインストールできなくてハマったところです。

Gtk.jlを動かすためにはCairo.jlが必要なのですが、macOS 10.14.5でCairo.jlがインストールできないという問題が発生しました。このことについてはGitHubの以下のIssueに報告されていますが、日本語でなく、また結構長いので、このエントリにも書き残しておきたいと思います。

github.com

僕の環境はmacOS 10.14.5、Julia 1.1.1で、Homebrewを使って様々なソフトを入れています。PkgモードでCairo.jlをaddするところまでは行ったものの、以下のような(上記Issueのここからコピペ)エラーが出てビルドができない状態でした。Cairo.jlがインストールできないのでGtk.jlも使えず、それに依存するImageView.jlもMakie.jlもprecompileできない状態です。

(v1.1) pkg> build Cairo
  Building LibCURL ─→ `~/.julia/packages/LibCURL/khRkY/deps/build.log`
  Building WinRPM ──→ `~/.julia/packages/WinRPM/Y9QdZ/deps/build.log`
  Building Homebrew → `~/.julia/packages/Homebrew/s09IX/deps/build.log`
  Building Cairo ───→ `~/.julia/packages/Cairo/CXPG1/deps/build.log`
┌ Error: Error building `Cairo`:
│
│ signal (11): Segmentation fault: 11
│ in expression starting at ~/.julia/packages/Cairo/CXPG1/deps/build.jl:165
│ _platform_strcmp at /usr/lib/system/libsystem_platform.dylib (unknown line)
│ Allocations: 19355595 (Pool: 19352771; Big: 2824); GC: 43
└ @ Pkg.Operations /Users/osx/buildbot/slave/package_osx64/build/usr/share/julia/stdlib/v1.1/Pkg/src/Operations.jl:1075

Issue上では原因究明のためにいろいろな試行錯誤がされていましたが、その結果、Cairo.jl 0.5.6やCairo.jl 0.6.0はglib 2.60とは相性が悪いというのが直接の原因だと判明しました。解決策としては、使用しているglibを0.58.3にすればよいとのこと。ただ、僕の環境ではすでにシステム側にもHomebrewが入っているので、Juliaが.juliaディレクトリ内にインストールするHomebrewと両方のglibのバージョンを下げる必要がありました。

https://github.com/JuliaGraphics/Cairo.jl/issues/271#issuecomment-480203117に書いてあるとおり、まずはJulia内Homebrewにglib 2.58.3が入っているか確認します。

julia> using Homebrew
julia> Homebrew.brew(`switch glib 2.58.3`)

ここで「Error: glib does not have a version "2.58.3" in the Cellar.(2.58.3はないよ)」と言われたので、2.58.3のソースコードをダウンロードしてコンパイル&インストールします。

julia> Homebrew.brew(`install --verbose --build-from-source https://raw.githubusercontent.com/Homebrew/homebrew-core/05871cb0394f78ef25a5c1c071456d0f1e4be4fe/Formula/glib.rb`)
julia> Homebrew.brew(`switch glib 2.58.3`)

すると、以下のように2.58.3がデフォルトで使用されるバージョンに指定されていることが確認できました。

julia> Homebrew.brew(`info glib`)
glib: stable 2.60.3
Core application library for C
https://developer.gnome.org/glib/
/Users/xxxxx/.julia/packages/Homebrew/s09IX/deps/usr/Cellar/glib/2.58.3 (451 files, 18.8MB) *
  Built from source on 2019-06-01 at 11:24:55
/Users/xxxxx/.julia/packages/Homebrew/s09IX/deps/usr/Cellar/glib/2.60.3 (429 files, 15.3MB)
  Poured from bottle on 2019-06-01 at 11:12:15
From: https://github.com/Homebrew/homebrew-core/blob/master/Formula/glib.rb
==> Dependencies
Build: meson ✘, ninja ✘, pkg-config ✔
Required: gettext ✔, libffi ✔, pcre ✔, python ✔
Process(`/Users/xxxxx/.julia/packages/Homebrew/s09IX/deps/usr/bin/brew info glib`, ProcessExited(0))

これでglibバージョンについては解決したと思いきや、さらにシステム側のHomebrewにインストールされているglib 2.60.3からも影響があるようで、まだビルドが通りません。そこでターミナルからも上記と同様のことをやります。

$ brew install --verbose --build-from-source https://raw.githubusercontent.com/Homebrew/homebrew-core/05871cb0394f78ef25a5c1c071456d0f1e4be4fe/Formula/glib.rb
$ brew switch glib 2.58.3

無事にシステム側のHomebrewでもglib 0.58.3が使える状態になりました。

$ brew info glib
glib: stable 2.60.3 (bottled)
Core application library for C
https://developer.gnome.org/glib/
/usr/local/Cellar/glib/2.58.3 (451 files, 18.8MB) *
  Built from source on 2019-06-01 at 11:38:14
/usr/local/Cellar/glib/2.60.3 (429 files, 15.3MB)
  Poured from bottle on 2019-05-22 at 16:33:03
From: https://github.com/Homebrew/homebrew-core/blob/master/Formula/glib.rb
==> Dependencies
Build: meson ✘, ninja ✘, pkg-config ✔
Required: gettext ✔, libffi ✔, pcre ✔, python ✔
==> Analytics
install: 267,265 (30 days), 780,837 (90 days), 2,035,504 (365 days)
install_on_request: 12,612 (30 days), 36,289 (90 days), 89,127 (365 days)
build_error: 0 (30 days)

その後、Julia上で(v1.1) pkg> build Cairoを行うと無事にCairoのインストールが完了し、次いでGtk.jlもインストールできました。

Juliaプロンプト上でGtk.jlのGetting Startedにあるサンプルを実行すると、無事にウィンドウとボタンが表示されました。よかった。

using Gtk
win = GtkWindow("My First Gtk.jl Program", 400, 200)
b = GtkButton("Click Me")
push!(win, b)
showall(win)

f:id:amarui:20190608132252p:plain