チャープZ変換 — DFT を好きな場所にズームする

DFT は単位円上の等間隔評価にすぎない。Bluestein のトリックで畳み込みに変えると、任意の周波数帯域を高分解能で評価できる。


DFT(離散フーリエ変換)は、単位円上のN個の等間隔な点でZ変換を評価する操作だ。zk=ej2πk/Nz_k = e^{j2\pi k/N}(k=0,1,…,N−1k = 0, 1, \ldots, N-1)、これだけだ。

ではもし「特定の周波数帯域だけを高分解能で見たい」としたら? 大きなNのFFTを走らせて不要な部分を捨てる、というのが素朴な答えだが、もっと直接的な方法がある。**チャープZ変換(CZT)**だ。

Z変換を螺旋上で評価する

CZT は Z平面の螺旋上でZ変換を評価する。始点 z0z_0 と1ステップの進み方を表す複素数 WW を決めて、

zk=z0⋅W−k,k=0,1,…,M−1z_k = z_0 \cdot W^{-k}, \quad k = 0, 1, \ldots, M-1

という点列上でZ変換 X(z)=∑nxnz−nX(z) = \sum_{n} x_n z^{-n} を計算する:

Xk=∑n=0N−1xn z0−n WnkX_k = \sum_{n=0}^{N-1} x_n \, z_0^{-n} \, W^{nk}

∣W∣=1|W| = 1 なら単位円上の評価、∣W∣≠1|W| \neq 1 にすれば螺旋だ。DFTは z0=1z_0 = 1、W=e−j2π/NW = e^{-j2\pi/N}、M=NM = N の特殊ケースにすぎない。評価点の個数Mと入力長Nを独立に選べる点も DFT にはない自由度だ。

Bluestein のトリック

そのまま計算すると O(NM)O(NM) かかる。N点からM点を個別に計算するから当然だ。

だが Bluestein は1968年に気づいた——これは畳み込みに化ける、と。

鍵は恒等式:

nk=−(k−n)22+n22+k22nk = -\frac{(k-n)^2}{2} + \frac{n^2}{2} + \frac{k^2}{2}

これを使うと Wnk=W−(k−n)2/2⋅Wn2/2⋅Wk2/2W^{nk} = W^{-(k-n)^2/2} \cdot W^{n^2/2} \cdot W^{k^2/2} と分解できる。代入すると、

Xk=Wk2/2∑n=0N−1[xn z0−n Wn2/2]⏟gnW−(k−n)2/2⏟hk−nX_k = W^{k^2/2} \sum_{n=0}^{N-1} \underbrace{\left[ x_n \, z_0^{-n} \, W^{n^2/2} \right]}_{g_n} \underbrace{W^{-(k-n)^2/2}}_{h_{k-n}}

括弧内 gng_n は nn だけに依存し、hk−nh_{k-n} は差 (k−n)(k-n) だけに依存する。これは畳み込み (g⋆h)[k](g \star h)[k] の形だ:

Xk=Wk2/2⋅(g⋆h)[k]X_k = W^{k^2/2} \cdot (g \star h)[k]

畳み込みはFFTで O(Nlog⁡N)O(N \log N) で計算できる。CZTは「前処理 → FFTベースの畳み込み → 後処理」という構成になり、計算量は DFT と同じオーダーに収まる。

「チャープ」という名前は Wn2/2W^{n^2/2} の二次的な位相項が離散チャープ信号の形をしているところから来ている。

ズームFFTとして使う

RF の文脈で一番よく使われる形はズームFFTだ。

fs=100 MHzf_s = 100\,\text{MHz} のサンプリングで、1010〜12 MHz12\,\text{MHz} の帯域だけを 10241024 点で分析したい、という状況を考える。素朴にやるなら全帯域をカバーできる巨大なFFTが必要になる。

CZT なら設定は単純だ:

  • z0=ej2π⋅10 MHz/fsz_0 = e^{j2\pi \cdot 10\,\text{MHz}/f_s}(分析開始周波数)
  • W=e−j2π⋅2 MHz/(1024⋅fs)W = e^{-j2\pi \cdot 2\,\text{MHz}/(1024 \cdot f_s)}(周波数ステップ)
  • M=1024M = 1024

これで 1010〜12 MHz12\,\text{MHz} の帯域を 10241024 点、周波数分解能 ≈2 kHz\approx 2\,\text{kHz} で評価できる。スペクトラムアナライザーの「ズーム表示」の内部動作はこれに近い。

実装上の注意

CZT は評価点列の自由度が高い分、いくつか気をつける点もある。

出力の M 点は等間隔だが、DFT とは異なる周波数グリッドなので後段の処理と軸を合わせる必要がある。窓関数も入力側に素直にかければよい。NumPy には標準実装がないが、scipy.signal.czt が使える。あるいは Bluestein の展開を自分で書くと、アルゴリズムへの理解が深まるかもしれない。


一言でいうと「DFT は単位円上の評価、CZT はどこでも好きな場所の評価」ということだと思う。Bluestein のトリックで計算量は O(Nlog⁡N)O(N \log N) のままに抑えられる——シンプルな発想だが、きれいだなと感じる。

— ランキン

コメント

まだコメントはないよ。最初のひとことをどうぞ。