フェーズボコーダーによるタイムストレッチ

はじめに

フェーズボコーダーについての備忘録.
数式も英語も碌に読めない素人の独学なので内容の正確性は保証しない.
例として挙げるコードは C# で書く.
また、コード中で複素数を表す型 Complex は System.Numerics.Complex である.
コード中に出てくる Fourier.STFT、Fourier.ISTFT は私の前の記事 (短時間フーリエ変換とその逆変換) で書いたものである.

spectrum に対応する日本語としては スペクトラム を使用する.
一般的にはフランス語の spectre 由来で スペクトル と呼ばれることが多いのだが、コード中は spectrum であるし、参考文献も読んで精々英語であるためまあいいかなって.

タイムストレッチとは

音声を単純に2倍速にする(=音声の長さを1/2倍にする)と、音高はオクターブ上がる(=周波数が2倍になる).
逆に1/2倍速にする(=音声の長さを2倍にする)と、音高はオクターブ下がる(=周波数が1/2倍になる).

タイムストレッチは、音高が変わらないようにしつつ音声の長さを変更する処理である.
最近 YouTube をはじめとする動画サイトなどで再生速度を変更しても音高が変わらないのはこのおかげ(YouTube で使用しているのはフェーズボコーダーではないらしいが).

この記事ではタイムストレッチを行うための技法の一種であるフェーズボコーダーについて見ていく.

フェーズボコーダーの気持ち

フェーズボコーダーですることを簡単にまとめると次のようになる:

1. ずらし量 Ha で STFT を行う. この部分の処理を指して分析<analysis> とも.
2. スペクトログラムを編集する. この部分の処理を指して編集<modify> や処理<process> とも.
3. ずらし量 Hs で iSTFT を行う. この部分の処理を指して(再-)合成<(re-)synthesis> とも.

このとき、Hs が Ha の S 倍であれば、出力信号の長さも入力信号の (だいたい) S 倍になるという算段である.

スペクトログラムの編集の意義

上の手順で2番目の、編集処理が必要な理由について考える.

k 番目、k+1 番目のフレームについて、ずらし量 Ha の場合と Hs の場合を図示する(スペクトラムは見づらいので波形で示す):

k 番目のフレーム(前側)とk+1 番目のフレーム(後ろ側)について、重なっている部分は同じ波形になっていることが望ましい.
入力波形をずらし量 Ha で分割したスペクトログラムについてずらし量を Hs に変更しただけでは当然そうならないため、この辻褄合わせに編集が必要なのである.

スペクトログラムによって各周波数成分の大きさと(初期)位相が分かるわけだが、位相の方だけ編集すればおおよそ問題ないため、フェーズボコーダーの中心は位相 (phase) の編集となる.
それゆえ、PHASE vocoder の名前がついている.

このことから、フェーズボコーダーのテンプレートとなるコードは以下のように書ける:

コード例

/// <summary>
/// フェーズボコーダーによるタイムストレッチの実装
/// </summary>
/// <param name="input">波形データ</param>
/// <param name="frameSize">フーリエ変換を行う窓の長さ</param>
/// <param name="analysisHop">分析時のフレームずらし量(Ha)</param>
/// <param name="synthesisHop">合成時のフレームずらし量(Hs)</param>
/// <param name="window">窓関数の値(長さが frameSize である配列)</param>
/// <returns>タイムストレッチされた波形データ</returns>
public static Complex[] TimeStretch(Complex[] input, int frameSize, int analysisHop, int synthesisHop, double[] window) {
  // ずらし量 analysisHop で STFT を行う
  Complex[][] spectrogram = Fourier.STFT(input, frameSize, analysisHop, window);

  // スペクトログラムを編集する
  Complex[][] newSpectrogram = PhaseVocoderProcess(spectrogram, analysisHop, synthesisHop);

  // ずらし量 synthesisHop で STFT を行う
  Complex[] output = Fourier.ISTFT(newSpectrogram, frameSize, synthesisHop, window);
  return output;
}// TimeStretch(Complex[] input, int frameSize, int analysisHop, int synthesisHop, double[] window)

private static Complex[][] PhaseVocoderProcess(Complex[][] spectrogram, int analysisHop, int synthesisHop) {
  int numFrames = spectrogram.Length;
  int frameSize = spectrogram[0].Length;

  // 計算のため、複素スペクトログラムを振幅スペクトログラムと位相スペクトログラムに分ける
  (double[][] magnitudeSpectrogram, double[][] phaseSpectrogram) = AnalyseComplexSpectrogram(spectrogram);

  // 新しい位相スペクトログラムを求める
  double[][] newPhaseSpectrogram = /* TODO */;

  // 振幅スペクトログラムと位相スペクトログラムから複素スペクトログラムを構築する
  // 振幅スペクトログラムにはノータッチ
  Complex[][] newSpectrogram = ReconstructSpectrogram(magnitudeSpectrogram, newPhaseSpectrogram);
  return newSpectrogram;
}

/** vvvvvvvvvvvvvvvvvvvvvvvvvvvv ** 以下、補助関数 ** vvvvvvvvvvvvvvvvvvvvvvvvvvvv **/

private static (double[][] magnitudeSpectrogram, double[][] phaseSpectrogram) AnalyseComplexSpectrogram(Complex[][] spectrogram) {
  int numFrames = spectrogram.Length;
  int frameSize = spectrogram[0].Length;

  double[][] magnitudeSpectrogram = new double[numFrames][];
  double[][] phaseSpectrogram = new double[numFrames][];

  for (int i = 0; i < numFrames; ++i) {
    magnitudeSpectrogram[i] = new double[frameSize];
    phaseSpectrogram[i] = new double[frameSize];
    for (int j = 0; j < frameSize; ++j) {
      magnitudeSpectrogram[i][j] = spectrogram[i][j].Magnitude;
      phaseSpectrogram[i][j] = spectrogram[i][j].Phase;
    }
  }
  return (magnitudeSpectrogram, phaseSpectrogram);
}

private static Complex[][] ReconstructSpectrogram(double[][] magnitudeSpectrogram, double[][] phaseSpectrogram) {
  int numFrames = magnitudeSpectrogram.Length;
  int frameSize = magnitudeSpectrogram[0].Length;

  Complex[][] spectrogram = new Complex[numFrames][];
  for (int i = 0; i < numFrames; ++i) {
    spectrogram[i] = new Complex[frameSize];
    for (int j = 0; j < frameSize; ++j) {
      // System.Numerics.Complex には関数があるので利用しているが、
      // 大きさ r、偏角 θ の極座標表示された複素数を直交座標にする場合: r * exp(i θ)
      spectrogram[i][j] = Complex.FromPolarCoordinates(magnitudeSpectrogram[i][j], phaseSpectrogram[i][j]);
    }
  }
  return spectrogram;
}

フェーズボコーダー

では、どのようにして k+1 番目のフレームの b 番目の周波数成分の位相を求めればよいかが課題となる.
単純なフェーズボコーダーでは、これは次の手順で行われる.

1. k+1 番目のフレームの b 番目のビンについて瞬間の角周波数を求める.
   1. b 番目の周波数ビンについて、中心となる角周波数を求める (これを ω とする).
   2. k 番目のフレームから k+1 番目のフレームまでに、b 番目の位相がどれだけ進んだかを求める (これを φ_a とする[phase actual]).
   3. 角周波数 ω で幅 Ha 進んだ場合の角度を求める (これを φ_e とする[phase expected]).
   4. ω + (φ_a - φ_e) / Ha が求めるべき角周波数である.
2. k 番目のフレームの初期位相から幅 Hs だけ進んだ位相を求める.
   1. その差は、上で求めた角周波数に Hs を乗ずればよい.

以上から、新しい位相スペクトログラムを求める処理は次のように書ける:

コード例

private static Complex[][] PhaseVocoderProcess(Complex[][] spectrogram, int analysisHop, int synthesisHop) {
  /*
   * 前略
   */

  // 各周波数ビンに対応する周波数の中心
  double[] omega = new double[frameSize];
  for (int i = 0; i < frameSize; ++i) {
    omega[i] = WrapToPi(2 * Math.PI * i / frameSize);
  }
  // フレームひとつの間に実際に進んだ位相量と omega の差(φ_a - φ_e)
  double[][] deltaPhi = new double[numFrames][];
  deltaPhi[0] = new double[frameSize];
  for (int j = 0; j < frameSize; ++j) {
    deltaPhi[0][j] = 0;
  }
  for (int i = 1; i < numFrames; ++i) {
    deltaPhi[i] = new double[frameSize];
    for (int j = 0; j < frameSize; ++j) {
      deltaPhi[i][j] = WrapToPi(phaseSpectrogram[i][j] - phaseSpectrogram[i - 1][j]) - analysisHop * omega[j];
    }
  }
  // 各フレーム、各ビンの瞬間角周波数
  double[][] instantaneousFrequencies = new double[numFrames][];
  for (int i = 0; i < numFrames; ++i) {
    instantaneousFrequencies[i] = new double[frameSize];
    for (int j = 0; j < frameSize; ++j) {
      instantaneousFrequencies[i][j] = omega[j] + deltaPhi[i][j] / analysisHop;
    }
  }

  // 角周波数から位相スペクトログラムを組み立てる
  double[][] newPhaseSpectrogram = ReconstructPhaseSpectrogram(phaseSpectrogram[0], instantaneousFrequencies, synthesisHop);

  /*
   * 後略
   */
}

/** vvvvvvvvvvvvvvvvvvvvvvvvvvvv ** 以下、補助関数 ** vvvvvvvvvvvvvvvvvvvvvvvvvvvv **/

/// <summary>
/// 引数を [-PI .. PI] の範囲にラッピングする
/// </summary>
/// <param name="phase">角度</param>
/// <returns>[-PI .. PI] の範囲にラッピングされた角度</returns>
private static double WrapToPi(double phase) {
  const double m = 2 * Math.PI;
  double q = Math.Round(phase / m);
  return phase - q * m;
}

/// <summary>
/// 与えられた角周波数から位相スペクトログラムを求める
/// </summary>
/// <param name="initialPhaseSpectrum">初期位相</param>
/// <param name="instantaneousFrequencies">各成分の角周波数</param>
/// <param name="hopSize">フレームずらし量</param>
/// <returns>位相スペクトログラム</returns>
private static double[][] ReconstructPhaseSpectrogram(double[] initialPhaseSpectrum, double[][] instantaneousFrequencies, int hopSize) {
  int numFrames = instantaneousFrequencies.Length;
  int frameSize = instantaneousFrequencies[0].Length;

  double[][] newPhaseSpectrogram = new double[numFrames][];
  newPhaseSpectrogram[0] = initialPhaseSpectrum;
  for (int i = 1; i < numFrames; ++i) {
    newPhaseSpectrogram[i] = new double[frameSize];
    for (int j = 0; j < frameSize; ++j) {
      newPhaseSpectrogram[i][j] = newPhaseSpectrogram[i - 1][j] + instantaneousFrequencies[i][j] * hopSize;
    }
  }
  return newPhaseSpectrogram;
}

改良されたフェーズボコーダー

上述のフェーズボコーダーでは、k 番目のフレームと k + 1 番目のフレームにおける b 番目の周波数成分の関係、つまり時間軸における周波数成分の整合性だけを考慮していた.
この時間軸における成分の整合性は水平コヒーレンスと呼ばれ、たしかに重大な意味を持つのだが、スペクトログラムは二次元データであるから、もうひとつ考えるべき関係性が存在する.
その関係性、周波数領域における周波数成分の整合性 (垂直コヒーレンス) はたしかに存在し、そちらについても考慮することによってフェーズボコーダーの性能を高めるべく研究が行われている.

スペクトラム漏れ

周波数成分間の整合性を考慮すべきケースとして、周波数成分の漏れを見ていく.

サンプルレート=256、周波数=16 Hz の余弦波と、それを DFT した結果を下に示す:

スペクトラムのピーク付近の値:

bin	Magnitude( / 256)	Phase( / 2 π)
11	0.000000	0.314736
12	0.000000	0.335939
13	0.000000	0.492690
14	0.000000	0.440294
15	0.000000	0.460864
16	0.500000	0.000000
17	0.000000	0.026738
18	0.000000	0.068382
19	0.000000	0.029771
20	0.000000	0.122023
21	0.000000	-0.035124

整数番目のビンに丁度収まる (=サンプルレートを割り切れる周波数の成分) はそのビンの成分だけが大きく、近隣のビンの成分はほとんど 0 になる.

対して、サンプルレート=256、周波数=16.1 Hz の余弦波の場合は以下のようになる:

スペクトラムのピーク付近の値:

bin	Magnitude( / 256)	Phase( / 2 π)
11	0.011456	0.034103
12	0.013744	0.037214
13	0.017553	0.040328
14	0.025052	0.043444
15	0.046290	0.046564
16	0.493346	0.049687
17	0.053162	-0.447186
18	0.024445	-0.444055
19	0.015560	-0.440919
20	0.011250	-0.437780
21	0.008713	-0.434636

この例では 256 [サンプル] / 16.1 [Hz] = 15.900... となり、収まるべきビンが整数番目でない.
この場合 DFT ではこの成分が近隣のビンに漏れ出てしまう.
こうして、本来その波形に含まれない周波数成分が出てくることをスペクトラム漏れ (spectral leakage) とか smearing (にじみ、汚れ、の意) などと呼ぶ.

このため、フェーズボコーダー(無印)で考慮していた水平コヒーレンスだけではなく、垂直コヒーレンスについても考慮すべきなのである.

Loose phase-locking vocoder

隣り合う周波数成分の位相 (phase) をある程度同期させる (locking) ことによって垂直コヒーレンスに対応するものである.
元の論文 (Puckette, 1995) では単に phase-locked vocoder だが、後に "Rigid" phase-locking vocoder が登場したためか、こちらは "loose" と付いている.

フェーズボコーダー(無印)に基づいて位相を求めた後各周波数ビンについて、その隣の成分が大きければそれの影響を受けた位相に修正する、という具合だ.

コード例

private static Complex[][] PhaseVocoderProcess(Complex[][] spectrogram, int analysisHop, int synthesisHop) {
  /*
   * 前略
   */

  double[][] tmpPhaseSpectrogram = /* フェーズボコーダー(無印) における newPhaseSpectrogram */

  double[][] newPhaseSpectrogram = new double[numFrames][];
  newPhaseSpectrogram[0] = tmpPhaseSpectrogram[0];
  for (int i = 1; i < numFrames; ++i) {
    newPhaseSpectrogram[i] = new double[frameSize];
    for (int j = 0; j < frameSize; ++j) {
      // 0 番目のビンにとっての「前」、最後のビンにとっての「後ろ」のビンをどうすべきかはよくわからない
      // とりあえず巡回するものとして扱う.
      int prevIdx = j - 1 >= 0? j - 1 : frameSize - 1;
      int currIdx = j;
      int nextIdx = j + 1 < frameSize? j + 1 : 0;

      // 現在のビンと、その前後のビンについて複素数としての周波数成分を得る
      Complex prev = Complex.FromPolarCoordinates(magnitudeSpectrogram[i][prevIdx], tmpPhaseSpectrogram[i][prevIdx]);
      Complex curr = Complex.FromPolarCoordinates(magnitudeSpectrogram[i][currIdx], tmpPhaseSpectrogram[i][currIdx]);
      Complex next = Complex.FromPolarCoordinates(magnitudeSpectrogram[i][nextIdx], tmpPhaseSpectrogram[i][nextIdx]);

      // 現在 - (前 + 後ろ) の位相成分を採用する
      newPhaseSpectrogram[i][j] = (curr - prev - next).Phase;
    }// loop for j, in [0 .. frameSize)
  }// loop for i, in [0 .. numFrames)

  /*
   * 後略
   */
}

Rigid phase-locking vocoder

Loose phase-locking vocoder ではすぐ隣の成分に影響を受ける/与えるという程度だった.
それに対して、Rigid phase-locking vocoder はスペクトラム中の山の頂点となるビンが周囲に影響を及ぼすという形を取る.

領域の分割

以下のような領域分割を行うことを考える:

ここで、赤丸はピーク、ピンクと青がそれぞれのピークが支配する領域である.
それぞれの領域にはただひとつのピークが存在している.
(領域を示すのにピンクと青を交互に並べたが、単に四角で囲むだけの方が見やすいか？)

このピーク検出、およびそれが支配的である領域の取得を行うコードの例を示す:

コード例

/// <summary>
/// データを領域に分割する.
/// すべてのピークは、それが影響を与える領域をただひとつ持つ.
/// この実装では各領域の境界をピーク同士の間で最小の値をもつデータとしている.
/// </summary>
/// <param name="data">入力データ</param>
/// <param name="threshold">閾値</param>
/// <returns>{ 領域のリスト, その領域に対して影響を持つピークのインデックスのリスト }</returns>
public static (List<(int leftInclusive, int rightExclusive)> regions, List<int> peakIndices) DivideIntoRegions(double[] data, double threshold) {

  // ピーク検出を行う
  List<int> peakIndices = DetectPeaks(data, threshold);

  if (data.Max() < threshold || peakIndices.Count == 0) {
    // 最大値が最小のピーク値以下である場合、またはピークが見つからない場合、全体でひとつの領域として扱う
    List<(int leftInclusive, int rightExclusive)> regions = new List<(int leftInclusive, int rightExclusive)> { (0, data.Length) };
    peakIndices = new List<int> { Array.IndexOf(data, data.Max()) };
    return (regions, peakIndices);
  }

  // データを複数の領域に分割する
  var detectedRegions = new List<(int leftInclusive, int rightExclusive)>();

  int borderLeft = 0;
  for (int j = 0; j < peakIndices.Count; ++j) {
    int borderRight;
    if (j == peakIndices.Count - 1) {
      // 右端のピークの場合、この領域に最後のサンプルまで含める.
      borderRight = data.Length;
    } else {
      var left = peakIndices[j];
      var right = peakIndices[j + 1];
      var slice = Enumerable.Range(left, right - left + 1).ToArray();
      // data[left .. right) 内の最小値 min を求める.
      // これが、この領域と右隣の境界となる.
      var min = slice.Min(k => data[k]);
      var idx = Array.IndexOf(data, min, left, right - left + 1);
      borderRight = idx;

      // "Improved Phase Vocoder Time-Scale Modication of Audio" によれば、
      // 単純に次のような実装でも問題はないらしい:
      // borderRight = (left + right) / 2;
    }

    // 領域は、data[left .. right) であるように定める
    (int borderLeft, int borderRight) region = (borderLeft, borderRight);
    detectedRegions.Add(region);
    borderLeft = borderRight;
  }// loop for j; indices.Length

  return (detectedRegions, peakIndices);
}// GetRegions(double[] data, double threshold)

/// <summary>
/// 閾値以上である点の内ピークであるものを見つける.
/// cf. 参考文献 "Peak Detection in a Measured Signal"
/// 参考文献 "Improved Phase Vocoder Time-Scale Modication of Audio" では
/// 最も単純な実装として「近隣の4つのビンよりも大きな値を持つビンをピークとする」方法を挙げているが、
/// せっかくなのでそれらしいピーク検出関数を用意.
/// けれど、こちらは閾値を要求するので近隣のビンと比較するだけの方がいいかも.
/// </summary>
/// <param name="data">入力データ</param>
/// <param name="threshold">閾値</param>
/// <returns>ピークのインデックスのリスト. 存在しない場合空のリスト.</returns>
public static List<int> DetectPeaks(double[] data, double threshold) {
  List<int> peakIndices = new List<int>();
  int peakIndex = -1;
  double peakValue = double.NegativeInfinity;
  for (int i = 0; i < data.Length; ++i) {
    if (data[i] > threshold) {
      if (peakValue < data[i]) {
        peakIndex = i;
        peakValue = data[i];
      }
    } else if (data[i] < threshold && peakIndex != -1) {
      peakIndices.Add(peakIndex);
      peakIndex = -1;
      peakValue = double.NegativeInfinity;
    }
  }// loop for i, in [0 .. data.Length)
  if (peakIndex != -1) {
    peakIndices.Add(peakIndex);
  }
  return peakIndices;
}// DetectPeaks(double[] data)

Identity phase-locking vocoder

ピークについてはフェーズボコーダー(無印)同様に計算し、その後そのピークが影響力を持つビンについて計算する.
この影響による計算は以下のコード例中の synthesisPhasePeak + phaseSpectrogram[i][k] - phaseSpectrogram[i][peak]; 部分.
計算されたピークの位相から、そのビンとピークの分析時の位相の差を加える形になっている.

コード例

private static Complex[][] PhaseVocoderProcess(Complex[][] spectrogram, int analysisHop, int synthesisHop) {
  /*
   * 前略
   */

  // それぞれのスペクトラムを領域に分割する
  (List<(int leftInclusive, int rightExclusive)> regions, List<int> peakIndices)[] regionsOnSpectrogram = new (List<(int leftInclusive, int rightExclusive)> regions, List<int> peakIndices)[numFrames];
  for (int i = 0; i < numFrames; ++i) {
    double[] spectrum = magnitudeSpectrogram[i];
    double threshold = spectrum.Max() / 5;// ピーク検出用の閾値. 何を使うべき？
    regionsOnSpectrogram[i] = PeakFinder.DivideIntoRegions(spectrum, threshold);
  }

  // 各周波数ビンに対応する周波数の中心
  double[] omega = new double[frameSize];
  for (int i = 0; i < frameSize; ++i) {
    omega[i] = WrapToPi(2 * Math.PI * i / frameSize);
  }

  double[][] newPhaseSpectrogram = new double[numFrames][];
  newPhaseSpectrogram[0] = phaseSpectrogram[0];

  for (int i = 1; i < numFrames; ++i) {
    newPhaseSpectrogram[i] = new double[frameSize];
    (List<(int leftInclusive, int rightExclusive)> regions, List<int> peakIndices) = regionsOnSpectrogram[i];

    // 領域単位で位相を求める
    for (int j = 0; j < peakIndices.Count; ++j) {
      (int leftInclusive, int rightExclusive) region = regions[j];
      int peak = peakIndices[j];
      // phaseSpectrogram[i][peak] に対して newPhaseSpectrogram[i][peak] を求める.
      // cf. フェーズボコーダー(無印)のとき
      // deltaPhi[i][j] = WrapToPi(phaseSpectrogram[i][j] - phaseSpectrogram[i - 1][j]) - analysisHop * omega[j];
      // instantaneousFrequencies[i][j] = omega[j] + deltaPhi[i][j] / analysisHop;
      // newPhaseSpectrogram[i][j] = newPhaseSpectrogram[i - 1][j] + instantaneousFrequencies[i][j] * synthesisHop;
      double deltaPhiPeak = WrapToPi(phaseSpectrogram[i][peak] - phaseSpectrogram[i - 1][peak] - analysisHop * omega[peak]);
      double instantaneousFrequencyPeak = omega[peak] + deltaPhiPeak / analysisHop;
      double synthesisPhasePeak = newPhaseSpectrogram[i - 1][peak] + instantaneousFrequencyPeak * synthesisHop;
      // ピークの影響を周囲に割り振る
      for (int k = region.leftInclusive; k < region.rightExclusive; ++k) {
        newPhaseSpectrogram[i][k] = synthesisPhasePeak + phaseSpectrogram[i][k] - phaseSpectrogram[i][peak];
      }// loop for k; in region
    }// loop for j; peakIndices.Count
  }// loop for i; numFrames

  /*
   * 後略
   */
}// PhaseVocoderProcess_IdentityPhaseLocking(Complex[][] spectrogram, int analysisHop, int synthesisHop)

Scaled phase-locking vocoder

Identity phase-locking とほとんど同じだが、係数がひとつ増えている.

コード例

private static Complex[][] PhaseVocoderProcess(Complex[][] spectrogram, int analysisHop, int synthesisHop) {
  // 参考文献 "Improved Phase Vocoder Time-Scale Modication of Audio" によると
  // β ~ 2/3+α/3　程度がいいらしい？
  double alpha = (double) synthesisHop / analysisHop;
  double beta = (alpha + 2) / 3;

  /*
   * 前略
   */

  for (int i = 1; i < numFrames; ++i) {
    : : : 中略
    for (int j = 0; j < peakIndices.Count; ++j) {
      : : : 中略
      for (int k = region.leftInclusive; k < region.rightExclusive; ++k) {
        // Identity phase-locking との実質唯一の差
        newPhaseSpectrogram[i][k] = synthesisPhasePeak + beta * (phaseSpectrogram[i][k] - phaseSpectrogram[i][peak]);
      }// loop for k; in region
    }// loop for j; peakIndices.Count
  }// loop for i; numFrames

  /*
   * 後略
   */
}

参考文献

• Stephan Bernsee's Blog "Pitch Shifting Using The Fourier Transform" http://blogs.zynaptiq.com/bernsee/pitch-shifting-using-the-ft/
  smearing と瞬間の周波数の計算について参考にした
• Miller Puckette, 1995, "Phase-locked Vocoder"
  loose phase-locking の論文. フェーズボコーダー(無印)と loose phase-locking について.
• Florian Hammer, 2001, "Time-scale Modification using the Phase Vocoder"
  各フェーズボコーダーについて
• Johannes Grünwald, 2010, "Theory, Implementation and Evaluation of the Digital Phase Vocoder in the Context of Audio Effects"
  各フェーズボコーダーについて
• Jean Laroche and Mark Dolson, 1999, "Improved Phase Vocoder Time-Scale Modication of Audio"
  scaled phase-locking の論文. 各フェーズボコーダーについて.
• Baeldung "Peak Detection in a Measured Signal" https://www.baeldung.com/cs/signal-peak-detection
  rigid phase-locking で使用するピーク検出について