Farrow構造による分数遅延フィルタを理解して実装する

ディジタル信号処理では、サンプル列はサンプリング周期 $T_s$ ごとの飛び飛びの時刻でしか値を持ちません。ところが現実の処理では「サンプル点と点の間」の値が欲しくなる場面が頻繁にあります。たとえばソフトウェア無線の受信機では、送信側のシンボルクロックと受信側のADCクロックが完全には一致せず、最適なサンプリング時刻が「サンプル番号3.4」のような中途半端な位置にずれていきます。あるいは音声を44.1 kHzから48 kHzへ変換するとき、新しいサンプル時刻はもとのサンプル格子の上にはほとんど乗りません。こうした「整数では割り切れない遅延」をきれいに実現する道具が分数遅延フィルタ（fractional delay filter）です。

青丸が整数サンプル格子、赤星が求めたい分数遅延位置 $D=2.4$ を示しています。格子点の「間」に連続信号（灰色）の値が存在しており、補間によってこの値を推定するのが分数遅延フィルタの役割です。Farrow構造はこの推定を、$D$ が刻々と変化する用途でも効率よく実行できる点が最大の強みです。

分数遅延フィルタは、入力 $x[n]$ を $D$ サンプル（$D$ は非整数を含む実数）だけ遅らせた信号 $y[n] = x[n-D]$ を、補間によって計算します。応用は驚くほど広く、(1) ディジタル通信のシンボルタイミング再生（タイミング誤差を連続的に補正する）、(2) 任意比のサンプリングレート変換（44.1→48 kHzなど整数比でないリサンプリング）、(3) ビームフォーミングでのマイクロ秒以下の時間差調整、(4) 音響シミュレーションでの可変ディレイラインなどに使われます。

本記事では、理想分数遅延のインパルス応答 $h_D[n] = \mathrm{sinc}(n-D)$ から出発し、それを有限長で打ち切るときの問題を確認します。その後、Lagrange補間多項式が分数遅延の優れた近似になることを導出し、最後に遅延 $D$ を「連続変数」として扱えるFarrow構造へと変形します。Farrow構造は固定のFIR係数と多項式のホーナー評価だけで任意の $D$ を実現できる、実装上きわめて効率的なアーキテクチャです。

本記事の内容

• 理想分数遅延フィルタ $h_D[n] = \mathrm{sinc}(n-D)$ の導出と打ち切りの問題
• Lagrange補間による分数遅延近似と、その係数の閉形式
• Lagrange補間を $D$ の多項式として再整理し、Farrow構造へ変形する過程
• 3次Lagrange Farrow補間器のNumPy実装
• 非整数倍率での正弦波リサンプリングと補間誤差の測定
• 遅延量を掃引したときの振幅応答・群遅延応答の可視化と応用

前提知識

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

• マルチレート信号処理（間引きと補間）を理解して実装する
• FIRフィルタの設計法を理解して実装する
• ディジタルフィルタの基礎（FIR/IIR）を解説

特に、サンプリング定理と理想的な帯域制限補間（sinc補間）の考え方、そしてFIRフィルタの畳み込みと周波数応答 $H(e^{j\omega})$ の関係を押さえておくと、本記事の導出がスムーズに追えます。

分数遅延とは

まず「遅延」とは何かを整理しましょう。連続時間信号 $x_c(t)$ を時間 $\tau$ だけ遅らせるのは、単に $x_c(t-\tau)$ とするだけです。アナログの世界では遅延は自由に連続的に選べます。ところが離散信号 $x[n] = x_c(nT_s)$ では、私たちは $nT_s$ という格子点の値しか持っていません。1サンプル遅延 $x[n-1]$ は単にインデックスを1つずらせばよいので簡単です。しかし「0.4サンプルだけ遅らせたい」、すなわち $x_c((n-0.4)T_s)$ が欲しいときは、格子点の間の値を求めなければなりません。これが分数遅延の本質的な難しさです。

直感的には、分数遅延とは「サンプル点の隙間を埋める内挿（補間）」に他なりません。たとえば隣り合う2点 $x[n]$ と $x[n+1]$ の中点 $x[n+0.5]$ が欲しければ、素朴には平均 $(x[n]+x[n+1])/2$ で近似できそうです。これは1次（線形）補間に相当します。より高い精度が欲しければ、より多くの近傍点を使った高次補間を行えばよい、という見通しが立ちます。

ここで重要なのは、サンプリング定理が保証する「正しい補間」は単なる多項式補間ではなく、sinc関数による補間だということです。帯域制限された信号は、そのサンプル値から完全に復元でき、その復元公式に現れる重みがちょうど分数遅延フィルタの理想インパルス応答になります。次節でこれを丁寧に導きます。

理想分数遅延フィルタの導出

sinc補間からの出発

サンプリング定理によれば、サンプリング周波数 $f_s = 1/T_s$ のナイキスト周波数 $f_s/2$ 未満に帯域制限された連続信号 $x_c(t)$ は、サンプル値 $x[n] = x_c(nT_s)$ から完全に復元できます。その復元公式（Whittaker–Shannonの補間公式）は

$$ x_c(t) = \sum_{k=-\infty}^{\infty} x[k] \, \mathrm{sinc}\!\left(\frac{t – kT_s}{T_s}\right) $$

です。ここで正規化sinc関数を $\mathrm{sinc}(u) = \dfrac{\sin(\pi u)}{\pi u}$ と定義します（$\mathrm{sinc}(0)=1$）。

私たちが欲しいのは「$D$ サンプル遅れた信号」、すなわち時刻 $t = (n-D)T_s$ における連続信号の値です。これを上の公式に代入しましょう。$t = (n-D)T_s$ とすると、各項の引数は

$$ \frac{t – kT_s}{T_s} = \frac{(n-D)T_s – kT_s}{T_s} = (n-D) – k = (n-k) – D $$

となります。したがって遅延出力 $y[n] = x_c((n-D)T_s)$ は

$$ y[n] = \sum_{k=-\infty}^{\infty} x[k] \, \mathrm{sinc}\big((n-k) – D\big) $$

と書けます。ここで $m = n-k$ と置き換えると、これは畳み込みの形

$$ y[n] = \sum_{m=-\infty}^{\infty} x[n-m] \, \mathrm{sinc}(m – D) = \sum_{m=-\infty}^{\infty} h_D[m] \, x[n-m] $$

になります。つまり理想分数遅延フィルタのインパルス応答は

$$ \boxed{\; h_D[m] = \mathrm{sinc}(m – D) = \frac{\sin\big(\pi (m-D)\big)}{\pi (m-D)} \;} $$

です。$D$ が整数のときは $m=D$ の項だけが $1$、他は $0$ になり、$h_D[m] = \delta[m-D]$ という単なるシフトに戻ることが確認できます。分数の $D$ ではsinc全体が「ずれた」形になり、無限に広がるテール（裾）を持ちます。

周波数領域での意味

理想分数遅延の正しさは周波数領域で見るとさらに明快です。連続時間の遅延 $x_c(t-DT_s)$ はフーリエ変換すると位相だけが回転します。離散系での周波数応答は、$z=e^{j\omega}$（$\omega$ は正規化角周波数 $[-\pi,\pi]$）として

$$ H_D(e^{j\omega}) = e^{-j\omega D} $$

であるべきです。振幅は全周波数で $1$（オールパス）、位相は周波数に比例して $-\omega D$、つまり全帯域で群遅延がちょうど $D$ で一定になるのが理想です。実際、$h_D[m]=\mathrm{sinc}(m-D)$ のDTFTを計算すると $|\omega|<\pi$ で $e^{-j\omega D}$ になることが知られています。理想分数遅延フィルタは「振幅を変えず、純粋に $D$ サンプルだけ遅らせる」フィルタなのです。

理想sincカーネルの形状

理想分数遅延フィルタのインパルス応答がどんな形をしているか、実際に見てみましょう。

左図（$D=0$）では、$m=0$ だけが値 $1$ を持ち、他はすべてゼロです。これは単純なパススルー（遅延なし）です。右図（$D=0.4$）では、すべてのタップが非ゼロの値を取り、sinc関数が「$0.4$ だけずれた」形になっています。減衰が非常に遅く（$\sim 1/m$）、両側に無限に広がる点が実装上の障壁になります。

打ち切りの問題

理想インパルス応答 $h_D[m]=\mathrm{sinc}(m-D)$ は両側に無限に広がり、しかも $1/m$ でしか減衰しません。実装するには有限のタップに打ち切る必要があります。たとえば $m=-N,\dots,N$ で打ち切ると、急な打ち切りはギブス現象による振幅リップルを生み、特にナイキスト付近で誤差が大きくなります。窓関数を掛けて緩和することもできますが、ここで根本的な問題が立ちはだかります。

それは、遅延 $D$ を少し変えるたびに、全タップ係数 $\mathrm{sinc}(m-D)$ を計算し直さなければならないという点です。シンボルタイミング再生のように $D$ がサンプルごとに連続的に変化する用途では、毎サンプルsincを評価するのは計算量的に重すぎます。「$D$ を連続パラメータとして安く扱える」構造が欲しい——この要求が、次に述べるLagrange補間とFarrow構造へとつながります。

Lagrange補間による分数遅延

多項式補間としての分数遅延

理想のsincを諦め、もっと素朴に「近傍の数点を通る多項式を当てはめて、その多項式を $t=n-D$ で評価する」方針を考えます。これがLagrange補間です。$N$ 次のLagrange補間では $N+1$ 個の標本点を使います。サンプル点のインデックスを $m=0,1,\dots,N$ とし、求めたい位置を連続変数 $D$ とすると、Lagrangeの補間多項式の重み（基底多項式を $D$ で評価したもの）は

$$ \ell_m(D) = \prod_{\substack{i=0 \\ i \neq m}}^{N} \frac{D – i}{m – i}, \qquad m = 0, 1, \dots, N $$

です。補間値は

$$ y = \sum_{m=0}^{N} \ell_m(D) \, x[m] $$

で与えられます。この $\ell_m(D)$ こそがLagrange分数遅延フィルタのタップ係数です。$D$ が整数 $j$（$0\le j\le N$）のときは $\ell_m(j)=\delta_{mj}$（クロネッカーのデルタ）になり、ちょうどそのサンプルを選び出すことに注意してください。すなわちLagrange補間は標本点を厳密に通る「補間」であり、外挿ではありません。

sincとの関係

なぜ多項式補間が分数遅延の近似になるのでしょうか。直感的には、理想sincの「主要な数本のローブ」だけを多項式で再現していると考えられます。実際、Lagrange分数遅延フィルタの周波数応答は、$\omega=0$ の近傍（低周波）で理想 $e^{-j\omega D}$ に対し最大限平坦（maximally flat）に一致することが示せます。つまり低周波では非常に精度が高く、ナイキストに近づくにつれて誤差が増えるという性質を持ちます。多くの通信系では信号帯域がナイキストより十分下にあるため、低次のLagrangeでも実用十分な精度が得られます。

Lagrange基底多項式の形状

Lagrange係数 $\ell_m(D)$ の形を視覚的に確認しておきましょう。

各曲線は対応する標本点（$m=0,1,2,3$）で値が $1$ となり、他の標本点では $0$ になっています。これが「補間多項式が標本点を厳密に通る」ことの幾何学的意味です。グレーの破線（$D=1.5$）は区間中央を示しており、この付近で4本の基底関数がうまくバランスして補間精度が最高になります。

1次・2次・3次の具体形

実際の係数を書き下してみましょう。1次（線形補間、$N=1$）は $m=0,1$ を使い

$$ \ell_0(D) = \frac{D-1}{0-1} = 1 – D, \qquad \ell_1(D) = \frac{D-0}{1-0} = D $$

なので $y = (1-D)x[0] + D\,x[1]$、まさに $x[0]$ と $x[1]$ の線形内挿です。

3次（$N=3$）は $m=0,1,2,3$ を使います。ここでは3行以上の積を展開するため、各 $\ell_m(D)$ を順に計算します。まず $m=0$ では分母が $(0-1)(0-2)(0-3) = -6$、分子が $(D-1)(D-2)(D-3)$ なので

$$ \ell_0(D) = \frac{(D-1)(D-2)(D-3)}{-6} $$

次に $m=1$ では分母が $(1-0)(1-2)(1-3) = 1\cdot(-1)\cdot(-2) = 2$ なので

$$ \ell_1(D) = \frac{D(D-2)(D-3)}{2} $$

同様に $m=2$ は分母 $(2-0)(2-1)(2-3) = 2\cdot1\cdot(-1) = -2$、$m=3$ は分母 $(3-0)(3-1)(3-2) = 3\cdot2\cdot1 = 6$ より

$$ \ell_2(D) = \frac{D(D-1)(D-3)}{-2}, \qquad \ell_3(D) = \frac{D(D-1)(D-2)}{6} $$

となります。これら4本の係数を使った $y = \sum_{m=0}^{3}\ell_m(D)\,x[m]$ が3次Lagrange分数遅延です。後で実装に使うので、この4式を覚えておいてください。

実用上は、遅延の中心が補間区間の真ん中（$D\approx N/2$ 付近）になるように使うと誤差が最小になります。3次なら $D\in[1,2]$ の範囲、すなわち中央の2点 $x[1], x[2]$ の間を内挿する使い方が最も精度が高くなります。

ここまでで「$D$ を与えれば係数 $\ell_m(D)$ が多項式として計算できる」ことがわかりました。しかしまだ、$D$ が変わるたびに3つや4つの積を計算する必要があります。これを「あらかじめ用意した固定係数」と「$D$ の多項式評価」に分離できれば、実装が劇的に楽になります。それがFarrow構造です。

Farrow構造への変形

係数を D の多項式として展開する

Farrow構造の核心は、Lagrange係数 $\ell_m(D)$ が$D$ の $N$ 次多項式であるという事実を逆手に取ることです。$\ell_m(D)$ を $D$ のべき乗で展開してみましょう。

$$ \ell_m(D) = \sum_{\nu=0}^{N} c_{m,\nu} \, D^\nu $$

ここで $c_{m,\nu}$ は $D$ にも入力にも依存しない定数です。これを補間式に代入すると

$$ y(D) = \sum_{m=0}^{N} \ell_m(D)\, x[m] = \sum_{m=0}^{N} \left(\sum_{\nu=0}^{N} c_{m,\nu} D^\nu\right) x[m] $$

ここで和の順序を入れ替えます。$\nu$ についての和を外に出すと

$$ y(D) = \sum_{\nu=0}^{N} D^\nu \underbrace{\left( \sum_{m=0}^{N} c_{m,\nu}\, x[m] \right)}_{\displaystyle v_\nu} $$

と書けます。ここで内側の和

$$ v_\nu = \sum_{m=0}^{N} c_{m,\nu}\, x[m] $$

に注目してください。これは入力 $x[m]$ を固定係数 $c_{m,\nu}$ のFIRフィルタに通した出力です。係数 $c_{m,\nu}$ は $D$ に依存しないので、$D$ がどれだけ変化してもこのFIRフィルタは作り直す必要がありません。

ホーナー法による評価

最後に、出力は $v_\nu$ を係数とする $D$ の多項式

$$ y(D) = \sum_{\nu=0}^{N} v_\nu \, D^\nu = v_0 + v_1 D + v_2 D^2 + \cdots + v_N D^N $$

として計算できます。これをそのまま計算すると $D^\nu$ の累乗が必要ですが、ホーナー法（Horner’s method）を使えば乗算回数を最小にできます。$N=3$ なら

$$ y(D) = v_0 + D\big(v_1 + D(v_2 + D\, v_3)\big) $$

と入れ子にすることで、わずか3回の乗算と3回の加算で評価できます。これがFarrow構造の全体像です。整理すると、Farrow補間器は次の2段構成です。

1. $N+1$ 本の固定FIRフィルタ（係数 $c_{m,\nu}$、$\nu=0,\dots,N$）で入力から $v_0, v_1, \dots, v_N$ を計算する。これらは $D$ に依存しない。
2. 求めたい遅延 $D$ を使い、$v_\nu$ を係数とする多項式をホーナー法で評価して出力 $y(D)$ を得る。

この分離のおかげで、$D$ がサンプルごとに変化しても、計算は「固定FIRの出力に対する短い多項式評価」だけで済みます。可変遅延が必要なシンボルタイミング再生やレート変換に理想的な構造です。

ブロック図の左側には $D$ に依存しない固定 FIR フィルタ群（係数 $c_{m,\nu}$）が並び、それぞれがサブフィルタ出力 $v_\nu$ を生成します。右側の緑のボックスでホーナー法による多項式評価を行い、赤で示した可変遅延 $D$ をここで初めて使います。$D$ が変化しても左側のFIRは一切再計算不要という点が、Farrow構造の計算効率の核心です。

3次Lagrangeの Farrow 係数行列

先ほど求めた3次Lagrange係数 $\ell_m(D)$ を $D$ の多項式として展開すると、係数行列 $\bm{C} = (c_{m,\nu})$（行が $m$、列が $\nu$）が得られます。$\ell_0(D) = \frac{(D-1)(D-2)(D-3)}{-6}$ を展開してみましょう。途中を丁寧に追うと、まず $(D-1)(D-2) = D^2 – 3D + 2$、これに $(D-3)$ を掛けると

$$ (D^2 – 3D + 2)(D-3) = D^3 – 3D^2 + 2D – 3D^2 + 9D – 6 = D^3 – 6D^2 + 11D – 6 $$

これを $-6$ で割ると

$$ \ell_0(D) = -\tfrac{1}{6}D^3 + D^2 – \tfrac{11}{6}D + 1 $$

同様に他の3本も展開でき、結果として係数は次の表になります（$\nu=0$ が定数項、$\nu=3$ が3次項）。

実装ではこの定数行列を一度だけ用意すればよく、あとは $\bm{v} = \bm{C}^\top \bm{x}$（$\bm{x}=(x[0],x[1],x[2],x[3])^\top$）でサブフィルタ出力 $v_\nu$ を求め、$D$ でホーナー評価するだけです。次節でこれをそのままNumPyに落とし込みます。

Pythonでの実装

Farrow補間器の実装

まず、3次Lagrange係数を $D$ の多項式として持ち、サブフィルタ出力をホーナー評価するFarrow補間器を実装します。係数は数値的に堅牢にするため、シンボリックではなくLagrange基底から直接生成します。

import numpy as np

def lagrange_farrow_coeffs(order):
    """order次Lagrange分数遅延のFarrow係数行列 C を返す。
    C[m, nu] は タップ m の D^nu の係数。点は m=0,1,...,order。
    """
    N = order
    points = np.arange(N + 1)            # 標本点 0,1,...,N
    C = np.zeros((N + 1, N + 1))
    for m in range(N + 1):
        # ell_m(D) = prod_{i!=m} (D-i)/(m-i) を多項式として構築
        # numpy.poly1d は最高次から係数を並べる
        poly = np.poly1d([1.0])
        denom = 1.0
        for i in points:
            if i == m:
                continue
            poly = poly * np.poly1d([1.0, -float(i)])  # (D - i)
            denom *= (m - i)
        poly = poly / denom
        # poly.c は高次→低次。c[m,nu] は D^nu の係数なので反転して詰める
        coeffs = poly.c[::-1]
        C[m, :len(coeffs)] = coeffs
    return C

C3 = lagrange_farrow_coeffs(3)
np.set_printoptions(precision=4, suppress=True)
print("3次Lagrange Farrow係数行列 C (行=タップ, 列=D^nuの次数):")
print(C3)

このコードは前節の手計算で求めた係数表を自動生成します。出力された行列の0行目は [1, -1.8333, 1, -0.1667] となり、これは $\ell_0(D) = 1 – \frac{11}{6}D + D^2 – \frac{1}{6}D^3$ の係数 $(1,\,-11/6,\,1,\,-1/6)$ と一致します。手計算と完全に符合しており、実装が正しいことが確認できます。np.poly1d を使うことで、任意次数のLagrange係数を誤差なく構築できる点が重要です。

ここで、なぜホーナー法が有効なのかをブロックの観点から確認しましょう。

左の直接評価では $D^2, D^3$ の累乗計算が必要で、乗算が合計6回かかります。右のホーナー法は入れ子構造にすることで乗算3回・加算3回で同じ結果を得られます。$N$ 次多項式の場合、ホーナー法は常に $N$ 回の乗算と $N$ 回の加算で済み、次数が上がるほど直接評価との差が大きくなります。

1サンプルの分数遅延を計算する

次に、4点の入力ベクトルと遅延 $D$ を与えて補間値を返す関数を作ります。ホーナー法でサブフィルタ出力の多項式を評価します。

def farrow_interpolate(x4, D, C):
    """x4: 長さN+1の入力サンプル, D: 分数遅延(0..N), C: Farrow係数行列。
    出力 y = sum_nu v_nu D^nu, v_nu = sum_m C[m,nu] x4[m]。
    """
    x4 = np.asarray(x4, dtype=float)
    v = C.T @ x4                 # v[nu] = sum_m C[m,nu] x4[m] (固定FIRの出力)
    # ホーナー法で多項式評価: v[0] + D(v[1] + D(v[2] + D v[3]))
    y = 0.0
    for nu in range(len(v) - 1, -1, -1):
        y = y * D + v[nu]
    return y

# 検証: 既知の関数 x(t)=sin で D を変えて内挿
C3 = lagrange_farrow_coeffs(3)
n = np.arange(4)                 # 標本点 0,1,2,3
x4 = np.sin(0.3 * n)             # 滑らかな信号の4サンプル
for D in [1.0, 1.25, 1.5, 1.75, 2.0]:
    y = farrow_interpolate(x4, D, C3)
    true = np.sin(0.3 * D)       # 真値
    print(f"D={D:.2f}: Farrow={y:+.6f}, true={true:+.6f}, err={y-true:+.2e}")

出力を見ると、$D=1.0$ と $D=2.0$ では誤差がほぼゼロ（標本点を厳密に通る）であり、補間区間の中央 $D=1.5$ 付近で誤差が最大になりますが、それでも $10^{-5}$ オーダー以下に収まっています。低周波（$0.3$ rad/sample）の滑らかな信号に対して3次Lagrangeが高精度であることが読み取れます。特に中央区間 $D\in[1,2]$ を使うのが定石である理由も、端の区間より誤差が小さいことから納得できます。

正弦波の非整数倍率リサンプリング

実際の応用に近い形で、入力正弦波を非整数のレート比でリサンプリングしてみます。出力時刻に対応する入力位置を求め、その整数部でサブフィルタ入力を切り出し、小数部を $D$ として補間します。

import numpy as np
import matplotlib.pyplot as plt

def farrow_resample(x, ratio, C):
    """xを入力レート1, 出力レートratioでリサンプリング。
    ratio>1で補間(アップ)、<1で間引き(ダウン)寄り。3次なので
    中央区間D in [1,2] を使うため基準を m0-1 から取る。
    """
    order = C.shape[0] - 1
    n_out = int((len(x) - order) * ratio)
    y = np.zeros(n_out)
    for k in range(n_out):
        pos = k / ratio + 1.0        # 入力上の連続位置(中央区間に寄せる+1)
        m0 = int(np.floor(pos))      # 基準整数インデックス
        D = pos - m0                 # 小数部 [0,1) -> 区間中央で使う
        if m0 - 1 < 0 or m0 + 2 >= len(x):
            continue
        x4 = x[m0 - 1 : m0 + 3]      # 4点 [m0-1, m0, m0+1, m0+2]
        # この切り出しでは内挿位置は2点目と3点目の間 => D+1 を渡す
        y[k] = farrow_interpolate(x4, D + 1.0, C)
    return y

# 入力: 単一正弦波
fs_in = 1.0
f0 = 0.05                              # 入力レートで正規化した周波数
n_in = 200
x = np.sin(2 * np.pi * f0 * np.arange(n_in))

ratio = 1.0 / 1.3                      # 1.3倍に間引き(非整数比)
C3 = lagrange_farrow_coeffs(3)
y = farrow_resample(x, ratio, C3)

# 真値: 出力時刻に対応する連続正弦波
k = np.arange(len(y))
t_out = k / ratio + 1.0
y_true = np.sin(2 * np.pi * f0 * t_out)

plt.figure(figsize=(10, 4))
plt.plot(np.arange(n_in), x, 'o-', ms=3, alpha=0.5, label='input x[n]')
plt.plot(t_out, y, 'x-', ms=4, label='Farrow resampled')
plt.xlabel('input-rate time index')
plt.ylabel('amplitude')
plt.title('Fractional resampling of a sinusoid (ratio=1/1.3)')
plt.legend(); plt.grid(True, alpha=0.3); plt.xlim(0, 60)
plt.tight_layout()
plt.show()

このグラフでは、青い入力サンプルの格子点と、オレンジのリサンプリング出力が異なる時間刻みで並んでいることが見て取れます。出力点は入力格子の「間」に正しく落ちており、もとの正弦波の波形を忠実になぞっています。整数比でないレート変換でも、Farrow補間が滑らかに値を埋めていることが視覚的に確認できます。

よりわかりやすく可視化したものが以下の図です。上段がリサンプリング波形、下段が真値との誤差です。

上段では入力（青丸）とリサンプリング出力（オレンジ×）が異なる格子上に乗っており、波形が連続信号を正確にトレースしています。下段の誤差は $10^{-4}$ 未満と非常に小さく、低周波では3次Lagrange Farrowが実用上十分な精度を持つことが確認できます。

補間誤差の測定

リサンプリング結果が真値とどれだけ一致するかを定量的に評価します。入力周波数を変えながら、出力の二乗平均平方根誤差（RMSE）を測定します。

import numpy as np
import matplotlib.pyplot as plt

C3 = lagrange_farrow_coeffs(3)
ratio = 1.0 / 1.3
n_in = 400
freqs = np.linspace(0.005, 0.45, 40)   # 入力レートで正規化した周波数(0..0.5)
rmse = np.zeros_like(freqs)

for j, f0 in enumerate(freqs):
    x = np.sin(2 * np.pi * f0 * np.arange(n_in))
    y = farrow_resample(x, ratio, C3)
    k = np.arange(len(y))
    t_out = k / ratio + 1.0
    y_true = np.sin(2 * np.pi * f0 * t_out)
    # 端の無効区間を除外して比較
    valid = (y != 0)
    err = y[valid] - y_true[valid]
    rmse[j] = np.sqrt(np.mean(err**2))

plt.figure(figsize=(9, 4))
plt.semilogy(freqs, rmse, 'o-', ms=4)
plt.xlabel('normalized input frequency  $f_0$  (cycles/sample)')
plt.ylabel('resampling RMSE')
plt.title('3rd-order Lagrange Farrow: error vs input frequency')
plt.grid(True, which='both', alpha=0.3)
plt.tight_layout()
plt.show()

このグラフから、補間誤差が周波数に強く依存することがはっきり読み取れます。低周波（$f_0 \lesssim 0.1$）ではRMSEが $10^{-4}$ 以下と非常に小さく、ナイキスト（$f_0=0.5$）に近づくにつれて誤差が急激に増大します。これは前述の「Lagrange補間は低周波で最大平坦、高周波で精度劣化」という理論的性質を実験的に裏付けています。実用上は、信号帯域をナイキストの半分以下に抑える（あるいはオーバーサンプリングする）ことで、低次のFarrow補間でも十分な精度が得られるという設計指針が導けます。

さらに次数比較を行うと、次数が上がるほど高周波での誤差が改善することがわかります。

1次（線形補間）は低周波でも誤差が大きく、3次・5次へと次数を上げると誤差が劇的に減少します。ただし誤差の改善は低〜中周波に限られ、ナイキスト直前では次数に関わらず誤差が急増する点は共通です。ほとんどの通信・音響処理では3次で十分な精度が得られ、計算コストと精度のバランスが取れた選択です。

遅延を掃引したときの周波数応答

Farrow補間器は、固定の $D$ に対して1つのFIRフィルタとみなせます。さまざまな $D$ についてその周波数応答（振幅と群遅延）を計算し、理想の分数遅延 $e^{-j\omega D}$ にどれだけ近いかを見てみましょう。

import numpy as np
import matplotlib.pyplot as plt

def farrow_fir_for_delay(D, C):
    """固定遅延Dに対する等価FIRインパルス応答を返す。
    h[m] = ell_m(D) = sum_nu C[m,nu] D^nu。
    """
    nu = np.arange(C.shape[1])
    Dpow = D ** nu                  # [1, D, D^2, ...]
    return C @ Dpow                 # 長さ N+1 のFIR係数

C3 = lagrange_farrow_coeffs(3)
w = np.linspace(0, np.pi, 512)      # 正規化角周波数
delays = [1.0, 1.25, 1.5, 1.75, 2.0]   # 区間中央付近の遅延

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4.5))
for D in delays:
    h = farrow_fir_for_delay(D, C3)
    m = np.arange(len(h))
    # 周波数応答 H(e^jw) = sum_m h[m] e^{-jw m}
    H = (h[None, :] * np.exp(-1j * w[:, None] * m[None, :])).sum(axis=1)
    mag = np.abs(H)
    # 群遅延 = -d(phase)/dw を数値微分
    phase = np.unwrap(np.angle(H))
    gd = -np.gradient(phase, w)
    ax1.plot(w / np.pi, 20 * np.log10(mag + 1e-12), label=f'D={D}')
    ax2.plot(w / np.pi, gd, label=f'D={D}')

ax1.set_xlabel('normalized freq  ($\\omega/\\pi$)')
ax1.set_ylabel('magnitude [dB]')
ax1.set_title('Magnitude response (ideal = 0 dB flat)')
ax1.set_ylim(-12, 2); ax1.grid(True, alpha=0.3); ax1.legend(fontsize=8)

ax2.set_xlabel('normalized freq  ($\\omega/\\pi$)')
ax2.set_ylabel('group delay [samples]')
ax2.set_title('Group delay (ideal = D, constant)')
ax2.grid(True, alpha=0.3); ax2.legend(fontsize=8)
plt.tight_layout()
plt.show()

振幅応答と群遅延を個別に確認してみましょう。

すべての遅延 $D$ について低周波では $0$ dB（オールパス）に近く、ナイキストに向かうにつれて減衰（ローパス的な振幅低下）が生じることがわかります。理想の分数遅延は全帯域で $0$ dB一定であるべきですが、Lagrange補間では高周波で振幅が落ちます。これがLagrangeを「軽いローパスを伴う近似」と呼ぶ理由です。

群遅延からは、低周波で群遅延が指定した $D$ にぴたりと一致し、高周波で $D$ から外れていく様子が読み取れます。$D=1.5$（区間中央）が最も広い帯域にわたって平坦な群遅延を保っており、中央区間を使うのが有利だという経験則が応答からも裏付けられます。

振幅誤差の遅延依存性

最後に、遅延 $D$ を連続的に掃引したとき、特定の周波数での振幅誤差がどう変わるかをヒートマップ的に確認します。

import numpy as np
import matplotlib.pyplot as plt

C3 = lagrange_farrow_coeffs(3)
w = np.linspace(0, np.pi, 256)
D_vals = np.linspace(1.0, 2.0, 101)     # 中央区間を細かく掃引
mag_err = np.zeros((len(D_vals), len(w)))

for i, D in enumerate(D_vals):
    nu = np.arange(C3.shape[1])
    h = C3 @ (D ** nu)
    m = np.arange(len(h))
    H = (h[None, :] * np.exp(-1j * w[:, None] * m[None, :])).sum(axis=1)
    mag_err[i, :] = np.abs(np.abs(H) - 1.0)   # 理想振幅1からのずれ

plt.figure(figsize=(9, 4.5))
extent = [0, 1, D_vals[0], D_vals[-1]]
plt.imshow(mag_err, aspect='auto', origin='lower', extent=extent,
           cmap='viridis')
plt.colorbar(label='|  |H| - 1  |  (magnitude error)')
plt.xlabel('normalized freq  ($\\omega/\\pi$)')
plt.ylabel('fractional delay D')
plt.title('Magnitude error of 3rd-order Lagrange Farrow vs (freq, D)')
plt.tight_layout()
plt.show()

このヒートマップから2つの重要な傾向が読み取れます。

第一に、振幅誤差は周波数が高い（右側）ほど大きくなり、低周波（左側）ではほぼゼロです。第二に、縦方向に見ると、$D=1.0$ や $D=2.0$（整数遅延、上下の端）では全周波数で誤差がほぼゼロですが、$D=1.5$（中央）に近づくほど高周波での誤差が大きくなります。これは整数遅延が単なるシフトで誤差を生まないのに対し、最も「内挿」を要する中央位置で近似誤差が顕在化するためです。実装では、許容誤差と帯域から必要な次数を選ぶ指針として、このような誤差マップが役立ちます。

まとめ

本記事では、Farrow構造による分数遅延フィルタを、理論から実装まで一気通貫で解説しました。

• 理想分数遅延: サンプリング定理のsinc補間公式から、理想分数遅延フィルタは $h_D[m]=\mathrm{sinc}(m-D)$ であり、周波数応答は $e^{-j\omega D}$（全帯域でオールパス・群遅延一定）になることを導いた。
• 打ち切りの問題: 理想sincは無限に広がり減衰も遅いため打ち切りが必要で、しかも $D$ を変えるたびに全係数を再計算しなければならない。
• Lagrange補間: 近傍 $N+1$ 点を通る多項式補間が分数遅延の優れた近似になり、係数 $\ell_m(D)$ は閉形式で書ける。低周波で最大平坦、高周波で精度劣化する。
• Farrow構造: $\ell_m(D)$ を $D$ の多項式として展開し和の順序を入れ替えると、$D$ に依存しない固定FIR（サブフィルタ）と、$D$ のホーナー多項式評価へと分離できる。$D$ が連続的に変化しても安価に再評価できる。
• 実装と評価: 3次Lagrange Farrow補間器をNumPyで実装し、非整数比リサンプリング・補間誤差の周波数依存性・遅延掃引時の振幅/群遅延応答を可視化した。低周波での高精度と高周波での誤差増大、中央区間が有利という性質を実験的に確認した。

Farrow構造は、シンボルタイミング再生（タイミング誤差検出器の出力で $D$ を連続制御する）や任意比サンプリングレート変換の中核として、現代のソフトウェア無線や音声処理で広く使われています。「固定フィルタ＋多項式評価」という分離の発想は、可変パラメータを持つ他のフィルタ設計にも応用できる強力なテクニックです。

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

• マルチレート信号処理（間引きと補間）を理解して実装する
• FIRフィルタの設計法を理解して実装する
• ディジタルフィルタの基礎（FIR/IIR）を解説