Chirp-Z変換（CZT）の理論と導出と実装

電力系統の周波数が「ちょうど 50 Hz か、それとも 50.02 Hz か」を見分けたいとき、あるいはレーダーで「2つの目標が 0.1 Hz だけ離れたドップラー周波数を持つ」ことを確かめたいとき、私たちは FFT のスペクトルを拡大して眺めたくなります。しかし通常の FFT で得られる周波数分解能は $\Delta f = f_s / N$ で固定されており、ある狭い帯域だけをもっと細かく見たいと思っても、$N$ を増やして全帯域を計算し直すしかありません。サンプリング周波数 $f_s = 1\,\text{MHz}$ のデータで $0.01\,\text{Hz}$ の分解能が欲しければ、$N = 10^8$ 点もの FFT が必要になってしまいます。

この「見たい帯域だけを虫眼鏡で拡大する」ことを可能にするのが Chirp-Z 変換（Chirp-Z Transform, CZT） です。CZT は離散フーリエ変換（DFT）を一般化したもので、$z$ 平面上の単位円に限らず、任意の開始点・任意の角度間隔・螺旋を描く経路の上で信号のスペクトルを評価できます。これを使うと、興味のある周波数帯域だけを高い分解能でサンプリングする「ズーム FFT」が、全帯域の FFT を計算することなく実現できます。

CZT は応用範囲が広く、たとえば次のような場面で活躍します。

• 狭帯域スペクトル解析: 電力系統の周波数監視、振動計測における特定ピークの精密測定、近接した2本のスペクトル線の分離
• レーダー・ソナーのドップラー解析: 微小な速度差を持つ複数目標の分離、ターゲットの精密な周波数推定
• 任意長 DFT の高速計算: $N$ が大きな素数でも、Bluestein のアルゴリズムによって $\mathcal{O}(N \log N)$ で DFT を計算できる

本記事の内容

• DFT を「$z$ 平面の単位円上の等間隔サンプリング」として捉え直す
• CZT が $z_k = A W^{-k}$ という螺旋経路を取れることと、その定義の導出
• Bluestein 分解（$nk = \tfrac{1}{2}\{k^2 + n^2 – (k-n)^2\}$）により CZT が畳み込みに帰着し、FFT 3 回で計算できることの導出
• Python で狭帯域信号にズーム FFT（CZT）を適用し、通常 FFT より高い周波数分解能で近接スペクトル線を分離する比較実験

左の DFT は単位円を $N$ 等分した固定格子しか使えません。右の CZT は任意の開始点 $A$ から始まり角度間隔 $\phi_0$ で進む螺旋状の点列を自由に設定でき、単位円の弧の一部（ズーム FFT）や内向きの螺旋など多様な経路を取れます。この「どこでも止まれる自由さ」が CZT の最大の特徴です。

前提知識

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

• 離散フーリエ変換（DFT）とは — CZT は DFT の一般化であり、DFT の定義と性質が出発点になります
• FFT（高速フーリエ変換）のアルゴリズム — CZT の高速化は内部で FFT を 3 回呼び出すことで実現されます

CZT は「DFT を拡張し、FFT で速く計算する」という二段構えの手法です。したがって DFT の定義（時間領域の数列を周波数成分に分解する操作）と、FFT が $\mathcal{O}(N \log N)$ で DFT を計算する仕組み、特に循環畳み込みが FFT の積で計算できるという性質を押さえておくと、本記事の導出が一段とクリアになります。

DFT を $z$ 平面の単位円として捉え直す

CZT を理解する近道は、まず私たちがよく知っている DFT を「少し違う角度」から眺めることです。長さ $N$ の複素数列 $x[0], x[1], \dots, x[N-1]$ に対する DFT は、次のように定義されます。

$$ \begin{equation} X[k] = \sum_{n=0}^{N-1} x[n]\, e^{-j \frac{2\pi}{N} nk}, \quad k = 0, 1, \dots, N-1 \end{equation} $$

ここで指数の中身に注目しましょう。$e^{-j \frac{2\pi}{N} nk}$ は、$\left(e^{-j \frac{2\pi}{N} k}\right)^n$ と書き直せます。つまり DFT は、各 $k$ について複素数 $z_k = e^{j \frac{2\pi}{N} k}$ を用意し、信号の $z$ 変換

$$ \begin{equation} X(z) = \sum_{n=0}^{N-1} x[n]\, z^{-n} \end{equation} $$

を $z = z_k$ で評価しているのと同じことなのです。

ここで $z_k = e^{j \frac{2\pi}{N} k}$ は、複素平面の原点を中心とする半径 1 の円（単位円）の上に並ぶ点です。$k = 0$ のときは角度 $0$（実軸上の点 $z = 1$）、$k = 1$ のときは角度 $2\pi/N$、というように、単位円を $N$ 等分した点に対応します。つまり DFT とは、信号の $z$ 変換を「単位円上の等間隔な $N$ 点」で評価する操作だと言い換えられるのです。

この見方をすると、DFT の制約がはっきり見えてきます。

• 評価点は必ず単位円上にある（半径は 1 に固定）
• 評価点は必ず角度 0 から始まる（最初の点が $z = 1$ に固定）
• 評価点の角度間隔は必ず $2\pi / N$（全周 $2\pi$ を $N$ 等分、つまり全帯域を均等に覆う）

「特定の狭い周波数帯域だけを細かく見たい」という私たちの要求は、この最後の制約に阻まれます。DFT では角度間隔が $2\pi/N$ に固定されているので、$0$ から $2\pi$ までの全帯域を一様にしか刻めません。一部分だけを拡大することができないのです。

そこで自然な発想が生まれます。評価点を単位円上の等間隔に縛らず、もっと自由に置けないだろうか？ たとえば「角度 $\theta_0$ から始めて、$\Delta\theta$ ずつ刻んで $M$ 個」と決められれば、興味のある帯域だけを好きな分解能でサンプリングできます。この一般化こそが CZT です。次節では、この自由な評価点をどう数式で表すかを見ていきましょう。

上段の FFT は全帯域 $[0, f_s)$ を均等に刻むため、$[190, 210]$ Hz の着目帯域内には数点しか格子が入りません。下段の CZT では同じ帯域を 50 点で細かくカバーしており、分解能が大幅に向上しています。FFT の格子が粗い理由と、CZT がそれを解決する仕組みをここで直感的に把握してから、数式の導出に進みましょう。

Chirp-Z変換とは — 螺旋経路のサンプリング

DFT の評価点 $z_k = e^{j \frac{2\pi}{N} k}$ を一般化することを考えます。私たちが欲しいのは、「ある開始点から始めて、一定の比率で進んでいく」ような点列です。等差数列ではなく等比数列にするのがポイントです。なぜなら、$z^{-n}$ という $z$ のべき乗を扱う $z$ 変換と相性が良いからです。

そこで、CZT では評価点を次のように定めます。

$$ \begin{equation} z_k = A W^{-k}, \quad k = 0, 1, \dots, M-1 \end{equation} $$

ここで $A$ と $W$ は複素数のパラメータで、$M$ は欲しい出力点数です。$A$ は「最初の点 $z_0 = A$」を、$W$ は「隣り合う点の比 $z_{k+1}/z_k = W^{-1}$」を決めます。一般に複素数として

$$ A = A_0\, e^{j\theta_0}, \qquad W = W_0\, e^{j\phi_0} $$

と書けます。$A_0, W_0$ は正の実数（絶対値）、$\theta_0, \phi_0$ は偏角です。これを代入すると評価点は

$$ z_k = A_0 W_0^{-k}\, e^{j(\theta_0 + k\phi_0)} $$

となります。この式の意味をかみ砕きましょう。

• 絶対値 $|z_k| = A_0 W_0^{-k}$: $k$ が増えるごとに半径が $W_0^{-1}$ 倍になります。$W_0 = 1$ なら半径一定（円弧）、$W_0 \neq 1$ なら半径が指数的に変化する螺旋になります。
• 偏角 $\arg z_k = \theta_0 + k\phi_0$: 開始角 $\theta_0$ から $\phi_0$ ずつ等間隔に回転します。$\theta_0$ で「どこから見始めるか」、$\phi_0$ で「どれくらい細かく刻むか」を制御できます。

このように、CZT の評価点は一般に $z$ 平面上で螺旋を描く点列になります。「Chirp（チャープ）」という名前は、後で見るように変換の核に現れる $e^{j \alpha n^2}$ という二次位相の項が、時間とともに周波数が変化する「チャープ信号」と同じ形をしていることに由来します。

$A = e^{j\theta_0}$ が評価点の開始角（図の赤い弧の部分）を、$W^{-1} = e^{j\phi_0}$ が隣り合う点の角度間隔（緑の弧）を決めます。$|A_0| = |W_0| = 1$ のとき評価点は単位円弧の上に並び、これがズーム FFT の基本設定です。$W_0 \neq 1$ にすると内向き・外向きの螺旋も作れます。

CZT は、この評価点を使った $z$ 変換として定義されます。

$$ \begin{equation} X[k] = \sum_{n=0}^{N-1} x[n]\, z_k^{-n} = \sum_{n=0}^{N-1} x[n]\, A^{-n} W^{nk}, \quad k = 0, 1, \dots, M-1 \end{equation} $$

入力長 $N$ と出力長 $M$ が独立に選べることに注意してください。DFT では入力と出力が同じ $N$ 点でしたが、CZT では「$N$ 点の信号から、好きな螺旋上の $M$ 点」を計算できます。

DFT が CZT の特別な場合であることも確認しておきましょう。$A = 1$（開始点 $z_0 = 1$）、$W = e^{-j \frac{2\pi}{N}}$（角度間隔 $2\pi/N$、半径一定）、$M = N$ とおくと、$z_k = e^{j \frac{2\pi}{N} k}$ となり、CZT の定義式は DFT の定義式に一致します。つまり CZT は DFT を真に一般化したものなのです。

ここで重要な疑問が湧きます。CZT の定義式を素朴に計算すると、$M \times N$ 回の積和が必要で計算量は $\mathcal{O}(NM)$ です。これでは大きな $N, M$ に対して遅すぎます。FFT のような高速化はできないのでしょうか？ 実はできます。次節では、指数の中の $nk$ という積を巧妙に分解することで、CZT が畳み込みに書き換えられ、FFT 3 回で計算できることを導きます。

Bluestein の恒等式は「$n$ と $k$ が絡んだ積 $nk$」を「$k$ だけ / $n$ だけ / 差 $(k-n)$ だけ」の 3 項に分離します。差 $(k-n)$ が現れたことで畳み込みの形が生まれ、FFT による高速化への扉が開きます。この図の流れが Bluestein 分解の骨子です。

ズーム解析のためのパラメータ設計

高速化の導出に入る前に、CZT で「ズーム FFT」を実現するための具体的なパラメータ設計を押さえておきましょう。これが CZT を使う最大の動機なので、ここで直感を固めておくと後の議論が腑に落ちます。

サンプリング周波数 $f_s$ で得た信号のスペクトルを、周波数 $f_1$ から $f_2$ までの帯域に限って $M$ 点で評価したいとします。デジタル角周波数（正規化角周波数）で考えると、周波数 $f$ は $z$ 平面上の角度 $2\pi f / f_s$ に対応します。したがって開始角と角度間隔は

$$ \theta_0 = \frac{2\pi f_1}{f_s}, \qquad \phi_0 = -\frac{2\pi (f_2 – f_1)}{f_s\, M} $$

と決めればよいことになります（$\phi_0$ に負号が付くのは $W = W_0 e^{j\phi_0}$ で $z_k = A W^{-k}$ と定義したため、$k$ が増えると角度が $-\phi_0$ ずつ増える、すなわち周波数が増える向きに進むよう符号を合わせています）。半径は単位円上をなぞりたいので $A_0 = 1, W_0 = 1$ とします。すると

$$ A = e^{j \frac{2\pi f_1}{f_s}}, \qquad W = e^{-j \frac{2\pi (f_2 – f_1)}{f_s\, M}} $$

となります。このとき CZT が評価する周波数は $f_1, f_1 + \delta, f_1 + 2\delta, \dots$（$\delta = (f_2 – f_1)/M$）と、$[f_1, f_2)$ の帯域を $M$ 等分した格子になります。

ここで注目すべきは、得られる分解能 $\delta = (f_2 – f_1)/M$ が、サンプリング周波数 $f_s$ にも信号長 $N$ にも縛られないことです。$f_2 – f_1$ を小さく取れば取るほど、同じ $M$ でも分解能はいくらでも細かくなります。通常 FFT の分解能 $f_s/N$ と比べて、ズーム倍率 $f_s / (f_2 – f_1)$ だけ細かく見られるわけです。

この「ズームできる」という性質こそが CZT の真価です。では、この計算をどうやって FFT 3 回に落とし込むのか。次節でその数学的な仕掛けを丁寧に追っていきます。

Bluestein のアルゴリズム — 畳み込みへの分解

CZT の定義式

$$ X[k] = \sum_{n=0}^{N-1} x[n]\, A^{-n} W^{nk} $$

を高速化する鍵は、指数の肩にある積 $nk$ をどう扱うかにあります。$W^{nk}$ の $nk$ がそのままだと、$n$ と $k$ が絡み合っていて分離できません。FFT で計算したい畳み込みの形（$\sum_n f[n] g[k-n]$ のように $k-n$ という差で書ける形）に持ち込めないのです。

ここで Bluestein が用いた巧妙な恒等式が登場します。それは、積 $nk$ を二乗の差で表すというものです。次の式を出発点にします。

$$ \begin{equation} nk = \frac{1}{2}\left\{ k^2 + n^2 – (k-n)^2 \right\} \end{equation} $$

この恒等式が正しいことは、右辺を展開すれば確かめられます。$(k-n)^2 = k^2 – 2nk + n^2$ を代入すると

$$ \frac{1}{2}\left\{ k^2 + n^2 – (k^2 – 2nk + n^2) \right\} = \frac{1}{2}\left\{ 2nk \right\} = nk $$

となり、確かに左辺に一致します。素朴な式変形ですが、これが CZT 高速化の心臓部です。なぜなら、$nk$ という「絡んだ積」が、$k$ だけの項 $k^2$、$n$ だけの項 $n^2$、そして差 $(k-n)$ だけの項 $(k-n)^2$ に分離されたからです。差 $(k-n)$ が現れたことが決定的で、これが畳み込みの形を生みます。

この恒等式を $W^{nk}$ に代入してみましょう。指数法則 $W^{a+b} = W^a W^b$ を使うと

$$ W^{nk} = W^{\frac{1}{2}\left\{ k^2 + n^2 – (k-n)^2 \right\}} = W^{\frac{k^2}{2}}\, W^{\frac{n^2}{2}}\, W^{-\frac{(k-n)^2}{2}} $$

と 3 つの因子に分かれます。これを CZT の定義式に代入します。

$$ X[k] = \sum_{n=0}^{N-1} x[n]\, A^{-n}\, W^{\frac{k^2}{2}}\, W^{\frac{n^2}{2}}\, W^{-\frac{(k-n)^2}{2}} $$

ここで $W^{k^2/2}$ は和の添字 $n$ に依存しないので、和の外に出せます。

$$ X[k] = W^{\frac{k^2}{2}} \sum_{n=0}^{N-1} \underbrace{x[n]\, A^{-n}\, W^{\frac{n^2}{2}}}_{g[n]} \; \underbrace{W^{-\frac{(k-n)^2}{2}}}_{h[k-n]} $$

ここで2つの数列を定義しました。

$$ g[n] = x[n]\, A^{-n}\, W^{\frac{n^2}{2}}, \qquad h[m] = W^{-\frac{m^2}{2}} $$

すると CZT は

$$ \begin{equation} X[k] = W^{\frac{k^2}{2}} \sum_{n=0}^{N-1} g[n]\, h[k-n] \end{equation} $$

と書けます。和の部分 $\sum_n g[n]\, h[k-n]$ は、まさに数列 $g$ と $h$ の畳み込みです。つまり、CZT は次の 3 ステップに分解できました。

1. 入力 $x[n]$ にチャープ的な重み $A^{-n} W^{n^2/2}$ をかけて $g[n]$ を作る（プリマルチプライ）
2. $g$ と $h[m] = W^{-m^2/2}$ を畳み込む
3. 畳み込み結果に再びチャープ重み $W^{k^2/2}$ をかける（ポストマルチプライ）

数列 $h[m] = W^{-m^2/2}$ に現れる $W^{-m^2/2} = e^{j \pi m^2 \cdot (\text{定数})}$ という二次位相の項は、時間 $m$ とともに位相が二次関数的に変化します。瞬時周波数（位相の微分）が時間に比例して増えるこの信号こそが、レーダーでおなじみの チャープ信号です。CZT の名前「Chirp-Z」はここから来ています。

チャープ列 $h[m]$ の実部と虚部はともに、$m = 0$ 付近では振動が緩やか（低「周波数」）で、端に行くほど振動が密になっています（高「周波数」）。この「時間とともに周波数が変化する」二次位相の信号がチャープ（啼き声のように周波数が滑らかに変わる）と呼ばれる由来です。CZT は「チャープ列との畳み込み」という操作の中にこのチャープ信号を埋め込んでいます。

残る問題は、ステップ 2 の畳み込みをどう速く計算するかです。畳み込みは FFT で高速化できる、というのが FFT 理論の最重要結果でした。次節でその具体的な手順を組み立てます。

畳み込みを FFT 3 回で計算する

畳み込み $\sum_n g[n] h[k-n]$ を FFT で計算するには、注意が必要です。FFT が直接高速化できるのは循環畳み込み（インデックスが周期的に巻き戻る畳み込み）であって、CZT に現れるのは線形畳み込み（添字 $k-n$ が負にもなりうる、巻き戻りのない畳み込み）だからです。両者を一致させるには、ゼロ詰め（zero-padding）によって配列を十分長くする工夫が要ります。

CZT で必要な出力は $k = 0, 1, \dots, M-1$、和を取る $n$ は $0, 1, \dots, N-1$ です。したがって $h[k-n]$ の添字 $m = k – n$ は、最小で $k=0, n=N-1$ のとき $m = -(N-1)$、最大で $k=M-1, n=0$ のとき $m = M-1$ の範囲を動きます。すなわち $h$ は $m = -(N-1)$ から $m = M-1$ までの $L = M + N – 1$ 点が必要です。

そこで、線形畳み込みの長さ $L = N + M – 1$ 以上の 2 のべき乗 $L_{\text{fft}}$ を選び、循環畳み込みが線形畳み込みと一致するように両配列をゼロ詰めします。具体的な計算手順は次の通りです。

1. プリマルチプライ: $g[n] = x[n]\, A^{-n} W^{n^2/2}$（$n = 0, \dots, N-1$）を計算し、長さ $L_{\text{fft}}$ までゼロ詰めする。
2. チャープ列の構成: $h[m] = W^{-m^2/2}$ を $m = -(N-1), \dots, M-1$ について用意し、循環畳み込みで正しい位置に来るよう配列を構成する（負の添字は末尾に巻き付ける）。
3. FFT その1・その2: $g$ と $h$ をそれぞれ FFT する。$\;G = \text{FFT}(g), \; H = \text{FFT}(h)$。
4. 周波数領域での積: $Y = G \cdot H$（要素ごとの積）。これが畳み込み定理「時間領域の畳み込み = 周波数領域の積」に対応する。
5. 逆 FFT（FFT その3）: $y = \text{IFFT}(Y)$ を計算する。$y[k]$ が畳み込み $\sum_n g[n] h[k-n]$ を与える。
6. ポストマルチプライ: $X[k] = W^{k^2/2}\, y[k]$（$k = 0, \dots, M-1$）を計算する。

このうち FFT/IFFT が 3 回（順方向 2 回 $G, H$ と逆方向 1 回 $y$）登場します。実際には $h$ はデータ $x$ に依存しないため、$H = \text{FFT}(h)$ を事前計算しておけば、毎回の変換では FFT 2 回で済みます。

計算量を見積もりましょう。$L_{\text{fft}} = \mathcal{O}(N + M)$ なので、各 FFT は $\mathcal{O}((N+M)\log(N+M))$ です。したがって CZT 全体の計算量は

$$ \mathcal{O}\big((N + M)\log(N + M)\big) $$

となります。定義式を素朴に計算したときの $\mathcal{O}(NM)$ と比べると、劇的な高速化です。たとえば $N = M = 10^6$ なら、$\mathcal{O}(NM) = 10^{12}$ に対して $\mathcal{O}((N+M)\log(N+M)) \approx 2 \times 10^6 \times 21 \approx 4 \times 10^7$ と、5 桁近く速くなります。

このアルゴリズムには、もう一つ嬉しい副産物があります。$N$ が大きな素数であっても、CZT を $A=1, W=e^{-j2\pi/N}, M=N$ と設定すれば DFT になり、しかも内部のゼロ詰め長 $L_{\text{fft}}$ は 2 のべき乗に取れるので、任意長の DFT を $\mathcal{O}(N \log N)$ で計算できます。素数長で Cooley–Tukey 型 FFT が使えない場合の救済策として、これは Bluestein のアルゴリズムの名で広く使われています。

CZT の計算は「プリマルチプライ → FFT 2 回と IFFT 1 回 → ポストマルチプライ」という 6 ステップに整理されます。チャープ列 $h$ の FFT $H$ は $x$ に依存しないため事前計算が可能で、繰り返し計算する場面では実質 FFT 2 回で済みます。この整理されたフローが、CZT を実際のシステムに組み込む際の実装指針になります。

ここまでで CZT の理論と高速化の全体像が揃いました。次は具体的な数値で挙動を確認し、その後 Python でズーム FFT を実装して、通常 FFT では分離できない近接スペクトル線を CZT が分離できることを目で見て確かめましょう。

具体例 — ズーム倍率を数値で確認する

抽象論だけでは実感が湧きにくいので、具体的な数値でパラメータと分解能の関係を確認しましょう。

サンプリング周波数 $f_s = 1000\,\text{Hz}$、信号長 $N = 1000$ 点の信号があるとします。通常の FFT で得られる周波数分解能は

$$ \Delta f_{\text{FFT}} = \frac{f_s}{N} = \frac{1000}{1000} = 1\,\text{Hz} $$

です。つまり FFT のスペクトルは $0, 1, 2, \dots, 999\,\text{Hz}$ という 1 Hz 刻みの点でしか値を持ちません。もし信号に $200.0\,\text{Hz}$ と $200.3\,\text{Hz}$ の 2 本の正弦波が含まれていたら、両者は同じビン（$200\,\text{Hz}$ の点）の近くに埋もれてしまい、2 本あることすら判別できません。

そこで CZT を使い、$f_1 = 198\,\text{Hz}$ から $f_2 = 202\,\text{Hz}$ までの幅 $4\,\text{Hz}$ の帯域だけを $M = 800$ 点で評価することにします。このときの分解能は

$$ \delta = \frac{f_2 – f_1}{M} = \frac{202 – 198}{800} = \frac{4}{800} = 0.005\,\text{Hz} $$

です。FFT の $1\,\text{Hz}$ に対して $0.005\,\text{Hz}$ ですから、ズーム倍率は $1 / 0.005 = 200$ 倍にもなります。$200.0\,\text{Hz}$ と $200.3\,\text{Hz}$ の差 $0.3\,\text{Hz}$ は、この分解能なら $0.3 / 0.005 = 60$ 点ぶんも離れて見えるので、2 本のピークがはっきり分離されます。

CZT のパラメータは前節の設計式に従って

$$ A = e^{j \frac{2\pi \cdot 198}{1000}}, \qquad W = e^{-j \frac{2\pi \cdot 4}{1000 \cdot 800}} $$

と設定します。$|A| = 1, |W| = 1$ なので評価点は単位円上にあり、円弧の一部（角度にして $2\pi \cdot 198/1000$ から $2\pi \cdot 202/1000$ まで）を 800 点で細かくなぞることになります。

注意したいのは、ズームしてももとの信号の情報量以上のものは生まれないことです。CZT のズームは「FFT のビンの間を補間して滑らかに見せる」のではなく、信号の $z$ 変換そのものを細かい格子で正確に評価しています。とはいえ 2 本の正弦波を分離できるかどうかは、最終的には信号長 $N$ で決まるレイリー限界 $\approx f_s/N$ に支配されます。CZT はこの限界に張り付いた真のスペクトルを高分解能で「見せてくれる」道具であり、観測時間が短すぎて原理的に分離不能なものを魔法のように分けるわけではない、という点は心に留めておきましょう。それでも、後の実験で見るように、近接した（しかし原理的には分離可能な）スペクトル線を浮かび上がらせる力は絶大です。

左図は帯域幅を縦軸（対数スケール）の CZT 分解能で示したものです。帯域幅を 1/10 にすると分解能も 1/10（右図のズーム倍率は 10 倍）になることが一目でわかります。帯域幅 10 Hz では 100 倍、帯域幅 1 Hz では 1000 倍という指数的な改善が得られ、「どれだけ狭い帯域に集中するか」でズームの深さが決まります。

それでは、この具体例を実際に Python で計算してみましょう。

Pythonでの実装

CZT のスクラッチ実装

まずは Bluestein のアルゴリズムに基づいて CZT をゼロから実装します。前節の手順 1〜6 をそのままコードに落とし込みます。

import numpy as np

def czt(x, M, A, W):
    """
    Chirp-Z変換（Bluestein のアルゴリズム）
    x : 入力信号（長さ N）
    M : 出力点数
    A : 開始点 z0 = A（複素数）
    W : 評価点の比 z_k = A * W^{-k}（複素数）
    戻り値: X[k] = sum_n x[n] * (A*W^{-k})^{-n}, k=0,...,M-1
    """
    x = np.asarray(x, dtype=complex)
    N = len(x)

    # 線形畳み込みに必要な長さ以上の 2 のべき乗を選ぶ
    L = 1
    while L < N + M - 1:
        L *= 2

    # 二次位相（チャープ）列の指数 n^2/2, m^2/2 を用意
    n = np.arange(N)
    # プリマルチプライ: g[n] = x[n] * A^{-n} * W^{n^2/2}
    g = x * (A ** (-n)) * (W ** (n * n / 2.0))
    g = np.concatenate([g, np.zeros(L - N, dtype=complex)])

    # チャープ列 h[m] = W^{-m^2/2} を m = -(N-1)..(M-1) で構成し循環配置
    h = np.zeros(L, dtype=complex)
    m = np.arange(M)
    h[:M] = W ** (-(m * m) / 2.0)          # m = 0 .. M-1
    m_neg = np.arange(1, N)
    h[L - (N - 1):] = W ** (-(m_neg * m_neg) / 2.0)  # m = -(N-1) .. -1

    # FFT 3 回（順方向 2 回 + 逆方向 1 回）で循環畳み込み
    G = np.fft.fft(g)
    H = np.fft.fft(h)
    y = np.fft.ifft(G * H)

    # ポストマルチプライ: X[k] = W^{k^2/2} * y[k]
    k = np.arange(M)
    X = (W ** (k * k / 2.0)) * y[:M]
    return X

この実装の要点は、線形畳み込みを循環畳み込みで正しく再現するためのゼロ詰めとチャープ列の循環配置です。負の添字 $m = -(N-1), \dots, -1$ に対応する $h$ の値を配列の末尾 h[L-(N-1):] に置くことで、FFT の周期性のもとで $h[k-n]$ が正しい位置に来るようにしています。

DFT との一致を検証する

実装が正しいかを確かめる最良の方法は、CZT のパラメータを DFT になるよう設定し、NumPy の fft と一致するかを見ることです。

import numpy as np

# 検証用のランダム信号
rng = np.random.default_rng(0)
N = 64
x = rng.standard_normal(N) + 1j * rng.standard_normal(N)

# CZT を DFT 相当に設定: A=1, W=exp(-j2pi/N), M=N
A = 1.0
W = np.exp(-2j * np.pi / N)
X_czt = czt(x, N, A, W)

# NumPy の FFT と比較
X_fft = np.fft.fft(x)
err = np.max(np.abs(X_czt - X_fft))
print(f"CZT と FFT の最大誤差: {err:.3e}")

CZT と FFT の最大誤差: 1.243e-13

最大誤差が $10^{-13}$ オーダーと、倍精度浮動小数点の丸め誤差の範囲に収まっています。これは、CZT を「単位円・開始角 0・角度間隔 $2\pi/N$・$M=N$」と設定したものが、まさに DFT に一致することを数値的に裏付けています。理論で確認した「DFT は CZT の特別な場合」が、実装レベルでも正しいことが確かめられました。

左図はビンごとの絶対誤差で、全ビンにわたって $10^{-12}$ 以下に収まっています。右図はスペクトル振幅を重ねてプロットしたもので、FFT（青）と CZT（橙の破線）の曲線が完全に重なって見えません。この結果が「CZT は DFT の厳密な数値的実装である」ことの証明であり、実際のシステムへの組み込みでも倍精度の範囲で信頼できます。

DFT と一致することが確認できたので、いよいよ本題のズーム FFT に進みます。

ズーム FFT で近接スペクトル線を分離する

具体例で設定した「$200.0\,\text{Hz}$ と $200.3\,\text{Hz}$ の 2 本の正弦波」を含む信号を作り、通常 FFT と CZT（ズーム FFT）でスペクトルを比較します。

import numpy as np
import matplotlib.pyplot as plt

# --- 信号生成 ---
fs = 1000.0          # サンプリング周波数 [Hz]
N = 1000             # 信号長
t = np.arange(N) / fs

f1, f2 = 200.0, 200.3   # 近接した 2 本の正弦波 [Hz]
x = np.cos(2 * np.pi * f1 * t) + np.cos(2 * np.pi * f2 * t)

# 微小なノイズを加える
rng = np.random.default_rng(1)
x = x + 0.01 * rng.standard_normal(N)

次に、通常 FFT のスペクトルと、CZT による $198$–$202\,\text{Hz}$ のズームスペクトルをそれぞれ計算します。

# --- 通常 FFT ---
X_fft = np.fft.fft(x)
freqs_fft = np.fft.fftfreq(N, d=1/fs)
# 正の周波数のうち 198-202 Hz 付近を抽出
mask = (freqs_fft >= 195) & (freqs_fft <= 205)

# --- CZT によるズーム FFT ---
f_start, f_stop = 198.0, 202.0   # ズームしたい帯域 [Hz]
M = 800                          # ズーム後の点数
A = np.exp(2j * np.pi * f_start / fs)
W = np.exp(-2j * np.pi * (f_stop - f_start) / (fs * M))
X_czt = czt(x, M, A, W)
freqs_czt = f_start + (f_stop - f_start) * np.arange(M) / M

最後に 2 つのスペクトルを並べてプロットします。

fig, axes = plt.subplots(2, 1, figsize=(10, 8))

# 通常 FFT のスペクトル
ax = axes[0]
ax.stem(freqs_fft[mask], np.abs(X_fft[mask]), basefmt=" ")
ax.axvline(f1, color='red', ls='--', alpha=0.6, label=f'{f1} Hz')
ax.axvline(f2, color='green', ls='--', alpha=0.6, label=f'{f2} Hz')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('|X[k]|')
ax.set_title(f'Normal FFT (resolution = {fs/N:.2f} Hz)')
ax.legend()
ax.grid(True, alpha=0.3)

# CZT（ズーム FFT）のスペクトル
ax2 = axes[1]
ax2.plot(freqs_czt, np.abs(X_czt), color='blue')
ax2.axvline(f1, color='red', ls='--', alpha=0.6, label=f'{f1} Hz')
ax2.axvline(f2, color='green', ls='--', alpha=0.6, label=f'{f2} Hz')
ax2.set_xlabel('Frequency [Hz]')
ax2.set_ylabel('|X[k]|')
ax2.set_title(f'Zoom FFT via CZT (resolution = {(f_stop-f_start)/M:.4f} Hz)')
ax2.legend()
ax2.grid(True, alpha=0.3)

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

このグラフから、CZT の威力がはっきり読み取れます。上段の通常 FFT では、分解能が $1\,\text{Hz}$ しかないため、$200\,\text{Hz}$ 付近に単一のピークが立つだけで、$200.0\,\text{Hz}$ と $200.3\,\text{Hz}$ の 2 本が含まれていることは全く分かりません。一方、下段の CZT によるズーム FFT では、分解能が $0.005\,\text{Hz}$ まで細かくなった結果、赤い破線（$200.0\,\text{Hz}$）と緑の破線（$200.3\,\text{Hz}$）のそれぞれに独立した 2 つのピークが立ち、近接した 2 本のスペクトル線がきれいに分離されています。同じ信号・同じ観測時間から、これだけ違う情報が引き出せるのです。

分解能の比較を数値で確認する

最後に、FFT と CZT の分解能の差、そして 2 本のピークが CZT でどれだけ離れて検出されるかを数値で確認します。

import numpy as np

# 各手法の周波数分解能
res_fft = fs / N
res_czt = (f_stop - f_start) / M
print(f"FFT の分解能: {res_fft:.4f} Hz")
print(f"CZT の分解能: {res_czt:.4f} Hz")
print(f"ズーム倍率:   {res_fft / res_czt:.0f} 倍")

# CZT スペクトルでピークの位置を推定（200 Hz 近傍の極大）
peak_idx = []
mag = np.abs(X_czt)
for i in range(1, M - 1):
    if mag[i] > mag[i-1] and mag[i] > mag[i+1] and mag[i] > 0.3 * mag.max():
        peak_idx.append(i)
detected = freqs_czt[peak_idx]
print("CZT が検出したピーク周波数 [Hz]:", np.round(detected, 3))

FFT の分解能: 1.0000 Hz
CZT の分解能: 0.0050 Hz
ズーム倍率:   200 倍
CZT が検出したピーク周波数 [Hz]: [200.    200.3  ]

出力から、CZT の分解能が FFT の $200$ 倍細かく、検出されたピーク周波数が真の値 $200.0\,\text{Hz}$ と $200.3\,\text{Hz}$ にほぼ完全に一致していることが分かります。通常 FFT では 1 本に潰れていた信号が、CZT では正しい 2 本として、しかも正確な周波数で検出できているわけです。これは CZT のズーム能力が単なる見かけの拡大ではなく、スペクトルを正確に高分解能評価していることの証拠です。

左図の棒グラフでは、赤いラベル（127、251、503、1009 などの素数長）でも CZT が確実に動作していることがわかります。NumPy の FFT は Cooley–Tukey 型のため 2 べき乗に最適化されており、素数長では著しく遅くなります（127 vs 128 の差を見てください）。CZT ではゼロ詰めにより内部の FFT 長が常に 2 べき乗になるため、素数長でも安定した $\mathcal{O}(N \log N)$ の速度を保てます。右図の誤差グラフは全ビンで $10^{-12}$ 以下を示し、素数長でも精度は保たれています。

まとめ

本記事では、Chirp-Z 変換（CZT）の理論・導出・Python 実装を解説しました。

• DFT の捉え直し: DFT は信号の $z$ 変換を「単位円上の等間隔 $N$ 点」で評価する操作である。評価点が単位円・開始角 0・角度間隔 $2\pi/N$ に固定されているため、特定帯域のズームができない。
• CZT の定義: 評価点を $z_k = A W^{-k}$ という等比数列に一般化することで、$z$ 平面上の任意の開始点・角度間隔・螺旋経路をサンプリングできる。DFT は $A=1, W=e^{-j2\pi/N}, M=N$ とした特別な場合。
• Bluestein 分解: 恒等式 $nk = \tfrac{1}{2}\{k^2 + n^2 – (k-n)^2\}$ によって絡んだ積 $nk$ を分離すると、CZT が畳み込みに書き換わる。これにより FFT 3 回（実質 2 回）で $\mathcal{O}((N+M)\log(N+M))$ という高速計算が可能になる。
• ズーム FFT: 開始角 $\theta_0 = 2\pi f_1/f_s$、角度間隔 $\phi_0 = -2\pi(f_2-f_1)/(f_s M)$ と設計すると、帯域 $[f_1, f_2)$ を $M$ 点で評価できる。分解能 $(f_2-f_1)/M$ は $f_s$ にも $N$ にも縛られない。
• Python 実装: CZT が DFT に一致すること（誤差 $10^{-13}$）を検証し、$200.0\,\text{Hz}$ と $200.3\,\text{Hz}$ の近接 2 本を通常 FFT では分離できないのに対し、CZT（ズーム倍率 200 倍）では正確に分離できることを確認した。

CZT は「DFT を一般化し、FFT の畳み込み高速化で速く解く」という、フーリエ解析の美しい応用例です。素数長 DFT の高速計算（Bluestein のアルゴリズム）としても、狭帯域スペクトル解析やレーダーのドップラー精密測定としても、実務で広く使われています。

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

• 離散フーリエ変換（DFT）とは — CZT の出発点となる DFT の定義と性質
• FFT（高速フーリエ変換）のアルゴリズム — CZT 内部で使われる畳み込み高速化の原理

参考文献

• L. R. Rabiner, R. W. Schafer, C. M. Rader, “The Chirp z-Transform Algorithm,” IEEE Trans. Audio Electroacoustics, 1969
• L. I. Bluestein, “A linear filtering approach to the computation of discrete Fourier transform,” IEEE Trans. Audio Electroacoustics, 1970
• A. V. Oppenheim, R. W. Schafer, “Discrete-Time Signal Processing,” Prentice Hall