| 給油日 | オドメーター (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の24日目です。
僕は日頃、おもにオーディオ信号を扱うのにJuliaを使用しています。今日はそのうち楽器音の分析のために作成したプログラムを紹介します。
楽器の音を聞いたとき、トランペットとピアノとでは同じ大きさ・同じ高さで鳴らされた音であっても異なる楽器だと区別できます。これは音色が異なるからです。音色は「音の大きさと高さをそろえても、なお異なる聞こえ方をするとき、その違いに相当するもの」というような定義がされています。その音色に対応する要因としては、含まれる成分音のレベルや時間変化などが考えられます。
楽器音に含まれている振動を（フーリエ変換などを利用して）周波数ごとに取り出すと「それぞれの周波数において、どのくらいの強さの音が含まれているのか」が分かります。これを周波数スペクトル（縦軸をどのようにとるかで振幅スペクトル、パワースペクトルなどとも）と言います。成分音の中で低い周波数に突出しているのが基音で、ピッチ感のある音については基音周波数の整数倍の周波数に倍音（あるいは高調波成分）があることがほとんどです。多くの場合は基音の周波数に対応したピッチを感じます（基音が最も強いエネルギーを持っている必要はありません）。打楽器のようにピッチ感が曖昧な音は、整数倍の音列になっていないことが原因のようです。
同じような音が続く（定常的な音の）場合は周波数分析をして倍音がどのように含まれているかを見れば良いのですが、楽器の音の場合は弾き始めの音の立ち上がりや減衰部分など、様々な時間変化があります。そのとき、時間に応じて倍音の構成なども変わる場合があります。そのため、時間波形の概形を見るためにエンベロープ（包絡線）を計算したりします。
時間波形全体の変化ではなく周波数成分の時間変化を見たいという時もあります。そこで登場するのがスペクトログラムです。主に短時間フーリエ変換を使用し、その名の通り短い時間範囲をフーリエ変換してその範囲のスペクトルを求めることを信号全体にわたって行うものです。下の図にその様子を描きました。横軸が時間、縦軸がパワー、手前から奥に向かって周波数としています。各時刻において、一定の時間幅についての周波数分析を行っており、その1フレームだけを取り出すと右側に取り出したようにスペクトルになっています。
フーリエ変換では分析にある程度の時間区間を必要とし、その時間の長さが分析できる周波数分解能に関わってきます。たとえば、4 Hzごとの波の振幅が知りたいときには、1 ÷ 4 = 0.25 秒の時間が必要です。一方で周波数成分の0.01秒単位の時間変化を見ようと思った場合には、 1 ÷ 0.01 = 100 Hzの周波数幅でしか分析ができません。
そういった弱点はありますが、以下ではスペクトログラムを見ていきましょう。
これ以降のコードは次の環境で実行したものです。
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を使用すると簡単です。引数は
hopsize）を定義しておいてwinsize-hopsizeと指定しています）の6つです。
# compute spectrogram winsize = 4096; hopsize = 256; S = spectrogram(x, winsize, winsize-hopsize, nfft=winsize*2, fs=fs, window=hanning);
spectrogramからの出力（S）にはpower、freq、timeが入っています。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
Juliaが1.0になって、徐々に環境が整ってきました。ここ数日、LibSndFile.jlやPortAudio.jlを使ってゴニョゴニョやっていたことをまとめておきます。
入出力の同期とか考えずに最低限の音再生ができればいいという場合は、
using FileIO: load, save import LibSndFile using PortAudio snd = load("hoge.wav"); stream = PortAudioStream(0, 2); write(stream, snd);
とすれば、ファイルを読み込んで、デフォルトの2チャンネル再生デバイス（内蔵スピーカとか）から再生してくれます。
続きを読む上記のJulia 0.6でやったことと同じですが、Julia 1.0がリリースされたので変化があったか見てみました。
前回と同じくadd1()を関数定義します。これはJulia 0.6でも1.0でも同じ。
julia> function add1(x)
return x+1
end
add1 (generic function with 1 method)
Julia 0.6では@code_llvmと@code_nativeの出力結果は以下のようになっていました。
julia> @code_llvm add1(10)
define i64 @julia_add1_62318(i64) #0 !dbg !5 {
top:
%1 = add i64 %0, 1
ret i64 %1
}
julia> @code_native add1(10)
.section __TEXT,__text,regular,pure_instructions
Filename: REPL[1]
pushl %ebp
decl %eax
movl %esp, %ebp
Source line: 2
decl %eax
leal 1(%edi), %eax
popl %ebp
retl
nop
nop
nop
nop
nop
nop
Julia 1.0では以下のとおり。
julia> @code_llvm add1(10)
; Function add1
; Location: REPL[1]:2
define i64 @julia_add1_35956(i64) {
top:
; Function +; {
; Location: int.jl:53
%1 = add i64 %0, 1
;}
ret i64 %1
}
julia> @code_native add1(10)
.section __TEXT,__text,regular,pure_instructions
; Function add1 {
; Location: REPL[1]:2
; Function +; {
; Location: REPL[1]:2
decl %eax
leal 1(%edi), %eax
;}
retl
nopw %cs:(%eax,%eax)
;}
コメントが挟まるようになり、何をしているかが分かりやすくなりました。LLVMコードは同じままですが、ネイティブコードの見た目がだいぶ変わりました。最近のアセンブリはわからないのですが、なんとなく雰囲気からEAX（レジスタ名）、DEC（1だけ引き算）、LEA（当該アドレスを返す）、RET（プロシージャから帰る）、NOP（何もしない）あたりだと考えられます。末尾にLがついているのは64 bit（Long）、Wは任意長（Words）でしょうか？ NOPを続けるのは一連の命令の合計長さを64 bitの倍数にしたいから？
ちゃんと理解しようと思うとIntelのマニュアルとか読むといいのかもしれませんね。
Intel® 64 and IA-32 Architectures Software Developer Manuals | Intel® Software
ヨーロッパLISPコンファレンスで発表された上記のプレゼンの中に、Juliaコードをネイティブコードに変換して中身を確認するマクロが紹介されていました。
公式ドキュメントではEssentials · The Julia Languageに書かれていますが、面白そうだったので、ちょっとだけ試してみました。
まずは簡単な関数を定義します。
julia> function add1(x)
return x+1
end
add1 (generic function with 1 method)
@code_llvmや@code_nativeマクロを使ってあげると、LLVMバイトコードや、ネイティブコードのアセンブリが見られます。
julia> @code_llvm add1(10)
define i64 @julia_add1_62362(i64) #0 !dbg !5 {
top:
%1 = add i64 %0, 1
ret i64 %1
}
julia> @code_native add1(10)
.section __TEXT,__text,regular,pure_instructions
Filename: REPL[5]
pushl %ebp
decl %eax
movl %esp, %ebp
Source line: 2
decl %eax
leal 1(%edi), %eax
popl %ebp
retl
nop
nop
nop
nop
nop
nop
@code_native add1()のようにしてもエラーが出てしまうのは
多重ディスパッチの影響があるので、add1(10)と具体的な値をつけてあげることで関数が適用する型を定めてあげるのだと思います。たとえばadd1(10.0)と浮動小数点数を入れたり、add1(BigInt(10))と多倍長整数にすると、生成されるコードも異なるものになります。
julia> @code_llvm add1(10.0)
define double @julia_add1_62439(double) #0 !dbg !5 {
top:
%1 = fadd double %0, 1.000000e+00
ret double %1
}
julia> @code_native add1(10.0)
.section __TEXT,__text,regular,pure_instructions
Filename: REPL[5]
pushl %ebp
decl %eax
movl %esp, %ebp
decl %eax
movl $200081256, %eax ## imm = 0xBECFF68
addl %eax, (%eax)
addb %al, (%eax)
Source line: 2
addsd (%eax), %xmm0
popl %ebp
retl
nop
nop
nop
nop
nop
nop
nop
nop
nop
nop
nop
nop
ちなみにCommon Lispにも言語仕様に似たものがありますね。（以下はSBCLの場合）
* (defun add1 (x)
(+ x 1))
ADD1
* (disassemble #'add1)
; disassembly for ADD1
; Size: 39 bytes. Origin: #x1001AC6434
; 34: 498B4C2458 MOV RCX, [R12+88] ; no-arg-parsing entry point
; thread.binding-stack-pointer
; 39: 48894DF8 MOV [RBP-8], RCX
; 3D: BF02000000 MOV EDI, 2
; 42: 488BD3 MOV RDX, RBX
; 45: 41BB8004B021 MOV R11D, #x21B00480 ; GENERIC-+
; 4B: 41FFD3 CALL R11
; 4E: 488B5DF0 MOV RBX, [RBP-16]
; 52: 488BE5 MOV RSP, RBP
; 55: F8 CLC
; 56: 5D POP RBP
; 57: C3 RET
; 58: 0F0B10 BREAK 16 ; Invalid argument count trap
NIL
Lauridsenのモノラル→ステレオのアップミックス・エフェクタをサクッと作ってみました。簡単ですねー。参考文献は以下。（注：最近使われている手法とは違ってかなり原始的です！）
desc: maruware Lauridsen Decorrelator (mono-to-stereo) // The ring buffer code was taken from // http://www.auriculaonline.com/wp/?p=232 slider1:20<1,100,1>Delay [ms] slider2:0<0,100,1>Mix [%] @init buf = 0; // buffer exists at offset 0 bufposR = 0; delay = srate * 0.010; bufposW = delay; buflength = srate * 0.100; mix = 0.0; @slider delay = srate * slider1 / 1000.0; // # of samples bufposR = bufposW - delay; bufposR < 0 ? bufposR = bufPosR + buflength; mix = slider2 / 100.0; @sample x = spl0; spl0 = x + buf[bufposR] * mix; spl1 = x - buf[bufposR] * mix; buf[bufposW] = x; bufposR = bufposR + 1 ; bufposR > buflength ? bufposR = 0; bufposW = bufposW + 1 ; bufposW > buflength ? bufposW = 0;
