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

この記事は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

Julia 1.0でオーディオ入出力(とインパルス応答測定)

Juliaが1.0になって、徐々に環境が整ってきました。ここ数日、LibSndFile.jlPortAudio.jlを使ってゴニョゴニョやっていたことをまとめておきます。

TLDR。入出力の同期とか考えずに最低限の音再生ができればいいという場合は、

using FileIO: load, save
import LibSndFile
using PortAudio
snd = load("hoge.wav");
stream = PortAudioStream(0, 2);
write(stream, snd);

とすれば、ファイルを読み込んで、デフォルトの2チャンネル再生デバイス(内蔵スピーカとか)から再生してくれます。

続きを読む

Juliaコードのディスアセンブル(Julia 1.0)

marui.hatenablog.com

上記の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

Juliaコードのディスアセンブル(Julia 0.6)

www.youtube.com

ヨーロッパ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

REAPERのJSFXでモノラル→ステレオのアップミックスやってみた

Lauridsenのモノラル→ステレオのアップミックス・エフェクタをサクッと作ってみました。簡単ですねー。参考文献は以下。(注:最近使われている手法とは違ってかなり原始的です!)

www.aes.org

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;