同じスペクトルを持つのにクレストファクタが違う波形（Julia Advent Calendar 2023）

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

今年は、以下の論文を参考にして「同じスペクトルを持つのにクレストファクタが違う波形」を作ってみます。

Schroeder, M. R.. “Synthesis of Low-Peak-Factor Signals and Binary Sequences with Low Autocorrelation.“ IEEE Transactions on Information Theory, pp.85–89. January 1970. DOI: 10.1109/TIT.1970.1054411.

クレストファクタ

クレストファクタ（波高率）は「振幅と実効値の比」で計算されます。振幅の幅はゼロから波形ピーク値で測られます（文献によっては負側のピークから正側のピークまでとしているものもあります）。一方で実効値は波形の時間平均に近いもので、以下の式 (root mean square) で計算できます。音の波形のような交流波において振幅値の算術平均を計算すると正負の値が打ち消しあってしまうので、それを防ぐために二乗してから平均をとり平方根でもとの単位に戻すという流れで計算しています。

$\displaystyle x_\mathrm{rms} = \sqrt{\frac{1}{N}\sum_{i=1}^{N}\left(x_i\right)^2}$

もし波形が±1の範囲に正規化された正弦波（サイン波）だったとすると、振幅値のピークは1、実効値は 1/√2 になります。クレストファクタは振幅と実効値の比なので 1 ÷ 1/√2 = √2 です。正規化波形に対しての実効値は必ず1以下になるので、クレストファクタは1以上の値になります。最もクレストファクタが小さくなる（1になる）のは矩形波のように-1と+1の信号しか含まない信号です。

一般に、定常的な持続音（オルガンやクラリネットなど）のほうがクレストファクタは小さくなりがちで、逆に、大きな振幅値が時々しか現れない音（太鼓のリズムとか）はクレストファクタが大きくなりがちです。正規化波形で考えるので実際に音を聞いたときに感じる音の大きさとは異なりますが、クレストファクタは音ファイルから音量感を計算する指標の一つとして使うこともできます（たとえば以下の論文参照）。

Emmanuel Deruty and Damien Tardieu. “About Dynamic Processing in Mainstream Music.” Journal of Audio Engineering Society, vol.62, no.1/2, pp.42–55. January 2014. http://www.aes.org/e-lib/browse.cfm?elib=17085

パルス列の生成

例として余弦波の重ね合わせでパルス列を作っていくことを考えます。

$\displaystyle r(t) = \sum_{k=1}^{N}\left(\frac{p_{k}}{2}\right)^{1/2} \cos\left(\frac{2\pi k t}{T} + \theta_k\right)$

このとき、各余弦波の位相をあわせて（つまり 𝜃ₖ をすべて同じ値にして）重ね合わせると、高いクレストファクタ（論文中ではpeak-factorと書いてある）になります。ここで、𝑝ₖ は 𝑘 番目の倍音のパワーで、∑ 𝑝ₖ = 1を満たす必要があります。以下のコードでは論文中の数式に従って、第16倍音まで計算しています。

# number of harmonics
N = 16

# relative power of each harmonics
p = zeros(N)
for n in 1:N
    p[n] = 1/8 * (sin(π * (2n - 1) / 32))^2
end
p ./= sum(p)

# phase angles
ϕ = zeros(N)
for n in 1:N
    ϕ[n] = -π/2
end

# signal generation (Eq.1)
fs = 1000
T = 1
t = range(-π/2, stop=π/2, step=1/fs)
r = zeros(size(t))
for k in 1:N
    r .+= (p[k] / 2)^(1/2) * cos.(2π * k * t / T .+ ϕ[k])
end

plot(t, r, xlabel="Time", ylabel="Amplitude", legend=false, ylim=(-2.5, 2.5))

位相を調整する

測定信号や通信用の信号を生成したいときには、一般にクレストファクタは低ければ低いほどよいです。許された振幅を十分に使えるということは、信号・雑音比（S/N）を低くすることができるし、たとえばFM変調で伝送するときを考えると周波数帯域幅をまんべんなく使えるということですからね。

できるだけ全体にわたってピーク値が均等になるように、各倍音の位相を倍音の強さに応じて調整します。

$\displaystyle \phi_n = \phi_1 - 2\pi\sum_{l=1}^{n-1}\left(n-l\right)p_l$

それをプログラムにすると以下のようになります。

# phase angles
ϕ = zeros(N)
ϕ[1] = -π/2
for n in 2:N
    tmp = 0
    for l in 1:(n-1)
        tmp += (n - l) * p[l]
    end
    ϕ[n] = ϕ[1] - 2π * tmp
end

# signal generation (Eq.1)
fs = 1000
T = 1
t = range(-π/2, stop=π/2, step=1/fs)
r = zeros(size(t))
for k in 1:N
    r .+= (p[k] / 2)^(1/2) * cos.(2π * k * t / T .+ ϕ[k])
end

plot(t, r, xlabel="Time", ylabel="Amplitude", legend=false, ylim=(-2.5, 2.5))

パワースペクトルを調べると、各倍音の強さは両方とも同じになっています。位相がずれているので異なる信号なのですが、定常音についてヒトは位相に鈍感なので、実際にスピーカから再生すると同じ音に聞こえます。波形はまったく違うのに、同じ音に聞こえるって不思議ですね。