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

(2020年8月15日:Julia 1.5で、REPLではlocalスコープ内からglobal変数を参照できるようになりました。Julia 0.6以前の挙動に戻ったことになります)

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

ピンクノイズJSFXからの出力のスペクトル

以下、JSFXコードです。

続きを読む

Gtk.jlの準備

(2020年2月15日追記:ここ半年ほどGtk.jlがインストールできない、あるいはインストールできてもウィンドウが表示されない状況が続いていました。さきほどJulia 1.3.1で再挑戦したところ、以下の対策をしなくても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年間ほとんど乗っていなかったことになる……。去年度はむちゃくちゃ忙しかったからな。今年は乗りたいが、どうなることか。

助言 —若さのように、若き日に無駄にしてしまうだけだろうが—

『助言 —若さのように、若き日に無駄にしてしまうだけだろうが—』というタイトルの以下の文章は、もともとMary SchmichによってChicago Tribune紙のコラムとして書かれたものです。Baz Lurhmannによって1999年に「Everybody's Free (to wear sunscreen)」(通称「Sunscreen Song」あるいは「Wear Sunscreen」)として発表されて有名になったもので、僕は留学先の寮の自室でMTVだかVH1だかで流れていたこの曲を聞いて感動して、すぐにCDを買ったけど歌詞カードがついてなくて、でも原作者のメールアドレスが載っていたので検索かけて、ようやく原作である『Advice, like youth, probably just wasted on the young - Chicago Tribune』にたどりついたのでした。

もちろん原文のほうが味わい深いし、当時の僕の語学力での翻訳なので怪しい部分もけっこうあります。今だったらこんな訳はしないだろうと思うところも沢山あります。だけど、僕が学生だった1999年5月3日(月)の作業を記録するためだけに掲載しておきます。

助言 —若さのように、若き日に無駄にしてしまうだけだろうが—

メアリ・シュミック/シカゴ・トリビューン紙 1997年6月1日

 97年卒業の皆さん、

 日焼けどめをつけなさい。

 君たちにひとつだけ将来へのちょっとした秘訣を教えるのであれば、日焼けどめがそれになるであろう。日焼けどめの長期間にわたる効果は科学者たちによって証明されているが、これから話す残りのことは私のとりとめのない経験によるものでしかない。いくつかの助言をしたいと思う。

 若さの力と美しさを存分に楽しみなさい。いや・・・若さの力と美しさなどというものは、それを失うまで理解できるものではない。しかし二十年後、自分の写真を見て、今ではもう手にすることのできない眼前に広がる途方もない可能性、そして自分の見た目がどんなに素晴らしかったか・・・そう思うことになる。君たちは、自分が想像するほど太ってはいない。

 将来のことを心配してはいけない。心配してもいいが、心配するのはガムをかむことで数式を解くのと同じ程度の効果しかないことを知っておいた方がいい。あなたがたの人生の大きな問題は、暇な火曜日の午後四時にいきなりやってくるような、頭をよぎりもしなかったようなことである場合が多い。

 毎日何かひとつ、自分を怯えさせることをやりなさい。

 歌いなさい。

 他人の心を傷つけてはいけない。あなたの心を傷つけるような人間と一緒にいることもない。

 歯を磨きなさい。

 嫉妬で時間を無駄にしてはいけない。前に出られるときもあるし、追いかけないといけないこともある。競争は長く、そして最後には、自分ひとりだけになる。

 受けた賛美の言葉を忘れないように。受けた侮辱は忘れなさい。うまく出来る方法を見つけたなら、私にも教えてほしい。

 昔の恋文はとっておきなさい。古い預金通帳は捨てなさい。

 伸びをしなさい。

 自分の人生をどうしたいかわからなくても、負い目を感じることはない。私が知っている、充実した人生を送っている人々は、二十二歳の時には自分の人生をどうしたいかなんてわかっていなかったし、私が知っている、とても充実した人生を送っている四十歳になる人々は、まだ人生をどうしたいかわかっていない。

 たっぷりとカルシウムをとりなさい。ひざをいたわりなさい。ひざが動かなくなったら、ひざが恋しくなる。

 あなたたちは結婚するかもしれないし、しないかもしれない。子供ができるかもしれないし、できないかもしれない。四十で離婚するかもしれないし、結婚七十五周年でいかしたダンスを踊るかもしれない。なにをするにせよ、自分を誉めすぎても叱りすぎてもいけない。自分のした選択はいつも半々だし、だれにとってもそうだ。

 自分のからだを楽しみなさい。できる限り使ってやりなさい。使うことを恐れることはないし、他人がどう思うかを心配することもない。あなたの体は、あなたが手に入れられる最高の道具なんだから。

 自宅の茶の間でしかできなくても、踊りなさい。

 書いてある通りにやらなくてもいいが、とにかく説明書は読みなさい。

 美容雑誌は自分自身を醜く感じさせるだけだから、読んではいけない。

 親のことを知りなさい。いつ親がいなくなるかなんて分からないから。兄弟とは仲良くしなさい。兄弟は君たちと過去をつなぐ最良の鎖で、きっとこれからも一緒にいることになるだろうから。

 友達は来て、そして去っていくものだと理解しなさい。でも、その中の大切な幾人かはしっかりつかんでおきなさい。地理的条件や生活習慣の違いを埋めるために、できるだけのことをしなさい。年をとるにつれ、若かったころのあなたを知る人が必要になってくる。

 一度ニューヨークに住むといい。でも自分が堅くなる前に離れること。北カリフォルニアにも住むといい。でも自分が軟らかくなる前に離れること。旅をしなさい。

 変わることのない真実は受け入れなさい。物価は騰がり、政治家は愛人をつくり、そして、あなたたちも歳をとる。そして歳をとったとき、自分が若かったころを美化するようになる。物価は妥当だったし、政治家は高潔で、そして子供たちは年上を尊敬していた。

 年上を尊敬しなさい。

 他人に支えてもらおうと期待してはいけない。信用基金が手に入るかもしれない。玉の輿に乗れるかもしれない。しかし、どちらもいつ打ち止めになるかは誰にも分からない。

 髪をいじりすぎてはいけない。さもないと君が四十になるころには、八十五のようになってしまう。

 だれの助言を受け入れるかは慎重に考えなさい。でも、助言してくれる人の話は辛抱強く聞くこと。助言は懐古のひとつの形なのだ。助言をすることは、ゴミの山から過去を拾ってきて丁寧に汚れを拭ってやり、汚い部分はペンキを塗りなおし、本来の価値以上に使うことに似ている。

 しかし、日焼けどめの効果についてだけは覚えておいてほしい。

楽器音の分析とスペクトログラム (Julia Advent Calendar 2018)

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

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

楽器音の分析

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

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

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

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

フーリエ変換では分析にある程度の時間区間を必要とし、その時間の長さが分析できる周波数分解能に関わってきます。たとえば、4 Hzごとの波の振幅が知りたいときには、1 ÷ 4 = 0.25 秒の時間が必要です。一方で周波数成分の0.01秒単位の時間変化を見ようと思った場合には、 1 ÷ 0.01 = 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では下のような再生可能なインターフェイスが表示されます(下は画像ですので再生できません)。

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)

スペクトログラムを作成してみましょう。以下のように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)")

これが音についての一般的なスペクトログラムで、横軸に時間、縦軸に周波数、色でパワーを表します。明るい部分に強いパワーがあるので、たとえば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