ウェーブテーブルのピッチベンド

ピッチベンドをかけても倍音の欠落やエイリアシングノイズが出ないウェーブテーブルを実装します。

ここで紹介するテクニックは Laurent de Soras さんによる The Quest For The Perfect Resampler の内容を基にしています。

エイリアシング低減の仕組み

ウェーブテーブルのピッチベンドの挙動を見ていきます。例としてサンプリング周波数 \(f_s = 48000\,\text{[Hz]}\) 、基本周波数 \(f_0 = 880\,\text{[Hz]}\) として帯域制限したウェーブテーブルを生成します。帯域制限の方法についてはウェーブテーブルの帯域制限と位相方向の補間を参照してください。

生成したウェーブテーブルを指定したピッチで再生します。

以下はウェーブテーブルを一定の周波数 \(f_0\) で再生したときのスペクトログラムです。明るい線がウェーブテーブルの倍音です。再生時の周波数が一定なので横一直線に伸びています。

以下は基本周波数より高くなるようにウェーブテーブルを \(f_0\) から \(2 f_0\) までピッチベンドしたときのスペクトログラムです。明るい線が上端で跳ね返っているように見えるのがエイリアシングノイズです。

以下は基本周波数より低くなるようにウェーブテーブルを \(f_0\) から \(0.5f_0\) までピッチベンドしたときのスペクトログラムです。高次の倍音が欠けています。

倍音の欠落とエイリアシングノイズの問題を解決するためにはオーバーサンプリングが使えます。オーバーサンプリングの倍率を \(L\) とするとナイキスト周波数を超えない最大の倍音の次数は \(\displaystyle N_h = \left\lfloor \frac{f_s}{2 L f_0} \right\rfloor\) です。例として 2 倍のオーバーサンプリングを行います。

以下は基本周波数より高くなるようにウェーブテーブルを \(f_0\) から \(2 f_0\) までピッチベンドしたときのスペクトログラムです。オーバーサンプリングによってできたマージンがあるのでエイリアシングが出ていません。図の左の周波数軸の値が 2 倍になっている点に注意してください。

マージンに現れた成分はダウンサンプリング時にかけるローパスフィルタで低減できます。以下はローパスフィルタをかけた後の信号のスペクトログラムです。

以下はフィルタをかけた信号をデシメーションした後のスペクトログラムです。オーバーサンプリングなしで高いほうに向かってピッチベンドしたときと比べると、エイリアシングノイズが消えていることが確認できます。また、ダウンサンプリングのフィルタによって 20000 Hz 以上の成分が低減されているので、図の上部が真っ黒になっています。

アップサンプリングの倍率を \(L\) とするとエイリアシングノイズが出ないピッチベンドの範囲は \([1, 2L - 1)\) です。あとは出力したい周波数の範囲に応じて複数のウェーブテーブルを用意して、入力されたピッチによって再生するウェーブテーブルを切り替えることで、ピッチベンドしてもエイリアシングノイズが出ないオシレータが作れます。

基本周波数 \(f_0\) より低くなるようにピッチベンドするとオーバーサンプリングの有無にかかわらず倍音の欠落が防げません。よって切り替えを行うときには、入力された周波数 \(f\) 以下のウェーブテーブルの中から、基本周波数 \(f_0\) が \(f\) に最も近いものを選びます。

以下は 1 オクターブ間隔で生成したウェーブテーブルと基本周波数、インデックスの対応表です。インデックスがそのまま基本周波数のオクターブ差になっています。表の半音 (semitone) は \(f_L\) を 0 とした相対的な値です。

以下のリンクは、ここまで掲載した図の作成に使ったウェーブテーブルの生成とピッチベンドのテストコードです。 Python 3 です。

• 1 つのウェーブテーブルのピッチベンドのテストコードを読む (github.com)

設計

ここでは以下のパラメータから他の必要な値を決めます。

• 周波数が \(f_L\) 以下になると倍音が欠落する。
• 1 つのウェーブテーブルを \(B\) 倍までピッチベンドする。
• 1 つのウェーブテーブルはナイキスト周波数の \(H\) 倍までの周波数成分を含む。
• \(L\) 倍のオーバーサンプリング。

このとき以下の不等式を満たせばエイリアシングは起こりません。

\[ BH \leq 2L - 1 \]

ピッチを上げる方向に向かってだけピッチベンドするので \(B \geq 1\) 、倍音の欠落が起こらないようにするので \(H \geq 1\) です。

\(2L - 1\) はオーバーサンプリングによってできる周波数方向のマージンを表しています。 \(B, H, L\) はすべて比率なのでサンプリング周波数は式に入っていません。

\(B\) と \(H\) の決め方は実装の節で実験しながら見ていきます。

ウェーブテーブルの作成

まずは倍音が欠けない最低の周波数 \(f_L\) を決めます。そして \(f_L\) を使ってウェーブテーブルの長さ \(N\) を求めます。インデックスの巻き戻しでビットマスクを使うために、ウェーブテーブルの長さは \(2^n\) に揃えます。

\[ N = 2^n, \quad n = \left\lfloor \log_2 \frac{f_s}{f_L} \right\rfloor \]

次にウェーブテーブルの数 \(M\) を求めます。 1 つのウェーブテーブルのピッチベンドの倍率を \(B\) とすると以下の式が立ちます。

\[ B^M = \frac{N}{2} \]

\(\dfrac{N}{2}\) は 1 つのウェーブテーブルが含むことのできる最大の倍音の数です。上の式は最低周波数から \(B\) 倍のピッチベンドを \(M\) 回繰り返したときにナイキスト周波数以下の倍音が 1 つだけになることを表しています。

\(M\) について解きます。 \(M\) を整数にしたいので床関数を加えています。

\[ M = \left\lfloor \log_B \frac{N}{2} \right\rfloor \]

最後にウェーブテーブルを作ります。ウェーブテーブルの作成は入力された波形を離散フーリエ変換で周波数領域に変換します。そしてエイリアシングが出ないように周波数成分を切り落としたあとに逆離散フーリエ変換で時間領域に戻すことで計算できます。ここではインデックス 0 のウェーブテーブルが最低周波数を再生するようにします。ウェーブテーブルのインデックスを \(i\) としたときに切り落とす倍音の次数 \(c\) は以下の式で計算できます。

\[ c = \left\lfloor \frac{N}{2} H^{-i} \right\rfloor + 1 \]

最後の \(+1\) は直流成分を表しています。インデックスが \(c\) 以上の周波数成分は 0 にします。

以下はウェーブテーブル作成の大まかな C++ のコードです。

/*
table は 2 次元配列。

sampleRate は出力のサンプリング周波数。
f_s はオーバーサンプリング中のサンプリング周波数。

waveform は元になる波形データ。
*/

float f_s = L * sampleRate;

size_t N = size_t(1) << size_t(std::log2(f_s / float(10))); // 2^n.
std::vector<std::complex<float>> spc = getSpectrum<float>(waveform);

size_t M = size_t(std::log(float(N / 2)) / std::log(B));
table.resize(M + 1); // 倍音を一つも含まない最後尾のウェーブテーブルを追加。

std::vector<std::complex<float>> tmp(spc.size()); // 一時変数。
for (size_t idx = 0; idx < table.size() - 1; ++idx) {
  size_t cutoff = size_t(N / 2 * std::pow(B, -float(idx))) + 1; // c の式。

  std::fill(tmp.begin(), tmp.end(), std::complex<float>(0, 0));
  std::copy(spc.begin(), spc.begin() + cutoff, tmp.begin());

  table[idx].resize(N);
  fft.complex2real(tmp.data(), table[idx].data()); // numpy.fft.irfft と同じ計算。
}

// 最後尾のウェーブテーブルは 0 で埋める。
table.back().resize(N);
std::fill(table.back().begin(), table.back().end(), float(0));

ウェーブテーブルの切り替え

前提としてオシレータのピッチが浮動小数点数の MIDI ノート番号で入力されることにします。

事前に最低周波数 \(f_L\) を最低 MIDI ノート番号 \(\mathtt{basenote}\) に変換しておきます。

\[ \mathtt{basenote} = 12 \log_2 \left( \frac{f_L}{440} \right) - 69 \]

ピッチベンドの倍率 \(B\) も半音 \(\mathtt{interval}\) に変換しておきます。

\[ \mathtt{interval} = 12 \log_2 B \]

すると MIDI ノート番号 \(\nu\) が入力されたときに再生するウェーブテーブルのインデックス \(i\) を以下の式で計算できます。

\[ i = \left\lfloor \frac{\nu - \mathtt{basenote}}{\mathtt{interval}} \right\rfloor \]

C++ では以下のように書けます。配列の範囲外へのアクセスを防ぐために念のため std::clamp を追加しています。

class Osc {
public:
  float basenote;
  float interval;

  // 他のメンバは省略。

  // 事前計算。
  void setup(/* ... */) {
    // ...

    float f_L = f_s / float(N);
    basenote = float(12) * std::log2(f_L / float(440)) - float(69);
    interval = float(12) * std::log2(bendRange);

    // ...
  }

  // 毎サンプルの計算。
  float process(/* ... */) {
    // ...

    size_t index = std::clamp(size_t((note - basenote) / interval), 0, N - 1);

    // ...
  }
};

実装

ここでは以下の 3 つ実装を試しています。

• オーバーサンプリングありの実装
• オーバーサンプリングなしの実装
• ミップマップを使う実装

以降ではウェーブテーブルの波形にのこぎり波を使っています。音のサンプルはこの記事の最後に掲載しています。

オーバーサンプリングありの実装

オーバーサンプリングを行うウェーブテーブルオシレータのブロック線図です。

インデックス方向の補間なし

まずは De Soras さんの資料に基づいて 1 オクターブ間隔で倍音を削ったウェーブテーブルを試します。 1 つのウェーブテーブルのピッチベンドの範囲は \([1, 2)\) で、 \(B=2, H=1\) です。

以下はコードへのリンクです。

• Python 3 によるオーバーサンプリングあり、インデックス方向の補間なしの実装を読む (github.com)
• C++ によるオーバーサンプリングあり、インデックス方向の補間なしの実装を読む (github.com)

以下は MIDI ノート番号で 0, 128, 0 とピッチベンドしたテスト出力のスペクトログラムです。エイリアシングは抑えられていますが、垂直な縦線が見えるのでポップノイズが出ています。

ポップノイズはウェーブテーブルの切り替えによって発生しています。以下のリンクから確認のために書いたコードを読むことができます。 Python 3 です。

• ウェーブテーブルの切り替えがポップノイズの原因になっていることを確認するコードを読む (github.com)

インデックス方向を線形補間

ウェーブテーブルをインデックス方向で補間することでポップノイズの低減を図ります。以下はインデックス方向の線形補間を加えたコード例です。

// ウェーブテーブルオシレータのクラスメソッドの抜粋。
//
// note は MIDI ノート番号。
// maxIdx = table.size() - 2 。
//
float processSample(float note)
{
  phase += std::clamp(midinoteToFrequency(note) / sampleRate, float(0), float(0.5));
  phase -= std::floor(phase);

  auto octFloat = std::clamp((note - basenote) / interval, float(0), maxIdx);
  auto iTbl = size_t(octFloat);
  auto yFrac = octFloat - float(iTbl);

  auto pos = float(size) * phase;
  auto idx = size_t(pos);
  auto xFrac = pos - float(idx);

  auto a0 = table[iTbl][idx];
  auto a1 = table[iTbl][idx + 1];
  auto s0 = a0 + xFrac * (a1 - a0); // 位相方向の線形補間。

  auto b0 = table[iTbl + 1][idx];
  auto b1 = table[iTbl + 1][idx + 1];
  auto s1 = b0 + xFrac * (b1 - b0); // 位相方向の線形補間。

  return s0 + yFrac * (s1 - s0); // インデックス方向の線形補間。
}

以下は完全な実装へのリンクです。

• Python 3 によるオーバーサンプリングあり、インデックス方向の補間ありの実装を読む (github.com)

以下はインデックス方向を線形補間した実装で MIDI ノート番号で 0, 128, 0 とピッチベンドしたテスト出力のスペクトログラムです。ポップノイズは消えましたが、図の上部で波打つように暗くなっている箇所で高次の倍音が弱くなっています。

最適な設計

インデックス方向を線形補間したときの高次の倍音の弱まりは、ウェーブテーブルを用意するときの倍音の削り方を工夫することで防ぐことができます。

最適な設計に関するパラメータを再掲します。最低周波数 \(f_L\) はここでは使いません。

• ウェーブテーブルを \(B\) 倍までピッチベンドする。
• ウェーブテーブルはナイキスト周波数の \(H\) 倍までの周波数成分を含む。
• \(L\) 倍のオーバーサンプリング。

このとき以下の条件を満たせばエイリアシングは起こりません。

\[ BH \leq 2L - 1, \quad B \geq 1, \quad H \geq 1. \]

\(B > H\) としてピッチベンドの幅を広げると倍音の欠落が起こります。逆に \(H > B\) として倍音の余剰を増やすとピッチベンドの幅が狭くなり、テーブルの数が増えます。倍音の欠落を防ぎつつ、テーブルの数が最小となるようにするには \(H = B\) として上の不等式の上限の値と等しくなるように設計すればよさそうです。つまり \(B\) と \(H\) の値は以下の式で求められます。

\[ B = H = \sqrt{2L - 1} \]

\(L = 2\) のとき \(B = H = \sqrt{3}\) です。

以下はコードへのリンクです。

• Python 3 によるオーバーサンプリングあり、最適な設計の実装を読む (github.com)
• C++ によるオーバーサンプリングあり、最適な設計の実装を読む (github.com)

以下は最適な設計による実装で MIDI ノート番号で 0, 128, 0 とピッチベンドしたテスト出力のスペクトログラムです。テーブルの数は増えましたが、エイリアシングも倍音の弱まりもなくなっています。

ここまでやってしまえば、あとはダウンサンプリングのローパスフィルタの設計くらいでしか品質は変わりません。

オーバーサンプリングなしの実装

ウェーブテーブルの数を増やして周波数方向の間隔を狭くすることで、オーバーサンプリング無しでもエイリアシングを低減できます。この方法は CubicPadSynth で使いました。ただし CubicPadSynth は倍音のカットオフの計算が異なっているので、正しい実装は以下のリンク先のコードを参照してください。

• Python 3 によるオーバーサンプリングなしの実装を読む (github.com)
• C++ によるオーバーサンプリングなしの実装を読む (github.com)

ここでの実装は MIDI ノート番号 0 から 136 まで、 1 半音ごとに 137 個のウェーブテーブルを作っています。 MIDI ノート番号 136 は周波数に直すと約 21100 Hz です。 MIDI ノート番号 137 の周波数は 22050 Hz よりも高いので、サンプリング周波数が 44100 Hz のときにエイリアシングが起こります。

以下は MIDI ノート番号で 0, 128, 0 とピッチベンドしたテスト出力のスペクトログラムです。オーバーサンプリングしていないのでエイリアシングも倍音の弱まりも起こっているはずですが、図を見る限りではどちらも確認できません。

オーバーサンプリングがない分だけ速いので使える場面はありそうです。ただし、メモリを大量に使うのでウェーブテーブルがあまりにも長くなるとパフォーマンスが落ちます。

ミップマップを使う実装

The Quest For The Perfect Resampler で紹介されているミップマップを使う方法も実装して試しました。

メモリを節約したいときはミップマップを使うことができます。以下はウェーブテーブルのミップマップを表した図です。図の左の数字はウェーブテーブルのインデックスです。

ミップマップを使うとウェーブテーブルの間隔が 1 オクターブに固定されるので、インデックス方向の補間による倍音の弱まりが避けられない点に注意してください。

以下はミップマップを位相方向に線形補間する実装で MIDI ノート番号で 0, 128, 0 とピッチベンドしたテスト出力のスペクトログラムです。全体に薄くノイズが乗っていることが見て取れます。

以下のリンク先でミップマップと線形補間を使ったオシレータのコードが読めます。

• Python 3 によるミップマップと線形補間を使ったオシレータの実装を読む (github.com)
• C++ によるミップマップと線形補間を使ったオシレータの実装を読む (github.com)

ミップマップの線形補間は質が悪すぎるので De Soras さんの資料にならってアップサンプラを使います。

アップサンプラのフィルタ設計

以下のパラメータを決めます。

• \(f_p\): パスバンドの上限の周波数。
• \(f_s\): ストップバンドの下限の周波数。
• \(N_{\mathrm{fir}}\): ポリフェイズ 1 つあたりのタップ数。
• \(M_{\mathrm{fir}}\): アップサンプリングの倍率。

フィルタ設計中のサンプリング周波数を \(M\) にします。

\(f_p\) と \(f_s\) の値は適当に決めています。ここでは出力のナイキスト周波数は 0.5 です。またローパスなので \(f_p < f_s\) を満たす必要があります。 \(f_s\) が 0.5 より大きいとエイリアシングが出ます。しかし \(f_p\) が 0.5 より小さいと高次の倍音が弱くなります。

\(N_{\mathrm{fir}}\) と \(M_{\mathrm{fir}}\) は増やすほどフィルタの質が良くなります。 De Soras さんの資料では \(N_{\mathrm{fir}} = 12, M_{\mathrm{fir}} = 64\) としていますが、ここでは \(N_{\mathrm{fir}} = 32, M_{\mathrm{fir}} = 64\) としました。

Python3 と SciPy を使ってフィルタを設計します。

import numpy as np
import scipy.signal as signal

f_p = 0.4
f_s = 0.5
N_fir = 32
M_fir = 64

bands = np.hstack((
    [0, f_p],
    [f_s, M_fir / 2],
))
desired = (1, 0)
weight = (1, 100)
fir = signal.remez(N_fir * M_fir - 1, bands, desired, weight, fs=M_fir, maxiter=1024)

設計したフィルタをポリフェイズ分解します。各フェイズの出力振幅を 1 に正規化したいので splitPolyPhase の呼び出しで M_fir * fir としています。

import json

def splitPolyPhase(fir, nPhase):
    padding = nPhase - len(fir) % nPhase
    fir = np.hstack((fir, np.zeros(padding)))
    polyphase = []
    for phase in range(nPhase):
        part = fir[phase::nPhase]
        polyphase.append(part.tolist())
    return polyphase

fir = splitPolyPhase(M * fir, M)[::-1]

得られた fir は以下のような 2 次元配列になっています。

[
  [/* ポリフェイズ 0 の FIR 係数 */],
  [/* ポリフェイズ 1 の FIR 係数 */],
  [/* ポリフェイズ 2 の FIR 係数 */],
  // ...
]

ミップマップを fir で補間するには以下のように計算します。

// phase は [0.0, 1.0) の範囲で正規化された現在の位相。
// octave は入力されたピッチから決めたウェーブテーブルのインデックス。
Sample processTable(const Sample phase, const size_t octave)
{
  if (octave >= table.size()) return Sample(0);
  const auto &tbl = table[octave]; // table はミップマップの 2 次元配列。

  auto tblPos = Sample(tbl.size()) * phase;
  auto tblIdx = size_t(tblPos);
  auto tblFrac = tblPos - Sample(tblIdx);

  auto firPos = Sample(InterpFir::nPhase) * tblFrac;
  auto firIdx = size_t(firPos);
  auto firFrac = firPos - Sample(firIdx);

  // `tbl.size()` は常に 2^n 。よって `tbl.size() - 1` は常に `0b0..1..` 。
  auto bitmask = tbl.size() - 1;

  auto sum = Sample(0);
  for (size_t j = 0; j < InterpFir::bufferSize; ++j) {
    sum += tbl[(tblIdx + j) & bitmask] // ビットマスクでインデックスの巻き戻し。
      * (InterpFir::coefficient[firIdx][j] + firFrac * InterpFir::diff[firIdx][j]);
  }
  return sum;
}

InterpFir::diff は事前に計算しておいたフィルタの差分です。 diff[firIdx][j] = fir[firIdx][j + 1] - fir[firIdx][j] となっています。この方法は De Soras さんの資料で紹介されていました。

以下のリンク先にミップマップとアップサンプラを使ったオシレータのコードを掲載しています。

• C++ によるミップマップとアップサンプラを使ったオシレータの実装を読む (github.com)

以下のリンク先にアップサンプラのフィルタを設計する Python 3 のコードを掲載しています。

• アップサンプラのフィルタを設計するコードを読む (github.com)

以下は位相方向の補間にアップサンプラを使ったミップマップオシレータで、 MIDI ノート番号 0, 128, 0 とピッチベンドしたテスト出力のスペクトログラムです。エイリアシングがウェーブテーブルの切り替え時に少し出ていますが、悪くない品質です。図の上部が波打つように暗くなっているのは、アップサンプラのローパスフィルタのストップバンドリップルと、ウェーブテーブルのインデックス方向の補間の組み合わせによるものです。

ミップマップを使う方法はアップサンプラのフィルタの質で大きく音が変わります。 \(f_p\) を 0.5 以下にすると今回のように倍音の弱まりが出ます。かと言って \(f_p\) を 0.5 にすると \(f_s\) を 0.5 より大きい値にする必要が出てくるので、今度はエイリアシングノイズが出ます。

ポリフェイズの 1 分岐あたりの FIR の長さ \(N_{\mathrm{fir}}\) の値を大きくすることで倍音の弱まりとエイリアシングの両方を改善できますが、今度は計算が重たくなります。さらにインデックス方向の線形補間を行うときは 1 サンプルあたりの計算量が約 2 倍になります。簡単なベンチマークを取ったところ、今回実装した \(N_{\mathrm{fir}} = 32\) でインデックス方向の線形補間ありのミップマップは、オーバーサンプリングありの実装と比べて約 1.5 倍ほどの計算時間がかかりました。

ミップマップを使う方法は計算量から見るとあまり魅力的ではないですが、アップサンプラのフィルタをローパス以外の特性に変えてエフェクトをかけられるという他の実装にない特長があります。

音のサンプル

各実装でのこぎり波をレンダリングした結果です。

「オーバーサンプリングを行う実装 - インデックス方向の補間なし」は音量を上げるとテーブルの切り替え時のポップノイズが聞こえます。

「オーバーサンプリングを行う実装 - インデックス方向を線形補間」と「ミップマップ - アップサンプラ」は高次の倍音が弱まっています。「ミップマップ - アップサンプラ」では倍音の弱まりが聞き取れます。

「オーバーサンプリングを行う実装 - 最適な設計」と「オーバーサンプリングを行わない実装」はエイリアシングも倍音の弱まりもないので理想的な音です。オーバーサンプリングを行う実装はダウンサンプリングのローパスに楕円フィルタを使っているので高域の位相特性があまりよくないのですが、耳で聴きとることは困難です。

「ミップマップ - 線形補間」はエイリアシングが大量に出ているので明らかに音が違います。

その他

ウェーブテーブルの標準的な長さは Serum で使われている 2048 サンプルのようです。サンプリング周波数が 48000 Hz のとき 48000 / 2048 = 23.4375 Hz までなら高次の倍音が欠けません。

• Wavetable file format? - DSP and Plug-in Development Forum - KVR Audio

参考文献

• The Quest For The Perfect Resampler - resampler-en.pdf
• Note names, MIDI numbers and frequencies
• Fourier Series–Sawtooth Wave – from Wolfram MathWorld

変更点

• 2021/02/24
  • 構成の変更。
  • 最適な設計のウェーブテーブルを追加。
  • De Soras のミップマップウェーブテーブルを追加。
  • オーバーサンプリングを行わない実装のバグ修正。倍音のカットオフの計算に間違いがあった。
• 2021/02/25
  • \(N_{\mathrm{fir}}\) と \(M_{\mathrm{fir}}\) が \(N\) と \(M\) になっていた箇所の修正。
  • 以前のウェーブテーブルの記事名を更新。
  • オーバーサンプリングなしの実装の説明を追加。
• 2021/03/10
  • ミップマップ実装の倍音の弱まりの説明を変更。
  • InterpFir::diff の説明を追加。
  • 文章の整理。
• 2021/04/03
  • ミップマップ実装の説明でストップバンドリップルをパスバンドリップルに間違えていた箇所を修正。