音のタイムストレッチ (Julia Advent Calendar 2024)

Julia Advent Calendar 2024の24日目の記事です。毎年Juliaで音に関係したことを書いていて、これまでにこのブログで書いたアドカレ記事は以下の通りです。

タイムストレッチとは

音のタイムストレッチ (time stretch) とは読んで字のごとく、音を時間方向に引き延ばすことです。録音した音をスロー再生すると再生時間が長くなりますよね。単純にスロー再生するだけだと音の高さも低くなってしまいますので、それを防いで音の高さはそのままにゆっくり再生したいときに使う技術です。タイムストレッチは音を扱うソフトの多くに実装されており、たとえばYouTubeでも動画の再生速度を変えても音の高さは変わらないようになっています。ただ、あまり高音質ではないものも散見されます。

今回は、ある程度の音質を確保しつつも簡単なアルゴリズムで実装する方法を紹介します。*1

まずは、このアルゴリズムで処理した例を聞いてみてください。

youtu.be

フーリエ変換

フーリエ変換を使うと、時間領域と周波数領域を行ったり来たりすることができます。音の場合は、波形をフーリエ変換すると周波数成分が得られます。1970年代にクーリーとテューキー*2が高速フーリエ変換のアルゴリズムを考案したため、コンピュータで手軽に実行できるようになり、こんにちでは当たり前のように使うことができるようになりました。

以下の図は高速フーリエ変換を使って波形を変換したときに、配列のどの部分にどんな情報が入るかを表しています。時間領域には1番目からn番目にかけて信号の振幅値が入っています。それを高速フーリエ変換を使って周波数領域に持っていくと、1番目には直流成分 (direct current) が入り、2番目からn/2番目までは各周波数ビンの振幅が入ります。n/2+1番目にはナイキスト周波数の成分が入ります。例えばサンプリング周波数が48000 HzでFFT点数が1024だとすると、ビン番号1〜513が0〜24000 Hzに対応することから各ビンの中心周波数を計算することができます。ビンの中に入っている値は複素数で、絶対値が振幅、偏角が位相を表します。

n/2+2番目からn番目までは負の周波数成分といって、正の周波数成分の複素共役の値が入っています。音響信号を扱うときには負の周波数成分は使わないことも多いので、FFTW.jlにはrfftとirfftという実数用の高速フーリエ変換関数が用意されています。以下では、それらを使ってタイムストレッチを実装していきます。

処理の概要

タイムストレッチは、入力波形を時間窓ごとに分割して処理を行なっていきます。フーリエ変換の前提として、波形は「無限に繰り返していること」「滑らかであること」「連続であること」があります。それを満たすために滑らかなハニング窓を使って処理をします。

いったん周波数領域に展開したあとでゼロ詰めを行います。FFTでは時間波形と周波数成分の配列が同じ大きさになるので、周波数領域でゼロ詰めをして時間波形に戻すと長い信号が得られます。ただし、各周波数ビンの位相成分をそのままにしておくと隣り合う時間窓で波形がうまくつながりません。そのため位相成分をランダムに並べ替える処理を入れてあります。

また、ステレオ信号に適用できるように、左と右信号を中央と側方の信号にしてから処理をかけるという工夫もしています。左右で位相の処理が変わると気持ち悪い広がり感が出てしまうので、位相の無作為化に使う乱数配列は同じものを使用します。

さらに、Distributed.jlのおかげで並列化も簡単でした。forループの外でSharedArrayを用意して、@syncと@distributedをつければ自動的にマルチコアを利用してくれます。

julia> @time include("time_stretch.jl")   #並列実行なし
 36.029736 seconds (1.18 M allocations: 81.857 GiB, 7.01% gc time, 0.79% compilation time)

julia> @time include("time_stretch.jl")   #並列実行あり
  9.220847 seconds (67.97 k allocations: 3.211 GiB, 2.48% gc time, 0.18% compilation time)

コード全体

using FileIO: load, save
import LibSndFile
using FFTW, DSP, SampledSignals
using Random
using Distributed, SharedArrays

# read sound file
snd = load("17.flac");   # Track 17 from EBU Tech 3287
fs = snd.samplerate;
x = snd.data;

# parameters
stretchratio = 15;   # how much you want to stretch
winduration = 0.5;   # window size in seconds
noverlap = 16;       # number of window overlap

winsize = round(Int, fs * winduration);
hopsize = winsize ÷ noverlap;
hann = hanning(winsize);
hannL = hanning(winsize * stretchratio);

# left-right to mid-side conversion
xM = zeros(ceil(Int, (ceil(size(x,1) ÷ winsize) + 1) * winsize));
xM[1:size(x,1)] .= convert.(Float64, x[:,1] .+ x[:,2]);
xS = zeros(ceil(Int, (ceil(size(x,1) ÷ winsize) + 1) * winsize));
xS[1:size(x,1)] .= convert.(Float64, x[:,1] .- x[:,2]);

# stretched signals go into these
xMy = SharedArray{Float64}(length(xM) * stretchratio);
xSy = SharedArray{Float64}(length(xS) * stretchratio);

@sync @distributed for n in 1:round(Int, length(xM) ÷ hopsize - (noverlap-1))
  # vector holding random permutation table
  perm = Random.randperm(stretchratio * winsize÷2 - 1) .+ 1

  # process time stretch on mid signal
  xMx = xM[(n-1)*hopsize+1 : (n-1)*hopsize+winsize];
  xMx = xMx .* hann;
  yM = FFTW.rfft([xMx ; zeros(winsize * (stretchratio-1))]);
  yMamp = abs.(yM);
  yMphs = angle.(yM);
  yMphs[2:end-1] = yMphs[perm];
  yMtmp = FFTW.irfft(yMamp .* exp.(-1.0im * yMphs), winsize * stretchratio);
  xMy[(n-1)*hopsize*stretchratio+1 : (n-1)*hopsize*stretchratio+winsize*stretchratio] += yMtmp .* hannL;

  # process time stretch on side signal
  xSx = xS[(n-1)*hopsize+1 : (n-1)*hopsize+winsize];
  xSx = xSx .* hann;
  yS = FFTW.rfft([xSx ; zeros(winsize * (stretchratio-1))]);
  ySamp = abs.(yS);
  ySphs = angle.(yS);
  ySphs[2:end-1] = ySphs[perm];
  yStmp = FFTW.irfft(ySamp .* exp.(-1.0im * ySphs), winsize * stretchratio);
  xSy[(n-1)*hopsize*stretchratio+1 : (n-1)*hopsize*stretchratio+winsize*stretchratio] += yStmp .* hannL;
end

# convert mid-side back to left-right
yL = (xMy .+ xSy) ./ 2;
yR = (xMy .- xSy) ./ 2;

# make it stereo, normalize, and save
y = [yL yR];
y = y ./ maximum(abs.(y)) / 256 * 255;
y = y[1:size(x,1)*stretchratio,:];
sndout = SampleBuf(convert.(PCM16Sample, y), fs);
save("output.wav", sndout);