本日の給油

給油日 オドメーター (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

本日の給油

給油日 オドメーター (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

過去5回分の給油だけ記載。2018年は1回だけ給油して、それから1年間ほとんど乗っていなかったことになる……。去年度はむちゃくちゃ忙しかったからな。今年は乗りたいが、どうなることか。

楽器音の分析とスペクトログラム

この記事はJulia Advent Calendar 2018の24日目です。

僕は日頃、おもにオーディオ信号を扱うのにJuliaを使用しています。今日はそのうち楽器音の分析のために作成したプログラムを紹介します。

楽器音の分析

楽器の音を聞いたとき、トランペットとピアノとでは同じ大きさ・同じ高さで鳴らされた音であっても異なる楽器だと区別できます。これは音色が異なるからです。音色は「音の大きさと高さをそろえても、なお異なる聞こえ方をするとき、その違いに相当するもの」というような定義がされています。その音色に対応する要因としては、含まれる成分音のレベルや時間変化などが考えられます。

楽器音に含まれている振動を(フーリエ変換などを利用して)周波数ごとに取り出すと「それぞれの周波数において、どのくらいの強さの音が含まれているのか」が分かります。これを周波数スペクトル(縦軸をどのようにとるかで振幅スペクトル、パワースペクトルなどとも)と言います。成分音の中で低い周波数に突出しているのが基音で、ピッチ感のある音については基音周波数の整数倍の周波数に倍音(あるいは高調波成分)があることがほとんどです。多くの場合は基音の周波数に対応したピッチを感じます(基音が最も強いエネルギーを持っている必要はありません)。打楽器のようにピッチ感が曖昧な音は、整数倍の音列になっていないことが原因のようです。

f:id:amarui:20181224002026p:plain

同じような音が続く(定常的な音の)場合は周波数分析をして倍音がどのように含まれているかを見れば良いのですが、楽器の音の場合は弾き始めの音の立ち上がりや減衰部分など、様々な時間変化があります。そのとき、時間に応じて倍音の構成なども変わる場合があります。そのため、時間波形の概形を見るためにエンベロープ(包絡線)を計算したりします。

時間波形全体の変化ではなく周波数成分の時間変化を見たいという時もあります。そこで登場するのがスペクトログラムです。主に短時間フーリエ変換を使用し、その名の通り短い時間範囲をフーリエ変換してその範囲のスペクトルを求めることを信号全体にわたって行うものです。下の図にその様子を描きました。横軸が時間、縦軸がパワー、手前から奥に向かって周波数としています。各時刻において、一定の時間幅についての周波数分析を行っており、その1フレームだけを取り出すと右側に取り出したようにスペクトルになっています。

f:id:amarui:20181224002056p:plain フーリエ変換では分析にある程度の時間区間を必要とし、その時間の長さが分析できる周波数分解能に関わってきます。たとえば、4 Hzごとの波の振幅が知りたいときには、1 ÷ 4 = 0.25 秒の時間が必要です。一方で周波数成分の0.01秒単位の時間変化を見ようと思った場合には、 1 ÷ 0.0.1 = 100 Hzの周波数幅でしか分析ができません。

そういった弱点はありますが、以下ではスペクトログラムを見ていきましょう。

Juliaでのスペクトログラムの計算

これ以降のコードは次の環境で実行したものです。

Julia Version 1.0.3
Commit 099e826241 (2018-12-18 01:34 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) i7-7920HQ CPU @ 3.10GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, skylake)

まずは必要になるパッケージ群をusingします。

# 信号処理に使う
using FFTW, DSP, SampledSignals

# WAVファイルの読み込みに使用
using FileIO: load
import LibSndFile

# グラフの作成はPyPlotを使う
using Plots
pyplot();

# 書式指定の文字式を作りたいときに使用
using Printf

LibSndFileを使ってWAVに保存された音ファイルを読み込みます。fileには拡張子まで入っていますが、拡張子を除いた部分をfilenameに入れて、これ以降で使用できるようにします。セミコロンなしでloadすると、音の概形が表示されます。コマンドラインでは実効値が、IJuliaでは下のような再生可能なインターフェイスが表示されます(下は画像ですので再生できません)。

f:id:amarui:20181224002150p:plain

file = "ching.wav";
filename = splitext(file)[1];
snd = load("$(filename).wav")

音オブジエクトには標本化周波数(samplerate)と音のサンプルデータ(data)が入っています。それぞれを取り出して、時間波形をプロットしてみます。0.2秒付近で立ち上がり、そのあと急速に減衰していく様子が見えます。一定間隔で振幅が波打って増減している部分がありますが、うなりが発生しているのでしょうか。

fs = snd.samplerate;
x = vec(snd.data[:,1]);
t = range(0, stop=length(x)-1) / fs;
plot(t, x, ylim=(-1, +1),
    xlabel="Time (sec)", ylabel="Amplitude", legend=false)

f:id:amarui:20181224002347p:plain

スペクトログラムを作成してみましょう。以下のようにDSP.jlのspectrogramを使用すると簡単です。引数は

  1. 分析対象の信号
  2. 分析窓長
  3. オーパーラップ量(ここでは窓の移動量(hopsize)を定義しておいてwinsize-hopsizeと指定しています)
  4. FFT長...窓あたりのフーリエ変換をどのくらいの長さで行うか(長く指定するとゼロパッドが行われます)
  5. 標本化周波数
  6. 窓の種類(今回はハニング窓)

の6つです。

# compute spectrogram
winsize = 4096;
hopsize = 256;
S = spectrogram(x, winsize, winsize-hopsize, nfft=winsize*2, fs=fs, window=hanning);

spectrogramからの出力(S)にはpowerfreqtimeが入っています。powerは、【freqの長さ】×【timeの長さ】の2次元行列で、各時刻・各周波数でのパワーが格納されています。$L=10\log_{10}\left(I\right)$を使ってデシベルに変換し、ヒートマップ表示してみます。

heatmap(S.time, S.freq, 10*log10.(S.power),
    xlabel="Time (sec)", ylabel="Frequency (Hz)")

f:id:amarui:20181224002457p:plain

これが音についての一般的なスペクトログラムで、横軸に時間、縦軸に周波数、色でパワーを表します。明るい部分に強いパワーがあるので、たとえば4000 Hz付近に横長に伸びる明るい線が基音を示しています。また、突出した成分音が整数次でない場所に出ていることも分かります(体鳴楽器に多く見られる特徴です)。

アニメーションを作る

上記のスペクトログラムを、一枚ずつ次々と表示していけばアニメーションさせることもできます。このほうが時間変化が直感的に把握しやすいという人もいるかもしれません。

以下で表示したい周波数範囲のみを取り出す処理をしています。楽器によって音域が違うので、分析対象の楽器にあわせてこの部分を調整します。今回は0 Hzから分析できる上限であるナイキスト周波数(標本化周波数の半分)までを分析するように指定しています。

# normalize the power (maximum will become 0 dBFS)
Spow = S.power ./ maximum(S.power);

# extract only the frequency range we need
fL = 0;
fH = fs/2;
Sf = S.freq[fL .<= S.freq .<= fH];
Spow = Spow[fL .<= S.freq .<= fH, :];
St = S.time;

アニメーションを作るのには@animateを使います。@animateつきのfor文の中で一枚一枚をプロットして、animに格納された画像をアニメGIFファイルとして保存します。gif()の引数にfps=15とあるのは1秒あたりのフレーム数(frames per second)で、あまり大きくない整数値であればOKのようです。

(実行すると結構時間がかかります! 自分の環境(MacBook Pro 2017, 3.1 GHz Intel Core i7)では2分ほどかかりました)

# create animation
@time anim = @animate for fr=1:length(St)
  timetext = @sprintf("%.2f", St[fr]);

  # plot spectrum
  plot(Sf, 10*log10.(Spow[:,fr]),
       xlim=(fL, fH), ylim=(-105, +5),
       title="$(filename), t=$(timetext) sec",
       xlabel="Frequency (Hz)", ylabel="Power (dB)",
       legend=false, size=(720, 480));
end

gif(anim, "tmp_$(filename).gif", fps=15);
run(`ffmpeg -loglevel quiet -y -r 24 -i tmp_$(filename).gif -r 24 -vcodec libx264 -pix_fmt yuv420p $(filename).mp4`);

なお、run命令を使って外部プログラムのffmpegを呼び出し、アニメーションGIFをMP4に変換しています。GIFだけでいい場合はrunの行を省略してください。完成したものを以下にアップロードしました。


蛇足

今回はDSP.jlのspectrogramを使用しましたが、スペクトログラムのプログラム(短時間フーリエ変換)を自作しようとすると以下のような感じになるでしょうか。ただし、簡易版なのでspectrogramのように機能豊富ではなく、必要最低限のものしか入っていません。参考までに。

function spect(x, fs, window_size=1024, hop_size=256)
  num_frames = ceil(Int, (length(x) + window_size-1) / hop_size);
  xx = [x ; zeros(2 * window_size)];
  N = nextpow(2, window_size);
  zeropad = zeros(N - window_size);
  y = zeros(Complex, round(Int, N/2+1), num_frames);
  for n=1:num_frames
    xxx = xx[hop_size*(n-1)+1 : hop_size*(n-1)+window_size];
    y[:,n] = rfft([xxx ; zeropad]) / N;
  end
  t = ((0:num_frames-1) * hop_size) / fs;
  f = range(0, stop=fs/2, length=round(Int, N/2+1));
  return(y, t, f);
end