# RでCQT（Constant-Q変換）をやってみる

ホクソエムサポーターの松本です。音楽を作ったり聴いたりするのが趣味なので、音楽分析に興味があります。音データの分析にはPythonだとlibrosaというとても便利なパッケージがあるのですが、Rにはそういった汎用的なパッケージがなくてちょっと不便です。
最近ふとRでCQT（Constant-Q変換）をしてみたいと思い、既存のパッケージを使ってできないか探してみたところ特に見つからなかったので、どのように実装すればいいのか調べてみました。

## スペクトログラムについて

音声や音楽データの分析を行う際には生の波形をそのまま扱うのではなく、スペクトログラム（時間周波数表現）に変換したものを特徴量として利用することがあります。下の画像は「あいうえお」という音声を録音したデータを表したものです。

<style>
.sbs {
display: flex;
justify-content: center;
}
</style>

<ul class="sbs">
<img src="https://i.imgur.com/sWQ6WdV.png" width=55%>
<img src="https://i.imgur.com/cU78NVh.png" width=55%>
</ul>

左図の波形データは横軸は時間、縦軸は振幅を表します。右図のスペクトログラムは横軸は時間、縦軸は周波数、色はその時間・周波数の成分の強度を表します。

このように視覚化した際、波形データは音の大きさの変化くらいしか分かりませんが、スペクトログラムは周波数ごとにどのように時間変化していくか等が分かります。視覚的に特徴を捉えやすく、画像のようにも扱うことができるので画像処理の手法を応用した分析や変換にも利用しやすいといった利点があります。

スペクトログラムはSTFT（短時間フーリエ変換）によって算出するのが一般的ですが、この方法の問題点として、音高の観点からすると、低域で周波数分解能が低くなるという点があります。人間の聴覚特性として低域では周波数の違いに敏感である（周波数分解能が高い）が、高域になるにつれて鈍感になる（周波数分解能が低くなっていく）という特性があります。音高の国際的な基準はラ（A）=440Hzとなっていますが、1オクターブ上のラは880Hzでちょうど二倍の周波数になります。1オクターブの間には12個の音階があり、これらは公比 $\sqrt[12]{2}$ の等比数列となっています。

STFTでは切り取った信号に対して全ての周波数で同じ窓幅でフーリエ変換を行いますが、CQTではこうした特性を考慮して対数周波数を利用し、かつ周波数ごとに窓幅を変化させて変換することで、音高に対応して周波数分解能を変化させています。

## 実装

以下ではピアノ演奏を録音した音源を使って、各種プロットを行う例を示します。

<iframe
frameborder="0"
width="500"
height="100"
src="https://drive.google.com/file/d/1QmvJmvQhF_zIM4cMdEbAAokyPPX9G8Fg/preview">
</iframe>

音源の読み込みには tuneR パッケージ、スペクトログラムのプロットには seewave パッケージを利用しています。

library(tuneR)
library(seewave)

# wavファイルの読み込み
wav <- readWave("pf.wav")

# 波形のプロット
plot((seq(wav@left) - 1) * 1 / [e-mail], wav@left, type="l", xlab="time", ylab="amplitude")

# STFTによるスペクトログラムのプロット
ggspectro(wav, ovlp=50) + 
  scale_y_continuous() +
  geom_tile(aes(fill=amplitude))

<ul class="sbs">
<img src="https://i.imgur.com/bsXIzVP.png" width=55%>
<img src="https://i.imgur.com/dPShHcn.png" width=55%>
</ul>

<br>

次にCQTを算出してみます。実装はこちらの記事を参考にさせて頂きました。
[音楽プログラミングの超入門（仮）: 【Python】 Constant-Q 変換 (対数周波数スペクトログラム)](https://yukara-13.hatenablog.com/entry/2013/12/01/222742)

library(signal) # ハミング窓を使うためにsignalパッケージを利用

# パラメータの設定
qrate <- 1
fmin <- 60
fmax <- 6000
bins_per_octave <- 12
hop_length <- 512
two_pi_j <- 2 * pi * 1i

# 対数周波数ビンの数とそれらに対応する周波数を算出
nfreq <- round(log2(fmax / fmin) / (1 / bins_per_octave)) + 1
freqs <- fmin * (2 ** ((seq(nfreq) - 1) * (1 / bins_per_octave)))

# CQTを行うフレーム数とそれらに対応する時間を算出
nframe <- ceiling(length(wav@left) / hop_length) + 1
times <- seq(nframe) * hop_length * (1 / [e-mail])

# 窓幅を設定するためのパラメータ
Q <- (1 / ((2 ** (1 / bins_per_octave)) - 1)) * qrate

# CQTスペクトログラムのための行列を初期化
ret <- matrix(0, nframe, nfreq)

# 周波数ごとに窓幅を変えてCQT
for (k in seq(nfreq)) {
  nsample <- round([e-mail] * Q / freqs[k])
  hsample <- round(nsample / 2)
  
  phase <- exp(-two_pi_j * Q * seq(nsample) / nsample)
  weight <- phase * hamming(nsample)
  
  # 各フレームごとに信号を切り取る→重みを掛けてretに代入
  for (iiter in seq(nframe)){
    iframe <- iiter
    t <- (iframe - 1) * hop_length * (1 / [e-mail])
    istart <- (iframe - 1) * hop_length - hsample
    iend <- istart + nsample
    sig_start <- min(max(1, istart + 1), length(wav@left))
    sig_end <- min(max(0, iend), length(wav@left))
    win_start <- min(max(1, sig_start - istart), nsample)
    win_end <- min(max(0, length(wav@left) - istart), nsample)
    win_slice <- weight[win_start : win_end]
    y <- sig[sig_start : sig_end]
    ret[iiter, k] <- sum(y * win_slice) / nsample
  }
}

<br>

次に算出したCQTスペクトログラムをggplotを使ってプロットしてみます。実装はseewaveパッケージのggspectroを参考にさせて頂きました。[https://github.com/cran/seewave](https://github.com/cran/seewave)

library(ggplot2)

P <- abs(t(ret)) # 絶対値に変換し、振幅を取り出す
P <- P/max(P) # ノーマライズ
P[P == 0] <- .Machine$double.eps # log10(0)で-Infにならないようにする
P <- 20*log10(P) # dBスケールに変換

# 外れ値で表示が変にならないように範囲を設定
limit_upper <- 0
limit_lower <- -50

# ggplotで表示するためのデータフレームを作成
frequency <- rep(freqs, times=ncol(P))
time <- rep(times, each=nrow(P))
amplitude <- as.vector(P)
amplitude[amplitude < limit_lower] <- limit_lower
df <- data.frame(time, frequency, amplitude)

tlab <- "Time (s)"
flab <- "Frequency (kHz)"
alab <- "Amplitude\n(dB)\n"

# CQTスペクトログラムのプロット
ggplot(df, aes_string(x="time", y="frequency", z = "amplitude")) + 
  xlab(tlab) + 
  ylab(flab) + 
  scale_y_log10(breaks=freqs[seq(1, length(freqs), 12)]) +
  geom_tile(aes(fill=amplitude))


CQTのスペクトログラムの算出とプロットを行うことができました。

プロットを見ると比較的はっきりとピアノで演奏された音が浮き上がっていることが分かります。CQTはメロディーや和音など音楽的な要素の分析によく使われる特徴量で、音楽分析においてはSTFTを使うよりも有用な場合が多いと思います。

## 参考文献
- 高校数学の美しい物語 - 音階と数学 https://manabitimes.jp/math/1353
- 音楽プログラミングの超入門（仮） - 【Python】 Constant-Q 変換 (対数周波数スペクトログラム) https://yukara-13.hatenablog.com/entry/2013/12/01/222742
- Wizard Notes - 定Q変換 (CQT: Constant-Q Transform) の 解説（音高・コード・メロディの分析向け） https://www.wizard-notes.com/entry/music-analysis/constant-q-transform
- seewave ~ an R package dedicated to sound analysis and synthesis ~ https://rug.mnhn.fr/seewave/