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

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

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

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

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

パッケージのインストール

LibSndFileはオーディオファイルの入出力をするためのライブラリ、PortAudioはリアルタイムにオーディオの入出力をするためのライブラリで、〜.jlは両者をラップするパッケージです。LibSndFile.jlもPortAudio.jlもまだJulia 1.0に正式対応しているわけではないので、インストールから戸惑う部分があります。

Julia 0.7以降ではパッケージのインストールはjulia>プロンプトで]をタイプしてpkgモードに入って行います(プロンプトが(v1.0) pkg>のようになります)。そこで、add SIUnits#masteradd SampledSignals#masteradd LibSndFile#masterの三つを導入します。#masterは「開発版を使用するぞ」というオマジナイです。PortAudio.jlについては、9月10日の時点で#masterではうまくいかず、add PortAudio#julia1とするとうまく行きました。こちらは#masterよりもさらに初期のバージョンを使用することになるかと思うので今後すぐに状況が変わると思います。

LibSndFileの使い方

Julia 0.6まではusing LibSndFileで良かったような気がしますが、いまは以下のようにFileIOからloadsaveを使用し、そこにLibSndFileを重ねます。load(ファイル名)で、WAV、FLACOGGから読み込みができます。

julia> using FileIO: load, save

julia> import LibSndFile

julia> snd = load("recorder.wav")
352800-frame, 1-channel SampleBuf{FixedPointNumbers.Fixed{Int16,15}, 2}
8.0s sampled at 44100.0Hz
▁▅▅▅▅▅▅▅▅▅▅▅▅▅▅▅▅▂▁▁▁▁▁▁▁▁▁▁▁▁▅▅▅▅▅▆▆▆▆▆▆▆▆▆▆▆▆▅▁▁▁▁▁▁▁▁▁▁▁▅▅▆▆▆▆▆▆▆▆▆▆▆▆▆▆▆▅▁▁▁

snd.dataにオーディオのサンプルが、snd.samplerateに標本化周波数が入っています。

一方、xにオーディオデータがあるとして、それを保存するときには

using SampledSignals
buf = SampleBuf(x, samplerate);
save("output.wav", buf);

とします。

PortAudioの使い方

準備

PortAudioでは、まず使用するオーディオ入出力機器を確認します。

julia> using PortAudio

julia> PortAudio.devices()
7-element Array{PortAudio.PortAudioDevice,1}:
 PortAudio.PortAudioDevice("Built-in Microphone", "Core Audio", 2, 0, 48000.0, 0)
 PortAudio.PortAudioDevice("Built-in Output", "Core Audio", 0, 2, 48000.0, 1)
 PortAudio.PortAudioDevice("USB PnP Audio Device", "Core Audio", 0, 2, 48000.0, 2)
 PortAudio.PortAudioDevice("USB PnP Audio Device", "Core Audio", 1, 0, 48000.0, 3)
 PortAudio.PortAudioDevice("Soundflower (2ch)", "Core Audio", 2, 2, 48000.0, 4)
 PortAudio.PortAudioDevice("Soundflower (64ch)", "Core Audio", 64, 64, 48000.0, 5)
 PortAudio.PortAudioDevice("Fireface UC Mac (xxxxxxxx)", "Core Audio", 18, 18, 48000.0, 6)

バイスリストの更新はusing PortAudioをしたときだけに行われるようなので、使用したいオーディオインターフェースはあらかじめ接続しておきます。このリストは、デバイス名、ホストAPI、最大入力チャンネル数、最大出力チャンネル数、標本化周波数のデフォルト値、デバイス番号の順になっています。(Firefaceはxxxxxxxxの部分にシリアルナンバーが入っています)

バイスとのやりとりをするストリームを開くことで、オーディオ信号のやりとりが可能になります。たとえば1入力・2出力を確保したければ以下のようにします。

julia> stream = PortAudioStream(1, 2)
PortAudioStream{Float32}
  Samplerate: 48000.0Hz
  Buffer Size: 4096 frames
  2 channel sink: "Built-in Output"
  1 channel source: "Built-in Microphone"

ただ、入力・出力を別々に指定したり、チャンネル数を指定したり、標本化周波数やバッファサイズ場合を変更したりなども可能です。(ただ、ホストAPIから「このタイミングで標本化周波数は変えられないよ」などのエラーが出てくるかもしれません)

julia> stream = PortAudioStream("Built-in Microphone", "Soundflower (2ch)", 1, 2;
                                samplerate=44100, blocksize=512, synced=true)
PortAudioStream{Float32}
  Samplerate: 44100.0Hz
  Buffer Size: 512 frames
  2 channel sink: "Soundflower (2ch)"
  1 channel source: "Built-in Microphone"

あるいは、入出力両方に対応しているデバイス(full duplex)の場合は以下のようにデバイス名は1回指定でOKです。

julia> stream = PortAudioStream("Fireface UC Mac (xxxxxxxx)", 1, 2;
                                samplerate=48000, blocksize=512, synced=true)
PortAudioStream{Float32}
  Samplerate: 48000.0Hz
  Buffer Size: 512 frames
  2 channel sink: "Fireface UC Mac (xxxxxxxx)"
  1 channel source: "Fireface UC Mac (xxxxxxxx)"

ここで、Fireface UCのように入出力両方にチャンネル数が表示されているものは入出力の同期がとれますが、Built-inのMicとOutは別デバイスなので入出力の同期が取れないようです。僕の用途では入出力の遅延がわかっていることが大切なので、full duplexなデバイスを使用します。

オーディオの入出力

read()write()でストリームから/へSampledSignals.SampleBuf型で読み書きできます。readはストリームと、そのブロックサイズを指定します。writeはストリーム、データ、ブロックサイズの3引数です。以下ではマイク入力を512サンプル取り込み、無音(zeros())を512サンプル出力しています。

julia> buf_in = read(stream, stream.blocksize)
512-frame, 1-channel SampleBuf{Float32, 2}
0.010666666666666666s sampled at 48000.0Hz
▃▃▃▃▂▁▁▂▂▂▁▁▂▂▂▃▃▃▃▃▃▂▁▁▃▃▃▂▂▁▁▂▂▁▁▂▂▂▁▁▁▂▂▁▂▃▃▃▃▂▂▂▂▂▁▁▁▁▃▃▃▃▃▂▂▁▁▁▁▂▂▂▂▂

julia> write(stream, zeros(stream.blocksize), stream.blocksize)
512

full duplexなデバイスで入出力を同期させるときには、上記のように入力と出力のデータ量を対応させる必要があります。たとえばインパルス応答測定では、測定用信号をスピーカ出力して、それに対して部屋の応答をマイク入力から受けますので、測定信号をブロックサイズぶん出力してはマイク入力からブロックサイズぶん信号を受けるということを繰り返します。

具体的な入出力部分を見てみましょう。(このあたり、アヲギリさんのNotebookがかなり参考になったはずなのですが、具体的にどのNotebookだったのかは失念……)

blksize = stream.blocksize;

sig_out = load("input.wav");
sig_in = zeros(round(Int, floor(length(sig_out)/blksize + 5) * blksize));

buf_in = read(stream, blksize);
for i in 0:ceil(Int, length(sig_out)/blksize)-1
  # get input from microphone
  buf_in = read(stream, blksize);
  sig_in[i*blksize+1 : (i+1)*blksize] = convert(Array{Float32}, buf_in.data[:,1]);

  # output to loudspeaker
  if (i+1)*blksize > length(sig_out)
    buf_out = sig_out[i*blksize+1 : end]
  else
    buf_out = sig_out[i*blksize+1 : (i+1)*blksize];
  end
  write(stream, buf_out, blksize);
end

recorded = SampleBuf(sig_in[floor(Int, blksize*4+1) : floor(Int, blksize*4)+length(sig_out)], sig_out.samplerate);
save("output.wav", recorded);
close(stream);

sig_outには出力したい信号(インパルス応答測定であればスイープサイン音など)を入れておきます。sig_inは最終的に録音されファイル保存される信号です。for文でsig_outの長さに相当する回数だけループし、各ループにおいてbuf_outbuf_inはブロックサイズ(blksize)ごとに分割された細切れの信号を入出力しています。他に色々と面倒なことをやっているのは、sig_outの信号長がブロックサイズの整数倍になっていないときを想定したりしているからです。

特筆しておきたいのは、PortAudio.jlでは入出力のバッファとしてブロックサイズ2つぶんの長さのリングバッファを持っているようであることです。そのため、合計ブロックサイズ4つぶんの遅延が想定されます。そのため、上記コードでは、(1)sig_outの信号長がブロックサイズの整数倍未満のとき用に1ブロック、(2)内部のリングバッファ用に4ブロック、の計5ブロックを考慮した処理を入れています。

それらの遅延などが考慮して、最終的にrecordedにはsig_outと同じ長さの録音信号が入ります。WAVファイルに保存し、使用が終わったストリームはclose(stream)で閉じます。

まとめ

一応上記のコードでいくつか試しに測定してみましたが、まぁまぁうまくいっているようです。今回はオーディオ信号の出力→オーディオ信号の入力、というシンプルな部分について書きましたが、ほぼ同じやり方を使って入力に対して処理をして出力するエフェクタのようなものも書けます。

ここで紹介した内容はPortAudio#julia1版でのものなので、今後(わりとすぐに)変わる可能性があります。そうは言っても、Juliaでオーディオを扱う人が増えるといいなということで、参考までにPinkTSP信号を用いたインパルス応答測定のコード全体をおいておきますね。play_and_record()内のオーディオデバイス名を自身の環境に合わせて変える必要があります。

※これは現時点でのラフなコードですので、絶対に本番環境で使ってはいけません!!!

なお、PinkTSPについては以下の論文を参照して下さい。

  • 藤本卓也, “低域バンドでのSN比改善を目的としたTSP信号に関する検討” 日本音響学会研究発表会講演論文集, p.433. (1999)
  • 藤本卓也, “低域バンドでのSN比改善を目的としたTSP信号に関する検討 —高調波歪の除去—” 日本音響学会研究発表会講演論文集, p.555. (2000)
using FileIO: load, save
using DSP, SampledSignals, LibSndFile, PortAudio, FFTW
using Printf, Primes


"""
    duration = minimum length of the signal (seconds)
          fs = number of samples per second (Hz)
       nbits = number of bits per sample (bit)
        reps = number of repeated measurements (times)
"""
function pinktsp_generate(duration=5.0, samprate=48000, nbits=32, reps=4)
  N = 2 ^ ceil(Int, log2(duration * samprate));
  m = prevprime(N >> 2) + 1;
  
  a = 2 * m * pi / ( (N/2) * log(N/2) );
  
  H = complex(zeros(N));
  H[1] = 1;
  H[2:convert(Int, N/2+1)] = exp.(1im * a .* (1:N/2) .* log.(1:N/2)) ./ sqrt.(1:N/2);
  H[convert(Int, N/2+2):N] = conj(H[convert(Int, N/2):-1:2]);
  h = real(ifft(H));
  mi = argmin(abs.(h));
  h = [h[mi:end] ; h[1:mi-1]];
  h = reverse(h);
  
  hh = [repeat(h, outer=reps) ; zeros(length(h))];

  # write to an audio file
  hh = hh / maximum(abs.(hh)) * sqrt(1/2);   # peak at -3 dBFS
  if nbits==16
    hh2 = map(PCM16Sample, hh);
  elseif nbits==24
    hh2 = map(PCM24Sample, hh);
  elseif nbits==32
    hh2 = map(PCM32Sample, hh);
  else
    nbits = 32;
    hh2 = map(PCM32Sample, hh);
    error("unsupported bitrate (falls back to 32bits)");
  end
  buf = SampleBuf(hh2, samprate);
  return(buf);
end



## playback and record
function play_and_record(snd)
  fs = snd.samplerate;

  stream = PortAudioStream("Fireface UC Mac (xxxxxxxx)", 1, 1; samplerate=fs, blocksize=512, synced=true);
  blksize = stream.blocksize;
  
  data_in = zeros(round(Int, floor(length(snd)/blksize + 5) * blksize));
  buf_in = read(stream, blksize);
  for i in 0:ceil(Int, length(snd)/blksize)-1
    # get input from microphone
    buf_in = read(stream, blksize);
    data_in[i*blksize+1 : (i+1)*blksize] = convert(Array{Float32}, buf_in.data[:,1]);
    
    # output to loudspeaker
    if (i+1)*blksize > length(snd)
      buf_out = snd[i*blksize+1 : end]
    else
      buf_out = snd[i*blksize+1 : (i+1)*blksize];
    end
    write(stream, buf_out, blksize);
  end

  close(stream);
  
  buf = SampleBuf(data_in[floor(Int, blksize*4+1) : floor(Int, blksize*4)+length(snd)], fs);
  return(buf);
end


## deconvolve recorded signal with tsp signal
function pinktsp_analyze(tspsignal, recorded, reps)
  k = convert(Array{Float64}, tspsignal.data[:,1]);
  k = sum(reshape(k, (round(Int, length(k)/(reps+1)), reps+1)), dims=2);
  x = convert(Array{Float64}, recorded.data[:,1]);
  x = sum(reshape(x, (round(Int, length(x)/(reps+1)), reps+1)), dims=2);
  y = real(ifft(fft(x) ./ fft(k)));
  outir = SampleBuf(y, snd.samplerate);
  return(outir);
end



dur = 5.0;
nbits = 32;
reps = 4;
fs = 48000;
snd = pinktsp_generate(dur, fs, nbits, reps);
save("sweep.wav", snd);

buf = play_and_record(snd);
save("recorded.wav", buf);

outir = pinktsp_analyze(snd, buf, reps);
save("ir.wav", outir);