定Q変換（CQT）の理論と導出と実装

STFT（左）では全周波数で帯域幅が固定なため、低音域（赤いゾーン）では帯域が広すぎて音程が潰れてしまいます。一方、定Q変換（右）はQ値（= 中心周波数 / 帯域幅）を一定にし、どの音域でも「相対的な細かさ」が均等になっています。この違いが、音楽信号解析でCQTが選ばれる根本的な理由です。

ピアノで一番低い「ラ」（A0, 約27.5 Hz）と、その1オクターブ上の「ラ」（A1, 約55 Hz）の音を区別したいとしましょう。2つの音は周波数にしてわずか27.5 Hzしか離れていません。一方、高い方の「ラ」（A7, 約3520 Hz）とそのオクターブ上（A8, 約7040 Hz）は、なんと3520 Hzも離れています。同じ「1オクターブ」という音楽的な隔たりなのに、周波数で測るとその間隔は128倍も違うのです。

ここで普通の短時間フーリエ変換（STFT）を使うと困ったことになります。STFTは全周波数で同じ周波数分解能（同じ刻み幅）を持つため、低音側の27.5 Hzの差を見分けようと分解能を細かくすると、高音側では無駄に細かすぎて時間方向の応答が鈍くなり、逆に高音に合わせると低音側のオクターブがまったく潰れて見えてしまいます。「低い音ほど周波数の差を細かく、高い音ほど時間の変化を素早く捉えたい」——この人間の聴覚や音楽の構造に沿った要求に応えるのが、本記事で扱う定Q変換（Constant-Q Transform, CQT）です。

定Q変換は、音楽情報処理（楽譜採譜、コード認識、楽器音解析）の標準ツールであるほか、対数的な周波数構造を持つあらゆる信号——たとえば動物の発声解析、機械の回転に伴う倍音構造の振動診断、地震波の周期解析——に広く使われています。STFTでは扱いにくかった「対数スケールでの一様性」を最初から組み込んだ変換であり、フーリエ解析を学ぶうえで「分解能とは何か」「窓と周波数の関係」を深く理解する格好の題材でもあります。

本記事の内容

• STFTの固定分解能という限界と、Q値（品質係数）という発想の直感的な理解
• 「Qを一定に保つ」という要請から、各周波数ビンの窓長 $N_k = Q \cdot f_s / f_k$ を導出
• カーネル行列を使った効率的な実装方法の数式と意味
• 対数掃引チャープにSTFTとCQTを適用し、対数周波数軸上での見え方の違いをPythonで可視化

前提知識

この記事を読む前に、以下の記事を読んでおくと理解が深まります。

• 短時間フーリエ変換（STFT）とスペクトログラムを理解する — 本記事はSTFTの「固定分解能」という限界を出発点とします
• 窓関数（ハニング窓・ハミング窓）の役割を理解する — CQTでも各ビンに窓をかけるため、窓と周波数分解能の関係が前提になります

STFTの限界とQ値という発想

まず、なぜ新しい変換が必要なのかを、STFTの仕組みに立ち返って確認しましょう。

STFTは、信号 $x[n]$ に有限長の窓 $w[n]$ をかけて切り出し、その区間でフーリエ変換を行う操作を、窓を少しずつずらしながら繰り返すものでした。窓の長さを $N$ サンプルとすると、得られる周波数の刻み幅（周波数分解能）は

$$ \Delta f = \frac{f_s}{N} $$

で与えられます。$f_s$ はサンプリング周波数です。ここで決定的に重要なのは、STFTでは窓長 $N$ が全周波数で共通だということです。つまり $\Delta f$ も全周波数で一定です。たとえば $f_s = 44100$ Hz、$N = 4096$ なら $\Delta f \approx 10.8$ Hz で、これは低音でも高音でも変わりません。

この「全周波数で一定の $\Delta f$」が、冒頭のピアノの例で問題を起こします。低音 A0（27.5 Hz）とその半音上 A#0（約29.1 Hz）の差はわずか1.6 Hzですから、これを見分けるには $\Delta f$ を1 Hz以下にする必要があり、窓長は数万サンプル（1秒近く）になります。ところが高音 A7（3520 Hz）とその半音上（約3729 Hz）の差は約209 Hzもあるので、こんなに細かい分解能はまったく不要で、むしろ長い窓のせいで「いつ音が鳴ったか」という時間情報がぼやけてしまいます。

ここで音楽家や聴覚の視点を取り入れます。私たちが「音程」として感じるのは周波数の差ではなく比です。どの高さでも1オクターブは周波数比2倍、半音は比 $2^{1/12}$ 倍です。だとすれば、分解能も「絶対的な刻み幅 $\Delta f$」ではなく、「中心周波数に対する相対的な細かさ」で測るべきではないか——この相対的な細かさを表すのがQ値（品質係数, quality factor）です。

Q値は、ある周波数フィルタについて、その中心周波数 $f$ と帯域幅 $\Delta f$ の比として定義されます。

$$ Q = \frac{f}{\Delta f} $$

Q値が大きいほど「中心周波数の割に帯域が狭い」=「鋭く選択的な」フィルタを意味します。STFTでは $\Delta f$ が一定だったので、$Q = f / \Delta f$ は周波数 $f$ に比例して増えていきます。つまりSTFTは低周波でQが小さく（甘く）、高周波でQが大きく（鋭く）なる——音楽の感覚とは正反対の振る舞いをしていたのです。

そこで発想を逆転させます。「$\Delta f$ を一定にする」のをやめ、「Q を一定にする」ことを設計目標にしたらどうなるでしょうか。これがまさに定Q変換の核心です。次のセクションで、この一文を数式に落とし込んでいきます。

このグラフは、STFTのQ値（青線）が周波数に比例して増加するのに対し、定Q変換のQ値（オレンジ破線）が16.8付近に一定に保たれることを示しています。赤い塗りつぶし（低音域）はSTFTのQ不足ゾーンで、音程が潰れやすい。青い塗りつぶし（高音域）はQ過剰ゾーンで、時間分解能が犠牲になっている。定Q変換はこの両問題を同時に解決する設計です。

定Q変換とは — Qを一定にするという設計

定Q変換を一言で言えば、「各周波数ビンごとに窓の長さを変え、どの周波数でもQ値（中心周波数と帯域幅の比）が同じになるようにしたフーリエ変換」です。STFTがすべて同じ長さの窓で信号を覗いていたのに対し、CQTは低い周波数を見るときは長い窓（時間をかけてじっくり）、高い周波数を見るときは短い窓（さっと一瞬で）を使い分けます。

イメージとしては、虫眼鏡のセットを思い浮かべると分かりやすいでしょう。STFTは倍率が固定の虫眼鏡1本で、近くの大きな文字（高周波）も遠くの小さな文字（低周波）も同じ倍率で見ようとします。一方CQTは、対象の大きさに応じて倍率を変える可変ズームのようなもので、「小さい対象（低周波の細かい音程差）は高倍率で、大きい対象（高周波の大まかな構造）は低倍率で」見る、という賢い覗き方をします。

この「Qを一定に」という要請には、もう一つ自然な帰結があります。Qが一定なら、隣り合うビンの中心周波数の比も一定になります。すなわちCQTの周波数ビンは、STFTのように $0, \Delta f, 2\Delta f, \dots$ と等差数列（線形）ではなく、$f_0, f_0 r, f_0 r^2, \dots$ と等比数列（対数的）に並びます。音楽でいえば、1オクターブを等分（たとえば12半音）した格子にビンを置くことに対応します。これが「CQTのスペクトログラムは対数周波数軸上で音楽的に等間隔になる」という大きな利点の源です。

それでは、この言葉による設計を数式に落とし込み、各ビンの窓長を具体的に求めていきましょう。

左図はビン番号 $k$ に対して中心周波数 $f_k$ を片対数軸で描いたもので、指数関数が直線に乗ることを確認できます。右図は対数軸上でのビン密度を比較したもので、CQT（橙）は対数軸で均等に並んでいるのに対し、STFT（青）は低音域でビンがスカスカになっていることが一目でわかります。この均等性こそ、CQTを音楽信号や倍音構造の解析に適した変換にしている性質です。

周波数ビンと窓長の導出

定Q変換を構成するには、(1) どの周波数にビンを置くか、(2) 各ビンの窓をどれだけの長さにするか、の2つを決める必要があります。順に導出します。

ステップ1：対数的な周波数ビンの配置

CQTのビンは対数的に等間隔、すなわち等比数列で並べます。1オクターブ（周波数比2倍）を $b$ 個のビンに分割するとしましょう。音楽用途では半音単位なら $b = 12$ がよく使われます。すると、隣り合うビンの周波数比 $r$ は、$b$ 個進むと2倍になる、という条件から決まります。

$$ r^{b} = 2 \quad \Longrightarrow \quad r = 2^{1/b} $$

最低周波数を $f_{\min}$ とすると、$k$ 番目（$k = 0, 1, 2, \dots$）のビンの中心周波数は

$$ \begin{equation} f_k = f_{\min} \cdot 2^{k/b} \end{equation} $$

となります。最高周波数 $f_{\max}$（ナイキスト周波数 $f_s/2$ 以下に取る）までに含まれるビンの総数 $K$ は、$f_{\min} \cdot 2^{K/b} \le f_{\max}$ という条件から

$$ K = \left\lceil b \cdot \log_2 \frac{f_{\max}}{f_{\min}} \right\rceil $$

と求まります。$\lceil \cdot \rceil$ は天井関数（切り上げ）です。たとえば $f_{\min} = 32.7$ Hz（C1）、$f_{\max} = 16000$ Hz、$b = 12$ なら、$\log_2(16000/32.7) \approx 8.93$ オクターブ分、$K \approx 108$ 個のビンとなり、これはピアノの鍵盤数（88鍵）に近い分割になります。

ステップ2：Q値の決定

次にQ値そのものを決めます。Qは「中心周波数 $f_k$ と、そのビンが受け持つ帯域幅 $\Delta f_k$ の比」 $Q = f_k / \Delta f_k$ でした。ここで自然な要請は、「隣のビンとちょうど接する程度の帯域幅にする」ことです。ビン $k$ と隣のビン $k+1$ の中心周波数の差は

$$ \Delta f_k = f_{k+1} – f_k = f_k (r – 1) = f_k \left(2^{1/b} – 1\right) $$

です。途中で $f_{k+1} = f_k \cdot r$ という等比の関係を使い、共通因子 $f_k$ をくくり出しました。この帯域幅を採用すると、Q値は

$$ \begin{equation} Q = \frac{f_k}{\Delta f_k} = \frac{f_k}{f_k (2^{1/b} – 1)} = \frac{1}{2^{1/b} – 1} \end{equation} $$

となります。ここで $f_k$ が分子と分母でみごとに打ち消し合うことに注目してください。結果として $Q$ は $f_k$ に依存せず、分割数 $b$ だけで決まる定数になりました。これこそ「定Q（Constant-Q）」という名前の由来です。$b = 12$ なら $Q = 1 / (2^{1/12} – 1) \approx 16.8$ となります。

ステップ3：窓長の導出

いよいよ核心、各ビンの窓長 $N_k$ を求めます。出発点は、STFTのところで出てきた「窓長 $N$ の窓を使うと周波数分解能は $\Delta f = f_s / N$ になる」という基本関係です。CQTでもビン $k$ の窓は、そのビンが受け持つ帯域幅 $\Delta f_k$ 相当の分解能を実現できる長さでなければなりません。そこで $\Delta f_k = f_s / N_k$ と置くと、

$$ N_k = \frac{f_s}{\Delta f_k} $$

が得られます。ここに $\Delta f_k = f_k / Q$（Qの定義を変形したもの）を代入すると、

$$ \begin{equation} N_k = \frac{f_s}{f_k / Q} = Q \cdot \frac{f_s}{f_k} \end{equation} $$

という、CQTの中心となる関係式が導かれます。言葉にすると「ビンの窓長は中心周波数 $f_k$ に反比例する」ということです。低周波（$f_k$ 小）では窓が長く、高周波（$f_k$ 大）では窓が短くなる——まさに冒頭で欲しかった性質が、Q一定という1つの要請から自動的に出てきました。

この式の意味をもう少し噛み砕きます。$N_k = Q \cdot f_s / f_k$ をサンプル数ではなく「周期の個数」で見ると、窓に入る中心周波数の波の数は

$$ \frac{N_k}{f_s / f_k} = \frac{Q \cdot f_s / f_k}{f_s / f_k} = Q $$

と、これも $Q$ そのものになります。つまりCQTはどの周波数でも「ちょうど $Q$ 周期分の波」を窓に収めて解析するのです。低周波だろうと高周波だろうと、観察する「波の個数」は同じ。これが「相対的な分解能（=Q）が一定」ということの最も直感的な言い換えです。

各ビンの窓長が決まったので、次にこの可変長窓を使って実際に変換係数を計算する式を組み立てましょう。

左の両対数グラフでは、窓長 $N_k$ が中心周波数 $f_k$ に対して傾き $-1$ の直線に乗っており、$N_k \propto 1/f_k$ という反比例関係（式3）が正確に実装されていることを確認できます。右図は窓の時間長（ms）で表示したもので、A1（55 Hz）では約306ms、A4（440 Hz）で約38ms、A7（3520 Hz）で約5msと、低音ほど長い時間をかけて「じっくり」測定していることがわかります。導出した $N_k = Q f_s / f_k$ が実際に正しく機能しています。

定Q変換の定義式

各ビン $k$ について、中心周波数 $f_k$、窓長 $N_k = Q f_s / f_k$ が決まりました。これを使って、信号 $x[n]$ のビン $k$ における定Q変換係数 $X^{CQ}[k]$ を定義します。

STFTの各周波数成分は「複素正弦波 $e^{-j 2\pi f n / f_s}$ を窓でかけたものとの内積」でした。CQTもまったく同じ形ですが、周波数 $f_k$ ごとに窓長 $N_k$ が変わる点だけが異なります。1フレーム（ある時刻周辺の切り出し）に対する定Q変換係数は

$$ \begin{equation} X^{CQ}[k] = \frac{1}{N_k} \sum_{n=0}^{N_k – 1} w_k[n] \, x[n] \, e^{-j 2\pi Q \frac{n}{N_k}} \end{equation} $$

と定義されます。各記号の意味を確認しましょう。

指数部 $e^{-j 2\pi Q n / N_k}$ がCQTの肝です。STFTでは指数部が $e^{-j 2\pi k n / N}$ で、$N$ が固定でした。CQTでは窓長 $N_k$ の中にちょうど $Q$ 周期の複素正弦波が収まるように作られています。実際、$N_k = Q f_s / f_k$ を $e^{-j 2\pi Q n / N_k}$ に代入すると

$$ e^{-j 2\pi Q \frac{n}{N_k}} = e^{-j 2\pi Q \frac{n}{Q f_s / f_k}} = e^{-j 2\pi f_k \frac{n}{f_s}} $$

となり、これは確かに物理的な周波数 $f_k$ の複素正弦波を表しています。途中で $Q$ が約分されて消えるところがポイントです。つまり $X^{CQ}[k]$ は「信号の中に周波数 $f_k$ の成分がどれだけ含まれているか」を、$f_k$ に応じた長さの窓で測った量だと解釈できます。

正規化係数 $1/N_k$ について補足します。窓長が長いビン（低周波）ほど和の項数が多くなるので、そのまま足すと低周波の係数だけ大きく出てしまいます。$1/N_k$ で割ることで、振幅の異なる各ビンを公平に比較できるようになります。

上段は各周波数ビンの時間領域カーネル（窓つき複素正弦波の実部）で、左（A1: 55 Hz）では数百msにわたって波が広がっているのに対し、右（A5: 880 Hz）では数ms程度に収まっています。下段は各カーネルの振幅スペクトルで、それぞれの中心周波数（赤破線）にエネルギーが集中しています。この「周波数領域でのスパース性」こそ、高速CQTがFFT後の行列積を軽量化できる根拠になっています。

この定義式をそのまま素直に計算することもできますが、フレームごと・ビンごとに長さの違う和を毎回計算するのは効率が悪いです。次のセクションで、行列演算1回にまとめる「カーネル行列」の考え方を導入します。

カーネル行列による効率的な実装

定義式（4）をフレームごとに愚直に計算すると、各フレームで $K$ 個のビン、各ビンで $N_k$ 回の積和、というループが必要で計算量が大きくなります。これを見通しよく・高速にするのがカーネル行列（spectral kernel）による定式化です。Brown と Puckette（1992年）が提案した古典的な手法で、現代の高速CQTの土台にもなっています。

アイデアはこうです。定義式（4）は、信号フレーム $\bm{x}$（長さ $N$、最大窓長 $N_{\max} = N_0$ に合わせて取る）と、ビン $k$ に対応する「窓つき複素正弦波」 $\bm{t}_k$ との内積になっています。

$$ X^{CQ}[k] = \sum_{n=0}^{N-1} \overline{T_k[n]} \, x[n], \qquad T_k[n] = \frac{1}{N_k} w_k[n] \, e^{j 2\pi Q n / N_k} $$

ここで $T_k[n]$ をビン $k$ の時間領域カーネル（temporal kernel）と呼びます。$\overline{(\cdot)}$ は複素共役です。窓長 $N_k$ より外（$n \ge N_k$）では $T_k[n] = 0$ とゼロ詰めし、すべてのカーネルを長さ $N$ にそろえます。すると $K$ 本のカーネルを縦に積んだカーネル行列 $\bm{T} \in \mathbb{C}^{K \times N}$ を作れて、変換は1回の行列・ベクトル積

$$ \begin{equation} \bm{X}^{CQ} = \bm{T}^{*} \bm{x} \end{equation} $$

で書けます（$\bm{T}^{*}$ は共役を取った行列、$\bm{x}$ はフレームの列ベクトル）。フレームをずらしながらこの積を繰り返せば、CQTスペクトログラムが得られます。

ここまでなら「ループを行列積にまとめただけ」ですが、Brown-Puckette のもう一つの工夫は計算を周波数領域に移すことです。時間領域カーネル $T_k[n]$ をあらかじめ離散フーリエ変換してスペクトルカーネル $\hat{T}_k = \mathrm{DFT}(T_k)$ を作っておきます。パーセバルの関係（内積はフーリエ変換しても保たれる）より、内積は周波数領域でも計算できます。

$$ X^{CQ}[k] = \sum_{n} \overline{T_k[n]}\, x[n] = \frac{1}{N}\sum_{m} \overline{\hat{T}_k[m]}\, \hat{X}[m] $$

ここで $\hat{X} = \mathrm{DFT}(x)$ です。決定的に重要なのは、各 $\hat{T}_k$ が周波数領域では中心周波数 $f_k$ の周りにごく少数の非ゼロ要素しか持たない（スパースになる）ことです。時間領域では窓つき正弦波が窓全体に広がっているのに対し、周波数領域ではエネルギーが $f_k$ 付近に集中するためです。したがって $\hat{\bm{T}}$ をスパース行列として保持すれば、1回のFFTで $\hat{X}$ を求めたあと、スパースな行列積でCQT係数が一気に求まります。これがBrownとPucketteの高速CQTの要点です。

まとめると、効率実装は次の流れになります。

1. ビン配置 $f_k$、Q、窓長 $N_k$ を決める（信号に依存しない準備）
2. 各ビンの時間カーネル $T_k$ を作り、DFTしてスペクトルカーネル $\hat{T}_k$ を得る。小さい要素を切り捨ててスパース化する
3. これらを行列 $\hat{\bm{T}}$ にまとめる（ここまで前計算。1度だけ）
4. 入力フレームを1回FFTして $\hat{X}$ を得る
5. $\bm{X}^{CQ} = \hat{\bm{T}}^{*} \hat{X} / N$ で全ビンの係数を一括計算する

本記事のPython実装では、理解しやすさを優先して時間領域カーネル（式5）で素直に計算しますが、librosa など実用ライブラリは上記のスペクトルカーネル方式や、オクターブごとにダウンサンプリングして窓長を一定に保つ「再帰的ダウンサンプリング法」で高速化しています。

仕組みが揃ったので、いよいよ具体的な信号にSTFTとCQTを適用し、両者の見え方の違いを目で確かめましょう。

この図は時間周波数平面を「タイル」として表現したものです。左のSTFTは全周波数で同じサイズのタイルが均等に並ぶ（固定窓）。右のCQTは低周波（下段）ほど縦長（周波数分解能が高い）、高周波（上段）ほど横長（時間分解能が高い）のタイルが並んでいます。同じ時間周波数平面を覆う「タイルの総面積」は同じでも、周波数帯ごとに形が最適化されているのがCQTの特徴です。

具体例：対数掃引チャープでの振る舞い

理論を確認するのにうってつけの信号が対数掃引チャープ（logarithmic sweep）です。これは時間とともに周波数が指数関数的に上がっていく（オクターブが一定の速さで上がる）信号で、たとえば $t=0$ で 50 Hz、$t=T$ で 8000 Hz まで掃引します。

なぜこの信号が良いかというと、対数掃引チャープは対数周波数軸の上では「まっすぐな斜め線」として現れるはずだからです。CQTは周波数ビンが対数的に並んでいるので、CQTスペクトログラムを対数軸で描けば、チャープは綺麗な直線になり、しかも線の太さ（=その瞬間の周波数分解能）が低音から高音まで一定になることが期待されます。一方STFTは周波数ビンが線形なので、対数軸で描くと低周波側でビンが粗くスカスカに、高周波側で過剰に密になり、線の太さも周波数によって大きく変わるはずです。

定性的な予想を整理しておきます。

• 低周波域：CQTは長い窓を使うので周波数方向にシャープ（線が細い）。STFTは窓が短いまま（相対的に）なので、低周波では周波数分解能が足りずぼやける
• 高周波域：CQTは短い窓を使うので時間方向に俊敏（時間のにじみが小さい）。STFTは窓長固定なので、高周波では周波数分解能が過剰で、時間応答が相対的に鈍い

この予想が実際に成り立つかを、次のPython実装で確かめます。

Pythonでの実装

ここからは、(1) CQTのビン配置と窓長を確認する図、(2) スクラッチ実装したCQTとSTFTを対数掃引チャープに適用して比較する図、の順に作っていきます。まずは必要なライブラリと共通設定です。

import numpy as np
import matplotlib.pyplot as plt

# 共通設定
fs = 22050          # サンプリング周波数 [Hz]
f_min = 55.0        # CQT最低周波数 [Hz]（音名A1付近）
bins_per_octave = 12  # 1オクターブあたりのビン数（半音単位）

# Q値（分割数だけで決まる定数）
Q = 1.0 / (2 ** (1.0 / bins_per_octave) - 1.0)
print(f"Q = {Q:.3f}（窓に入る波の周期数）")

実行すると Q = 16.817 と表示されます。$b=12$ のとき、CQTはどの周波数でも約16.8周期分の波を窓に収めて解析している、ということです。この「波の周期数が一定」という事実が、これから見る性質すべての根っこにあります。

CQTのビン配置と窓長の確認

次に、式(1)と式(3)に従ってビンの中心周波数 $f_k$ と窓長 $N_k$ を計算し、その関係を可視化します。

# 周波数ビンの配置（対数=等比数列）
f_max = fs / 2 * 0.9        # ナイキスト手前まで
n_bins = int(np.ceil(bins_per_octave * np.log2(f_max / f_min)))
k = np.arange(n_bins)
f_k = f_min * 2 ** (k / bins_per_octave)   # 式(1)

# 各ビンの窓長（中心周波数に反比例） 式(3)
N_k = np.ceil(Q * fs / f_k).astype(int)

print(f"ビン数 K = {n_bins}")
print(f"最低ビン f_0 = {f_k[0]:.1f} Hz, 窓長 N_0 = {N_k[0]} サンプル")
print(f"最高ビン f_K = {f_k[-1]:.1f} Hz, 窓長 N_K = {N_k[-1]} サンプル")

# ビン中心周波数と窓長の関係を図示
fig, axes = plt.subplots(1, 2, figsize=(12, 4.5))

axes[0].plot(k, f_k, 'o-', color='royalblue', ms=3)
axes[0].set_xlabel('bin index k')
axes[0].set_ylabel('center frequency $f_k$ [Hz]')
axes[0].set_yscale('log')
axes[0].set_title('Log-spaced center frequencies')
axes[0].grid(True, which='both', alpha=0.3)

axes[1].plot(f_k, N_k, 'o-', color='crimson', ms=3)
axes[1].set_xlabel('center frequency $f_k$ [Hz]')
axes[1].set_ylabel('window length $N_k$ [samples]')
axes[1].set_xscale('log')
axes[1].set_yscale('log')
axes[1].set_title('Window length $N_k = Q f_s / f_k$')
axes[1].grid(True, which='both', alpha=0.3)

plt.tight_layout()
plt.savefig('cqt_bins.png', dpi=150, bbox_inches='tight')
plt.show()

左のグラフでは、ビン番号 $k$ に対して中心周波数 $f_k$ が片対数軸で直線に乗っています。これは $f_k = f_{\min} 2^{k/b}$ が指数関数、つまり対数軸で直線になることの確認で、ビンが等比数列（対数的に等間隔）に並んでいることを意味します。右のグラフでは、両対数軸で窓長 $N_k$ が中心周波数 $f_k$ に対して傾き $-1$ の直線になっています。これは $N_k = Q f_s / f_k$、すなわち $N_k \propto 1/f_k$ という反比例関係そのもので、低周波ほど窓が長く（数千サンプル）、高周波ほど窓が短く（数十サンプル）なることが読み取れます。導出した式(3)が正しく実装されていることが目で確認できました。

CQTのスクラッチ実装

カーネル（式5）を素直に組んでCQTを計算する関数を作ります。理解しやすさ優先で、各ビンごとに窓つき複素正弦波を生成し、フレームと内積を取ります。

def make_cqt_kernels(fs, f_min, bins_per_octave, f_max):
    """各ビンの時間領域カーネル（窓つき複素正弦波）を作る"""
    Q = 1.0 / (2 ** (1.0 / bins_per_octave) - 1.0)
    n_bins = int(np.ceil(bins_per_octave * np.log2(f_max / f_min)))
    f_k = f_min * 2 ** (np.arange(n_bins) / bins_per_octave)
    N_k = np.ceil(Q * fs / f_k).astype(int)
    N_max = N_k[0]  # 最長窓（=最低周波数）にそろえる

    kernels = np.zeros((n_bins, N_max), dtype=complex)
    for ki in range(n_bins):
        Nk = N_k[ki]
        n = np.arange(Nk)
        win = np.hanning(Nk)                       # ハニング窓
        # 物理周波数 f_k の複素正弦波（指数部に Q n / Nk）
        sig = np.exp(2j * np.pi * Q * n / Nk)
        # 窓中央にそろえて配置（時間的中心を合わせる）
        start = (N_max - Nk) // 2
        kernels[ki, start:start + Nk] = (win * sig) / Nk
    return kernels, f_k, N_max

このカーネル生成関数は、ビンごとに長さ $N_k$ のハニング窓つき複素正弦波を作り、最長窓 $N_{\max}$ の中央に配置してゼロ詰めしています。中央そろえにするのは、ビンによって窓長が違っても「同じ時刻を中心に解析している」と揃えるためです。(win * sig) / Nk の割り算が定義式の正規化 $1/N_k$ に対応します。

def cqt_spectrogram(x, kernels, N_max, hop):
    """フレームをずらしながらCQT係数を計算"""
    n_bins = kernels.shape[0]
    n_frames = 1 + (len(x) - N_max) // hop
    spec = np.zeros((n_bins, n_frames), dtype=complex)
    for fi in range(n_frames):
        frame = x[fi * hop: fi * hop + N_max]
        # 全ビン一括の内積（カーネルの共役と内積）
        spec[:, fi] = kernels.conj() @ frame
    return spec

このスペクトログラム関数は、長さ $N_{\max}$ のフレームを hop サンプルずつずらしながら、カーネル行列との内積（式5の $\bm{T}^{*}\bm{x}$）を1回の行列積で計算しています。kernels.conj() @ frame の一行が、全ビンの係数を同時に求める核心部分です。

対数掃引チャープへの適用とSTFTとの比較

いよいよ対数掃引チャープを作り、CQTとSTFTの両方を適用して見比べます。まずテスト信号を生成します。

# 対数掃引チャープの生成（50 Hz -> 8000 Hz）
dur = 4.0
t = np.arange(int(dur * fs)) / fs
f0, f1 = 50.0, 8000.0
# 瞬時周波数が指数関数的に増える信号の位相
K_sweep = dur / np.log(f1 / f0)
phase = 2 * np.pi * f0 * K_sweep * (np.exp(t / K_sweep) - 1)
x = np.sin(phase)

この信号は瞬時周波数が $f(t) = f_0 (f_1/f_0)^{t/T}$ で指数関数的に上がる、すなわちオクターブが一定速度で上昇するチャープです。対数周波数軸の上では、時間に対して直線的に周波数が上がる「斜め直線」として現れるはずです。

# CQTを計算
hop = 512
kernels, f_k, N_max = make_cqt_kernels(fs, f_min, bins_per_octave, fs/2*0.9)
cqt = cqt_spectrogram(x, kernels, N_max, hop)
cqt_db = 20 * np.log10(np.abs(cqt) + 1e-6)

# STFTを計算（固定窓長）
def stft(x, n_fft, hop):
    win = np.hanning(n_fft)
    n_frames = 1 + (len(x) - n_fft) // hop
    out = np.zeros((n_fft // 2 + 1, n_frames), dtype=complex)
    for fi in range(n_frames):
        frame = x[fi*hop: fi*hop+n_fft] * win
        out[:, fi] = np.fft.rfft(frame)
    return out

n_fft = 2048
S = stft(x, n_fft, hop)
S_db = 20 * np.log10(np.abs(S) + 1e-6)
f_stft = np.fft.rfftfreq(n_fft, 1/fs)

# 比較プロット（両方とも縦軸を対数周波数に）
fig, axes = plt.subplots(1, 2, figsize=(13, 5), sharey=False)

# STFT
t_stft = np.arange(S.shape[1]) * hop / fs
axes[0].pcolormesh(t_stft, f_stft, S_db, shading='auto', cmap='magma',
                   vmin=S_db.max()-60, vmax=S_db.max())
axes[0].set_yscale('log')
axes[0].set_ylim(f_min, fs/2*0.9)
axes[0].set_xlabel('time [s]'); axes[0].set_ylabel('frequency [Hz]')
axes[0].set_title(f'STFT (fixed N={n_fft})')

# CQT
t_cqt = np.arange(cqt.shape[1]) * hop / fs
axes[1].pcolormesh(t_cqt, f_k, cqt_db, shading='auto', cmap='magma',
                   vmin=cqt_db.max()-60, vmax=cqt_db.max())
axes[1].set_yscale('log')
axes[1].set_ylim(f_min, fs/2*0.9)
axes[1].set_xlabel('time [s]'); axes[1].set_ylabel('frequency [Hz]')
axes[1].set_title('Constant-Q Transform')

plt.tight_layout()
plt.savefig('cqt_vs_stft.png', dpi=150, bbox_inches='tight')
plt.show()

2つのスペクトログラムから、理論で予想した違いがはっきり読み取れます。第一に、CQT（右）では掃引チャープがほぼ一定の太さの綺麗な斜め直線として現れます。対数周波数軸の上で、低音から高音まで線の太さ（その瞬間の周波数分解能の幅）が一定なのは、Q値が一定＝相対的な分解能が一定だからです。第二に、STFT（左）では同じチャープが、低周波側では太く（周波数方向にぼやけ）、高周波側では細いという不均一な線になります。これは窓長が固定で、絶対分解能 $\Delta f = f_s/N$ が一定のため、低周波では相対的に分解能が粗く、高周波では過剰に細かいことの表れです。第三に、STFTを対数軸で描くと低周波域でビンがスカスカに離れて見えるのに対し、CQTは全域でビンが対数的に均等に並んでおり、音楽的な「半音の格子」としてそのまま読めることがわかります。

低周波の音程分解能の比較

最後に、CQTが低周波で音程をよく分離できることを、近接する2つの低音トーンで確かめます。

# 近接する低音2トーン（A1=55Hz と A#1≒58.3Hz）を合成
f_tone1, f_tone2 = 55.0, 58.27   # 半音差
tt = np.arange(int(2.0 * fs)) / fs
x2 = np.sin(2*np.pi*f_tone1*tt) + np.sin(2*np.pi*f_tone2*tt)

cqt2 = np.abs(cqt_spectrogram(x2, kernels, N_max, hop)).mean(axis=1)
S2 = np.abs(stft(x2, n_fft, hop)).mean(axis=1)

fig, axes = plt.subplots(1, 2, figsize=(12, 4))
axes[0].plot(f_stft, S2, color='gray'); axes[0].set_xlim(40, 80)
axes[0].set_title(f'STFT spectrum (N={n_fft})')
axes[0].set_xlabel('frequency [Hz]'); axes[0].axvline(55, ls='--', c='b'); axes[0].axvline(58.27, ls='--', c='r')

axes[1].plot(f_k, cqt2, 'o-', color='crimson', ms=3); axes[1].set_xlim(40, 80)
axes[1].set_title('CQT spectrum')
axes[1].set_xlabel('frequency [Hz]'); axes[1].axvline(55, ls='--', c='b'); axes[1].axvline(58.27, ls='--', c='r')

plt.tight_layout()
plt.savefig('cqt_lowfreq_resolution.png', dpi=150, bbox_inches='tight')
plt.show()

このグラフでは、55 Hz と 58.27 Hz（半音差、約3.3 Hzしか離れていない）の2つの低音トーンに対するスペクトルを比べています。STFT（左、$N=2048$ では $\Delta f \approx 10.8$ Hz）では2本のピークが分離できず1つの山に潰れてしまうのに対し、CQT（右）では2つのピークがはっきり2本に分かれて見えます。これは低周波域でCQTが長い窓（高い周波数分解能）を使っているおかげで、まさに導出した $N_k = Q f_s / f_k$ が低周波で大きな値を取ることの実用上の御利益です。STFTで同じ分解能を得ようとすれば窓長を大幅に伸ばす必要があり、その分だけ高周波の時間応答を犠牲にしてしまいます。CQTは周波数ごとに窓長を変えることで、この二律背反を周波数帯ごとに最適化しているのです。

和音のCQTスペクトログラム — 音楽構造の可視化

最後に、CQTが和音（コード）の分析においていかに自然な表示になるかを示します。Aメジャーコード（A4: 440 Hz、C#5: 554 Hz、E5: 659 Hz とその倍音）を合成してCQTスペクトログラムを描くと、倍音の構造が対数軸上で等間隔に並ぶ音楽的なパターンが浮かび上がります。

スペクトログラムの横軸（シアン破線）は各音名の周波数に対応しています。A4（440 Hz）、A5（880 Hz）が対数軸上でちょうど同じ間隔に並んでいることや、倍音列がほぼ等間隔の縦縞として現れることが視覚的に確認できます。また、エンベロープを付けているため時間が進むにつれてスペクトルが薄くなる様子も見えており、音の減衰が自然に表現されています。STFTで同じ解析をすると、低音域のA4とA5の間隔が対数軸で密になり、倍音の「等間隔性」が視覚的に崩れてしまいます。

高周波バースト信号の時間分解能

CQTのもう一つの利点は、高周波では短い窓を使うため時間応答が鋭いことです。2000 Hz のバースト信号（50 msパルス）を使って、STFTとCQTの時間分解能を比べてみます。

STFT（左、窓長 4096 サンプル ≈ 186 ms）ではバースト信号のエネルギーが時間方向に大きく広がっており、50 ms のパルスが100 ms以上にわたって「にじんで」見えます。これは窓長が信号の持続時間より長いためです。一方、CQT（右、2000 Hz ビンの窓長 ≈ 8 ms）ではバースト区間（赤いゾーン）のエネルギーが鋭く立ち上がり、パルスの立ち上がりと終端が明確に識別できます。このように「周波数帯ごとに窓長を最適化する」というCQTの設計が、高周波における時間応答の向上として現れています。

以上の実験で、理論で導いた性質（低周波: 長窓で周波数精度↑、高周波: 短窓で時間精度↑）がすべて実際の信号で確認できました。

まとめ

本記事では、定Q変換（CQT）の理論・導出・Python実装を解説しました。

• 動機：STFTは全周波数で窓長が固定なため $\Delta f = f_s/N$ も一定で、Q $= f/\Delta f$ が周波数に比例してしまう。音楽や聴覚が感じるのは周波数の「差」ではなく「比」なので、相対的な分解能（=Q）を一定にしたい
• ビン配置：Q一定の帰結として、ビンは等比数列 $f_k = f_{\min} 2^{k/b}$ に対数的に並ぶ
• Q値：隣接ビンに接する帯域幅を採るとQは $f_k$ に依存せず $Q = 1/(2^{1/b}-1)$ という定数になる（定Qの名の由来）
• 窓長の導出：$\Delta f_k = f_s/N_k$ と $\Delta f_k = f_k/Q$ から $N_k = Q f_s / f_k$。窓長は中心周波数に反比例し、どの周波数でもちょうど $Q$ 周期分の波を窓に収める
• 効率実装：時間カーネルを束ねたカーネル行列、さらに周波数領域でスパースなスペクトルカーネルにすると行列積1回で全ビンを計算できる
• Python実装：対数掃引チャープがCQTでは一定の太さの直線、STFTでは不均一な線になること、低周波の半音差をCQTは分離できるがSTFTは潰してしまうことを確認

CQTは「分解能」「窓」「周波数」という時間周波数解析の三角関係を、対数スケールという視点で再設計したものです。同じ問題意識をより数学的に追究すると、解像度を連続的に変えるウェーブレット変換（CQTはウェーブレット変換の一種とみなせます）や、メル尺度に基づくメルスペクトログラムとMFCCへとつながります。

次のステップとして、以下の記事も参考にしてください。

• 短時間フーリエ変換（STFT）とスペクトログラム — CQTの出発点となる固定分解能の時間周波数解析
• 窓関数（ハニング窓・ハミング窓）の役割 — 各ビンの窓と周波数分解能の関係
• ウェーブレット変換の理論と実装 — 解像度を連続的に変える多重解像度解析

窓関数（ハニング・ハミング・ブラックマン）の定義と効果

ハニング窓・ハミング窓・ブラックマン窓の数式・周波数特性・スペクトルリーケージへの効果をPythonで比較。CQTでも各ビンに窓をかけるため必須の前提知識