音楽プログラミングの超入門(仮)

Python / 音楽情報処理 初心者が、初心者にも分かるような記事を書きたい。

【Python】 Constant-Q 変換 (対数周波数スペクトログラム)

関連記事

導入:対数周波数スペクトログラム

Pythonで短時間フーリエ変換(STFT)と逆変換 - 音楽プログラミングの超入門(仮)

上記の記事で,音響信号を周波数スペクトルの時間変化を表すスペクトログラムに変換する短時間フーリエ変換を扱いました.簡単にアルゴリズムを復習すると,音響信号を一定の幅で切り取ってフーリエ変換するという処理を少しずつずらして行っていくことでスペクトログラムを得ていました.

 ここで,音響信号を一定の幅で切り取ってフーリエ変換するということについて少し考えてみましょう.切り取られた信号が以下のようなものであることを考えます.

f:id:yukara_13:20131201171527p:plain

窓幅が 256 サンプル(fs=4[kHz])で、ここでは信号を 1000[Hz],400[Hz],100[Hz],20[Hz] のサイン波を足し合わしたものにしています。これを見ると、1000[Hz]・400[Hz] は十分過ぎるくらいサンプルが取れており、20[Hz] ではちょっとサンプルが足りなさそう、100[Hz] はちょうどいいくらいかなって感じになってますね。短時間フーリエ変換(STFT)の記事で述べた分解能に即して考えると、高周波数では十分な周波数分解能にも関わらず、低周波数では十分でない周波数分解能となっているのです。窓幅が全ての周波数で同じだからこのようなことが起こるのですね。

 じゃあ各周波数に適した分解能を使う、つまり、各周波数で窓幅を変えたらいいんじゃない?って思いますよね。これを実現するのが Constant-Q 変換 (定Q変換) です。例えば、各周波数について丁度 5 周期分とるように窓幅を決めてやると、以下の図のようになります。

f:id:yukara_13:20131201180159p:plain

サンプリング周波数が小さいので 1000[Hz] のサイン波がアレな感じになってますが、大体全部同じ形です.当然ながら、周波数ごとに窓幅が変わっているのがお分かりいただけると思います。

 Constant-Q 変換では、どの周波数でも含まれる周期数が同じになるように窓幅を変えるので、周波数と窓幅は反比例の関係になることが分かります。すると必然的に、周波数の下限を決めてやらなければならないことも分かりますね。このことからさらに、Constant-Q 変換は不可逆変換だということも分かってしまいました(泣)。

 ここで対数周波数というものを導入しましょう。といっても特に難しいものでもなく、文字通り、今まで線形で扱っていた周波数軸 (単位が Hz) をただ対数にしたもののことです。周波数 F1 [Hz] と周波数 F2 [Hz] があったときに、線形周波数では F1 と F2 の差を重視し、対数周波数では F1 と F2 の比率(何倍か)を重視した軸であると言えますね。対数周波数では、倍率が同じであればその距離が同じになるんですね。例えば、200 [Hz] と 400 [Hz]、300 [Hz] と 600[Hz] の周波数を考えたとき、線形周波数では当然それぞれ 200 [Hz] と 300 [Hz] の差ですが、対数周波数では、



\begin{align}
\displaystyle
 \log 400 - \log 200 &= \log (400/200) = \log 2 \\
 \log 600 - \log 300 &= \log (600/300) = \log 2
\end{align}

という感じで、同じになりますね。フーリエ変換で得られるような線形周波数を対数周波数に変換すると、低周波数部がびよーんと伸び分散の大きな不格好なものとなってしまいます。一方、Constant-Q 変換では対数周波数領域において高周波でも低周波でも綺麗なピークの立った対数スペクトルが得られます。

f:id:yukara_13:20140126202459p:plain

上図は、上からフーリエ変換フーリエ変換対数スケールにしたもの、Constant-Q 変換のスペクトルを表しています。下の二つは周波数を合わせていますが、フーリエ変換では周波数が低くなるに従って調波の分散が大きくなっているのに対し、Constant-Q 変換では全ての周波数において分散が一定である様子が分かります。

 元論文を読むと、Constant-Q 変換 は音楽信号の解析用に作られたアルゴリズムで、上記のように倍率が同じ時に距離が同じになることやシンプルなところがウェーブレット変換とかよりええんやで、みたいに書いてありますね。人間は音の変化を対数的に捉えるので(1オクターブで2倍に変化とか)、そういう意味で相性がいいんです。それでは、具体的にConstant-Q 変換の定式化・アルゴリズムを見ていきましょう。

Constant-Q 変換

以下が元論文です。
    Judith C. Brown : Calculation of a constant Q spectral transform, J. Acoust. Soc. Am., 89(1):425–434, 1991.
まずはぴゃぴゃっと定式化がどのようになっているか見ていきましょう。



\displaystyle
\begin{align}
X^{\mathrm{cq}}[k]=\frac{1}{N[k]}\sum_{n=0}^{N[k]-1} w[ n,k] \ x[ n]\ \exp(-2j\pi Q n/N[k])
\end{align}

なんやこれ・・。と思われた方もいるかもしれませんが、よく見ると全く難しくないことが分かります。これは、フーリエ変換の式とほとんど同じなのです。 $k$ がいわゆる対数周波数で、STFT では一定だった窓幅  $N$ が、この式では  $N[ k]$ というように周波数によって変わっていることが分かります。ここで、



\begin{align}
$\displaystyle
 N[ k] &= Q\frac{f_s}{f_k} \\
 Q &= (2^{\frac{1}{\mathrm{fratio}}} - 1)^{-1} \;\times\; \mathrm{qrate}
\end{align}

です。  $f_s$ $f_k$ はそれぞれサンプリング周波数、対数周波数が  $k$ のときの線形周波数を表しています。 $Q$ ってなんやねんと思った人も多いでしょうが、これは簡単にいうとどれくらいの周期数を含むように窓幅を設定するかを決めるパラメータで、大きいほど窓幅も全体的に大きくなります。 $\mathrm{fratio}$オクターブを幾つに分割するかを表したパラメータで結構重要です。 $\mathrm{qrate}$ は(確か)元論文には登場せず、通常 1 ですが、これだと窓幅が大きすぎて計算量的にも辛いです。自分がいつもどのように設定しているかは下の実装のとおりなので参考にしてください。その設定で困ったことはまだありません。

実装

それでは、これを実装していきましょう!!今回は、上に書いた式をそのまんま計算するだけなので、目で追っていけば簡単に理解できると思います。

# -*- coding: utf-8 -*-
"""
Constant-Q Transform
"""
from scipy import arange, complex128, exp, hamming, log2, zeros
from scipy import pi as mpi
from scipy.io.wavfile import read

from matplotlib import pylab as pl

# ========
#  HELPER
# ========

# 対数周波数ビン数を計算
def get_num_freq(fmin, fmax, fratio):
    return int(round(log2(float(fmax) / fmin) / fratio)) + 1

# 各対数周波数ビンに対応する周波数を計算
def get_freqs(fmin, nfreq, fratio):
    return fmin * (2 ** (arange(nfreq) * fratio))


# ======
#  MAIN
# ======
two_pi_j = 2 * mpi * 1j

### default parameters
fmin_default = 60 # min of frequency
fmax_default = 6000 # max of frequency
fratio_default = 1. / 24 # fratio (24 freq bins per 1 octave)
q_rate_def = 20. * fratio_default # qrate

def cq(sig, fs, q_rate = q_rate_def, fmin = fmin_default, fmax = fmax_default, 
       fratio = fratio_default, win = hamming):
    # 100 frames per 1 second
    nhop = int(round(0.01 * fs))

    # Calculate Constant-Q Properties
    nfreq = get_num_freq(fmin, fmax, fratio) # number of freq bins
    freqs = get_freqs(fmin, nfreq, fratio) # freqs [Hz]
    Q = int((1. / ((2 ** fratio) - 1)) * q_rate) # Q value
    
    # Preparation
    sig_len = len(sig)
    nframe = sig_len / nhop # number of time frames

    ret = zeros([nframe, nfreq], dtype = complex128) # Constant-Q spectrogram

    # Execution
    for k in xrange(nfreq):
        print k + 1, "/", nfreq 
        freq = freqs[k]
        nsample = int(round(float(fs * Q) / freq))
        hsample = int(nsample / 2)

        # Calculate window function (and weight).
        phase = exp(-two_pi_j * Q * arange(nsample, dtype = float) / nsample)
        weight = phase * win(nsample)

        # Perform Constant-Q Transform.
        for iiter in xrange(nframe):
            iframe = iiter
            istart = iframe * nhop - hsample   # i番目のフレームがどこから始まるか.
            iend = istart + nsample            # i番目のフレームがどこで終わるか.
            sig_start = min(max(0, istart), sig_len)
            sig_end = min(max(0, iend), sig_len)
            win_start = min(max(0, sig_start - istart), nsample)  # 最初らへん対策
            win_end = min(max(0, sig_len - istart), nsample)      # 最後らへん対策
            win_slice = weight[win_start : win_end]
            y = sig[sig_start : sig_end]
            ret[iiter, k] = (y * win_slice).sum() / nsample

    return ret, freqs


from scipy import interp
if __name__ == "__main__":
    wavfile = "../wav/aiueo.wav"
    fs, data = read(wavfile)

    ### Constant-Q Transform
    cq_spec, freqs = cq(data, fs)
    w, h = cq_spec.shape

    ### Plot
    fig = pl.figure()
    fig.add_subplot(111)
    pl.imshow(abs(cq_spec).T, aspect = "auto", origin = "lower")
    pl.yticks(arange(h)[::24], freqs[::24], fontsize = 15)
    pl.xticks(fontsize = 15)
    pl.ylabel("Frequency [Hz]", fontsize = 20)
    pl.xlabel("Time [10 msec]", fontsize = 20)
    pl.tight_layout()
    pl.show()

特に難しいところはないですね。モノラル信号しか動かないので、ステレオをやりたい場合は適当に改良してください。今回は元論文に従って、周波数ビンは1オクターブにつき 24 となっています。得られる対数スペクトログラムはこんな感じです。

f:id:yukara_13:20131201215844p:plain

あぁ~対数ぅ~。って感じが見てとれるでしょう。

高速化について

上のプログラムを実際に動かしてみると分かりますが、すごく遅いです。まあ式どおりに for 文で愚直に計算しているので当然です。

 実は FFT を用いて、Constant-Q 変換を高速に行うアルゴリズムが存在します。

    Judith C. Brown and Miller S. Puckette: An efficient algorithm for the calculation of a constant Q transform, J. Acoust. Soc. Am., 92(5):2698–2701, 1992.

このアルゴリズムは次回解説しようと思います。こっちのはかなり速いです。最初からこれを載せろと思うかもしれませんが、アルゴリズムを理解するためにも愚直に実装することも大事かなあと思いました(適当)。

 また、マルチレート信号処理みたいなことをして高速化をはかるアルゴリズムもあります。

    Christian Schorkhuber and Anssi Klapuri: Constant-Q Transform Toolbox for Music Processing, SMC Conference 2010, 2010.

こちらも実装してみたんですが、遅くなった上に、微妙に変な結果になってしまったので、とりあえず放置しています。多分、一つ目のアルゴリズムが十分に速いので、もう手を付けることはないでしょう・・。

逆変換について

最初の方に Constant-Q 変換は不可逆変換であると書きました。それはそうなのですが、何とかして近似的な逆変換をしようという試みは色々あります。

    Derry Fitzgerald and Matt Cranitch and Marcin T. Cychowski: Towards an Inverse Constant Q Transform, An Audio Engineering Society Convention Paper presented at the 120th. Convention, Paris, France, May 20th. - 23rd., 2006.

これは、一旦 STFT 領域に変換してから iSTFT することで元の信号に逆変換するそうです。かつて、実装を試みましたが、バグか何かで動かなかったので、それ以来放置しています・・。いつか成功させたいです・・。

 ちなみに、Brown さんは擬似逆行列を用いた逆変換を行っています。このアルゴリズムは非常に簡単ですが、質はかなり低いです。あまり実用的ではない気がします。

まとめ

Constant-Q 変換は色々と便利なので、是非遊んで見てください。