トゥルーピーク～デジタルデータのサンプル間は直線ではない～デジタルの誤解を解く

2023年12月10日 01:43

まえがき

オーディオ系のお話です。

この記事は Kindle 電子書籍
「MATLAB で簡単オーディオプラグイン開発 : エフェクター・プラグインが数分～で作れる！（VST/AU）」
の「サンプリング（標本化）定理」と「トゥルーピーク」の項を抜粋して再編集したものです。

サンプリング（標本化）定理

サンプリング周波数を Fs とすると、「最高周波数（広義には信号帯域）が Fs/2 未満の信号であれば、サンプリング周波数 Fs 以上の標本化によって情報の消失なしに離散化できる」、というのがサンプリング（標本化）定理です。（「シャノン-染谷の標本化定理」とも呼ばれる）

詳しくは、「折り返し雑音とは？～なぜ車のホイールは逆回転するのか～デジタルはアナログの近似ではない」をご参照ください。

アナログ <-> デジタル信号の変換で、サンプリング周波数を上げればいくらでも精度が上がりそうに思えるかもしれませんがそうではなく、元信号帯域がある範囲に収まっていれば、その倍を超える周波数でサンプリングを行えば、（原理的には）アナログ信号とデジタル信号は一意に変換できると言うことです。

別の言い方をすると、たまに見かける「デジタル信号はデータが時間的に飛び飛びなので、その間は適当に繋ぐしかなく、元のアナログ信号とは異なる」というのは誤解、ということになります。

間の値は、適当な値でも直線でもなく、$${sinc}$$ 関数 $${sin(x)/x}$$（の積分）で決まるただ一つの曲線であることが理論的に導出されます。

上のリンク先の記事内で詳しく説明していますが、$${sinc}$$ 関数は理想 LPF 特性の逆フーリエ変換で、D/A 変換により補間されます。

例えば 44.1 kHz サンプリングでの 18 kHz 正弦波のサンプル点を、plot 関数を用いてプロットすると以下のようになります。

このサンプルレートでは一周期にサンプル点が２～３点しかありません。plot 関数はサンプル間を直線で繋ぐため、あたかも正弦波が表現できていないかのように見えてしまいます。

しかし上で述べたように、サンプル間は直線ではなく $${sin(x)/x}$$ の積分です。

理論的には無限区間積分ですが、注目範囲前後で積分してみましょう。

Fs = 44100;
t0 = (0:1/Fs:1)';
y18 = sin(2*pi*18e3*t0);

st0 = 17228;  % 注目区間
en0 = 17236;  % 注目区間
st = st0-20;  % 積分区間
en = en0+20;  % 積分区間
x = st:en;

yc = y18(st:en);
ts = linspace(st-10,en+10,800);  % 補間サンプル点

[Ts,T] = ndgrid(ts,x);  % コンボルーション用2次元グリッド生成
y = sinc(Ts - T) * yc;  % sinc 補間

plot(x,yc,'square',MarkerFaceColor='b')
hold on
plot(ts,y,'b')
hold off
grid on
xlim([st0-0.5 en0+0.5])
ylim([-1 1])

このように、18 kHz 正弦波が表現できているのが分かります。

また、一部の DAW 等では、拡大時に補間されるものもあります。

現実的には、理想 LPF は無限タップ長となるため完全には作れませんしデバイスの精度等色々と問題もありますが、原理的な問題とは切り離して考える必要があります。

また、こういったサンプリングに関わる問題は、Fs の数倍のサンプリングレートを用いて変換するオーバーサンプリングにより影響を軽減することができます。

トゥルーピーク

サンプリング定理より、デジタル信号は基本的に Fs/2 以上の周波数成分を持つことができないことが分かります。

そのため、デジタルデータを直接操作すると、帯域制限をせずにサンプリングしたときに発生する折り返し雑音相当の成分が発生することがあります。

例えば振幅が 1.2 の 10 kHz 正弦波を ±1.0 でハードクリップさせた場合のスペクトルを見てみましょう。

MATLAB には簡単にスペクトル／スペクトログラムを表示することができる spectrumAnalyzer 関数が用意されています（DSP System Toolbox）。
（古いバージョンは、アップデート毎に波形が異なるので注意　古いリリースでは、高調波が多めに出るようです　最新のパッチを当てておきましょう）

Fs = 44100;
t = (0:1/Fs:1)';  % 1秒分の列ベクトル
% ' は行ベクトルからの転置
% 転置しなくてもよいが、オーディオ信号では
% 列はチャンネルとして扱うため
f = 10e3;
y = 1.2*sin(2*pi*f*t);  % 1秒分の信号生成

% スペクトルアナライザー・オブジェクト生成
scope1 = spectrumAnalyzer(SampleRate=Fs, ...
    SpectrumUnits="dBFS", ...
    YLimits=[-80 0], ...
    PlotAsTwoSidedSpectrum=false, ...
    Title='Non clip');
% スペクトログラムも表示する場合
% scope.ViewType = "spectrum-and-spectrogram";

scope1(y)  % スペクトルアナライザー表示
scope2 = clone(scope1);  % オブジェクトを設定ごとコピー
scope1.PeakFinder.Enabled = true;  % ピーク検出

% hard clip
yd = min(y,1);
yd = max(yd,-1);

scope2.Title='Hard clip';  % タイトルのみ変更
scope2(yd)
scope2.PeakFinder.Enabled = true;

このように、本来の周波数以外の成分が多く発生し、「ビ～」という濁った音に聞こえます。

spectrumAnalyzer のスクリプト中からのピーク検出は、タイミングによって出ないことがあります。その場合は実行行を変えてみるか、一旦 false にしてから true にし直すか、メニューから行ってください。

一旦閉じてしまった Window を再表示したい場合は、show() を使います。

>> show(scope2)

参考までに、Limiter、Compressor、簡易的なリミッターとして用いられることもある tanh を使用した場合のスペクトルも下に示しておきます。
Limiter、Compressor の設定は以下の通りです。

Limiter：Threshold=0dB, AttackTime=0s, ReleaseTime=0.2s
Compressor：Threshold=-1.4dB, Ratio=5, AttackTime=0.05s, ReleaseTime=0.2s

また、上で述べたようにデジタルのサンプル間は直線ではなく $${sin(x)/x}$$ の積分であるため、実際の波形は plot 関数でプロットした波形とは違っています。

**実際の波形**（区間積分）
（サンプル点間が sinc(x) の積分）
赤の横線が ± 1.0

これを見て分かるとおり、デジタルデータ自体は ±1 に収まっていますが実際の信号振幅のピークはそれを越えています。

つまり、デジタルデータでは ±1 に収まっていても、アナログ信号に戻したときに ±1 を超えることがあります。これをトゥルーピークと呼びます。

loudnessMeter 関数を用いると、内部でオーバーサンプリングしてトゥルーピークの推定値を求めることができます。

> loudMtr = loudnessMeter('SampleRate',Fs);
> [~,~,~,~,peakLDN] = loudMtr(yd);
> db2mag(peakLDN)

ans =
    1.2193

つまり、デジタルデータをクリップさせるとデジタルデータとして「歪む」だけでなく、アナログ信号に直したときにさらに歪む可能性が大きく、しかもその歪み方は再生系に依存しコントロールできないということです。

デジタル放送等に関する標準規格策定を行っている ARIB では、ラウドネス測定に関する国際規格 ITU-R BS.1770 に基づき、4倍オーバーサンプリングを行ってそのピークを -1 dB とする（48 kHzサンプリング時）などの運用規定が勧告（ARIB TR-B32）されています。

最後に、波形的にも見やすいように、1 kHz 正弦波の例も示しておきます。