【Python】高速な Constant-Q 変換 (with FFT)

Constant-Q 変換
【Python】 Constant-Q 変換 (対数周波数スペクトログラム) - 音楽プログラミングの超入門（仮）
疎行列の保存：カーネル行列を保存するのに必要です
Pythonで疎行列を保存する方法 - 音楽プログラミングの超入門（仮）

Constant-Q 変換が遅い件

上記の記事で、音響信号を対数周波数スペクトログラムに変換する Constant-Q 変換を扱いました。実装もしましたが、for 文をひたすら回すので、計算にとても時間がかかりました。
　そこで、今回はこれを高速化するという試みです。具体的には以下の論文のアルゴリズムを実装します。

An efficient algorithm for the calculation of a constant Q transform

この論文のアルゴリズムは、簡単に言うと、うまく数式を変換して高速フーリエ変換(FFT)を使えるようにしようというものです。FFT はかなり強力なアルゴリズムなので、何か速くなりそうですよね。

高速化アルゴリズム

パーセバルの等式

このアルゴリズムでは、パーセバルの等式というものを用いた数式変換が鍵となってくるので、まずこれを説明します。これは以下の等式を指します。

$\displaystyle \begin{align} \sum_{n=0}^{N-1}x[n]y^*[n] = \frac{1}{N}\sum_{k=0}^{N-1}X[k]Y^*[k] \end{align}$

ここで、 $X[k]$ と $Y[k]$ はそれぞれ、 $x[n]$ と $y[n]$ のフーリエ変換で、また、 $Y^*[k]$ は $Y[k]$ の複素共役です。

数式の変換

まず、Constant-Q 変換は下式で表せました。（詳しくは前の記事を読んでください。）

$\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}$

まず、この式を行列演算に変換してみましょう。Constant-Q 変換を音響信号 x に変換行列 T をかけ合わせる表現で表すイメージです。

$\displaystyle \begin{align} X^{\mathrm{cq}}[:] = x\cdot T^* \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; (T \ \ \in \ \ N * K) \end{align}$

ここで

$\displaystyle \begin{align} T_{nk}&= \left \{ \begin{array}{l} \frac{1}{N[k]}w[n,k]\exp(2j\pi Q n/N[k]) \;\;\;\;\; (if \ n< N_k) \\ 0 \;\;\;\;\; (otherwise) \end{array} \right. \end{align}$

こんな感じです。 $N$ は最大の $N[k]$ (最小の周波数に対応する)よりも大きい値であることに注意しましょう。積和の行列表現はよく使いますね。SciPy による高速計算を行う上でも重要です。
　さて、ここで先ほどの式をパーセバルの等式に代入すると以下の式が得られます。

$\displaystyle \begin{align} X^{\mathrm{cq}}[:] = \frac{1}{N}X\cdot S^* \end{align}$

ここで

$\displaystyle \begin{align} X &:=& \ \ DFT(x) \;\;\;\;\;\;\;\; (\in N)\\ S &:=& \ \ DFT(T) \;\;\;\;\;\;\;\; (\in N*K) \end{align}$

この $S$ をカーネル行列と名付けましょう。

行列のスパース性

さて、こんなわけで、Constant-Q 変換を FFT(DFT) を使って計算できる形に変形できました。しかし、数式をみると結局、長さ N のベクトルと、大きさ N*K の行列の乗算を計算しなければならず、それって結局計算重いんじゃないの？と思った人もいると思います。
　ここが（多分）このアルゴリズムのポイントとなってくるのですが、実は、上のカーネル行列は疎行列(スパース行列)になっているのです。疎行列とは、ほとんどの要素値が 0 で占められているような行列のことをいいます。単純な行列積では、全ての要素の組み合わせについて乗算を行わなければなりませんが、疎行列ではそのほとんどを計算せずに済むため、計算量が大幅に削減されます。そして、SciPy には scipy.sparse という、疎行列を扱うためのパッケージが含まれています！！
　本当にカーネル行列って疎行列になってるの？という人のために、実際に計算した一例を載せるとこんな感じです。

f:id:yukara_13:20140105044339p:plain,w400

この青い部分が全部 0 なので、ほんとに一部しか値が入っていないのが分かりますね。

カーネル行列は事前に用意できる

変形したあとの式を実装していけばいいのですが、ここで重要なのが、カーネル行列の計算が入力音響信号に依存していない点です。つまり、Q値や周波数分解能などのパラメータを変える必要がなければ、カーネル行列を事前に計算して保存しておけるということです。実際に計算する部分は、音響信号の FFT と最後の行列積のみとなります。

実装

では実装していきましょう。この実装にはカーネル行列の計算も含んでいますが、事前に保存してあるものを使う場合は、適当に改良してください。

# -*- coding: utf-8 -*-
"""
Efficient Constant-Q Transform (with FFT)
"""
from scipy import arange, ceil, complex128, dot, exp, float64, hamming, log2, zeros
from scipy import pi as mpi
from scipy.fftpack import fft
from scipy.sparse import lil_matrix, csr_matrix
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_fft(sig, fs, q_rate = q_rate_def, fmin = fmin_default, fmax = fmax_default, 
           fratio = fratio_default, win = hamming, spThresh = 0.0054):
    # 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
    L = len(sig)
    nframe = L / nhop # number of time frames

    # N  > max(N_k)
    fftLen = int(2 ** (ceil(log2(int(float(fs * Q) / freqs[0])))))
    h_fftLen = fftLen / 2
   
    # ===================
    #  カーネル行列の計算
    # ===================
    sparseKernel = zeros([nfreq, fftLen], dtype = complex128)
    for k in xrange(nfreq):
        tmpKernel = zeros(fftLen, dtype = complex128)
        freq = freqs[k]
        # N_k 
        N_k = int(float(fs * Q) / freq)
        # FFT窓の中心を解析部分に合わせる．
        startWin = (fftLen - N_k) / 2
        tmpKernel[startWin : startWin + N_k] = (hamming(N_k) / N_k) * exp(two_pi_j * Q * arange(N_k, dtype = float64) / N_k)
        # FFT (kernel matrix)
        sparseKernel[k] = fft(tmpKernel)

    ### 十分小さい値を０にする
    sparseKernel[abs(sparseKernel) <= spThresh] = 0
    
    ### スパース行列に変換する
    sparseKernel = csr_matrix(sparseKernel)
    ### 複素共役にする
    sparseKernel = sparseKernel.conjugate() / fftLen
 

    # ===========
    #  Execution
    # ===========
    ### New signal (for Calculation)
    new_sig = zeros(len(sig) + fftLen, dtype = float64)
    new_sig[h_fftLen : -h_fftLen] = sig
    
    ret = zeros([nframe, nfreq], dtype = complex128)
    for iiter in xrange(nframe):
        print iiter + 1, "/", nframe
        istart = iiter * nhop
        iend = istart + fftLen
        # FFT (input signal)
        sig_fft = fft(new_sig[istart : iend])
        # 行列積
        ret[iiter] = sig_fft * sparseKernel.T
 
    return ret, freqs


# ======
#  Test
# ======
if __name__ == "__main__":
    wavfile = "../wav/aiueo.wav"
    fs, data = read(wavfile)

    ### Constant-Q Transform
    cq_spec, freqs = cq_fft(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()