パルス整形と整合フィルタ

前のレッスンでは、ビットをコンステレーション上の点に対応付け、各シンボルを1つの数として扱いました。無線機は数をそのまま送信できません。送信するのは波形です。各シンボルを運ぶパルスの形は、2つのことを決めます。信号が占有するスペクトラムの広さと、受信側で隣り合うシンボルが互いに干渉するかどうかです。このレッスンは、numpyを実行でき、この2つの効果を実際に測定して確かめたい読者を対象としています。ルートレイズドコサインフィルタを式から作り、スペクトラムを比較し、受信側に整合フィルタを加えます。最後に、モデムのエンジニアが最初の確認に使うアイパターンを読みます。

矩形パルスの問題

シンボルを送る最も素直な方法は、その値を1シンボル周期Tの間保持することです。これが矩形パルスです。そのスペクトラムはsinc関数になります。メインローブの幅は2/Tで、サイドローブの減衰は緩やかです。第1サイドローブはピークから約13 dBしか下がらず、その後の包絡線も1オクターブあたり6 dBずつしか下がりません。チャネル割当ては固定された幅のスペクトラムなので、これらのサイドローブのエネルギーは隣接チャネルに漏れ込みます。規制当局や標準規格がスペクトラムマスクを定めているのは、まさにこれを防ぐためです。

反対の極端は、時間領域のsincパルスです。そのスペクトラムは幅1/Tの完全な矩形です。これは、シンボルレート1/Tの信号が情報を失わずに取りうる最も狭い帯域です。これより狭い帯域で、レート1/Tの独立なシンボルを干渉なしに運べる信号はありません。問題は、sincパルスが無限に続き、その裾が1/tでしか減衰しないことです。受信側にタイミング誤差があると多数の隣接シンボルからの寄与を拾い、その和が大きくなることがあります。

求めるものはその中間にあります。スペクトラムがコンパクトで、裾が速く減衰し、なおかつ他のすべてのシンボル時刻でゼロ交差をもつパルスです。この条件は符号間干渉ゼロのためのナイキスト基準で、Harry Nyquistにちなんで名付けられています。Nyquistによる1928年の電信伝送に関する論文は末尾に挙げています。サンプリングのナイキストレートとは名前が同じですが、別の考え方です。全体のパルスp(t)が0以外のすべてのkについてt = kTでゼロであれば、正しい時刻でサンプリングしたときに見えるのは現在のシンボルだけです。

レイズドコサインとルートレイズドコサイン

レイズドコサインパルスは、調整可能なトレードオフを持たせながらこの基準を満たします。そのスペクトラムは(1 − β)/(2T)まで平坦で、その後コサインの半周期の形で減衰し、(1 + β)/(2T)でゼロになります。ロールオフ率βは0から1の範囲をとります。β = 0ではパルスはsincになります。β = 1ではスペクトラムの幅は最小値の2倍になり、時間領域の裾は速く消えます。複素ベースバンド信号の占有帯域幅の合計は(1 + β)/Tです。したがって、β = 0.35の1 Mbaudの信号には1.35 MHzが必要です。

実際のリンクでは、レイズドコサインの応答を2つに分けます。送信機はルートレイズドコサイン（RRC）フィルタを適用します。その周波数応答はレイズドコサインの平方根です。受信機も同じRRCフィルタを適用します。両者の積はレイズドコサインになるので、縦続接続全体の符号間干渉（ISI）はゼロです。分割する理由は、後の整合フィルタの節で説明します。

このブロックは、RRCのインパルス応答を閉形式の式から作ります。1シンボルあたり8サンプル、スパンは12シンボルです。1つのRRCと、2つを縦続接続したものについて、シンボル時刻でのISIを確認します。次に16,384個のランダムなBPSKシンボルを生成し、矩形パルスとRRCパルスでそれぞれ整形して、両方のスペクトラムを推定します。

import numpy as np

def rrc_taps(beta: float, sps: int, span: int) -> np.ndarray:
    """Root-raised-cosine impulse response, span symbols long, unit energy."""
    t = np.arange(-span * sps // 2, span * sps // 2 + 1) / sps   # time in symbols
    h = np.zeros_like(t)
    for i, ti in enumerate(t):
        if np.isclose(ti, 0.0):
            h[i] = 1.0 - beta + 4 * beta / np.pi
        elif beta > 0 and np.isclose(abs(ti), 1 / (4 * beta)):
            h[i] = (beta / np.sqrt(2)) * ((1 + 2 / np.pi) * np.sin(np.pi / (4 * beta))
                                          + (1 - 2 / np.pi) * np.cos(np.pi / (4 * beta)))
        else:
            num = np.sin(np.pi * ti * (1 - beta)) + 4 * beta * ti * np.cos(np.pi * ti * (1 + beta))
            h[i] = num / (np.pi * ti * (1 - (4 * beta * ti) ** 2))
    return h / np.sqrt(np.sum(h**2))

sps, span, beta = 8, 12, 0.35          # 8 samples per symbol, 12-symbol filter, roll-off 0.35
h_rrc = rrc_taps(beta, sps, span)
h_rc = np.convolve(h_rrc, h_rrc)       # RRC at the transmitter, RRC at the receiver
centre = len(h_rc) // 2

def isi_at_symbol_instants(p: np.ndarray, c: int) -> float:
    """Largest |p(kT)| for k != 0, relative to p(0)."""
    k = np.arange(-(c // sps), c // sps + 1) * sps + c
    k = k[(k >= 0) & (k < len(p)) & (k != c)]
    return float(np.max(np.abs(p[k])) / abs(p[c]))

print(f"RRC filter: {len(h_rrc)} taps, peak tap {h_rrc.max():.4f}")
print(f"RRC alone, worst ISI sample: {isi_at_symbol_instants(h_rrc, len(h_rrc) // 2):.4f}")
print(f"RRC * RRC (raised cosine), worst ISI sample: {isi_at_symbol_instants(h_rc, centre):.5f}")

# Power spectral density of a random BPSK stream with a rectangular pulse and with the RRC pulse.
rng = np.random.default_rng(6)
n_sym = 16384
symbols = rng.choice([-1.0, 1.0], n_sym)
up = np.zeros(n_sym * sps)
up[::sps] = symbols
x_rect = np.convolve(up, np.ones(sps) / np.sqrt(sps))[: n_sym * sps]
x_rrc = np.convolve(up, h_rrc)[: n_sym * sps]

def psd_db(x: np.ndarray, nfft: int = 256) -> tuple[np.ndarray, np.ndarray]:
    segs = x[: len(x) // nfft * nfft].reshape(-1, nfft) * np.hanning(nfft)
    p = np.mean(np.abs(np.fft.fftshift(np.fft.fft(segs, axis=1), axes=1)) ** 2, axis=0)
    f = np.fft.fftshift(np.fft.fftfreq(nfft, d=1 / sps))   # frequency in units of the symbol rate
    return f, 10 * np.log10(p / p.max())

f, p_rect = psd_db(x_rect)
_, p_rrc = psd_db(x_rrc)

def occupied_bw(x: np.ndarray, frac: float = 0.99) -> float:
    X = np.abs(np.fft.fftshift(np.fft.fft(x))) ** 2
    fr = np.fft.fftshift(np.fft.fftfreq(len(x), d=1 / sps))
    keep = fr >= 0                                        # real signal: one side is enough
    cum = np.cumsum(X[keep]) / np.sum(X[keep])
    return 2 * float(fr[keep][np.searchsorted(cum, frac)])

print(f"99 % power bandwidth, rectangular: {occupied_bw(x_rect):.2f} x symbol rate")
print(f"99 % power bandwidth, RRC beta={beta}: {occupied_bw(x_rrc):.2f} x symbol rate (edge-to-edge width 1+beta = {1 + beta:.2f})")
print("f/Rs, rect dB, RRC dB")
for i in range(128, 194, 2):            # 0 to 2 times the symbol rate
    print(f"{f[i]:.3f}, {p_rect[i]:.1f}, {p_rrc[i]:.1f}")

RRC filter: 97 taps, peak tap 0.3874
RRC alone, worst ISI sample: 0.0773
RRC * RRC (raised cosine), worst ISI sample: 0.00120
99 % power bandwidth, rectangular: 6.88 x symbol rate
99 % power bandwidth, RRC beta=0.35: 1.16 x symbol rate (edge-to-edge width 1+beta = 1.35)
f/Rs, rect dB, RRC dB
0.000, -0.0, -0.2
0.062, -0.2, -0.5
0.125, -0.1, -0.3
0.188, -0.4, 0.0
0.250, -1.0, -0.4
0.312, -1.3, -0.3
0.375, -2.4, -0.3
0.438, -2.8, -1.2
0.500, -4.0, -3.3
0.562, -4.9, -6.2
0.625, -6.7, -12.9
0.688, -8.1, -28.5
0.750, -10.4, -42.1
0.812, -12.9, -45.3
0.875, -16.7, -48.8
0.938, -23.1, -52.7
1.000, -34.2, -54.8
1.062, -24.3, -54.4
1.125, -18.8, -54.3
1.188, -16.1, -56.3
1.250, -14.7, -59.4
1.312, -13.4, -58.7
1.375, -13.3, -57.6
1.438, -12.7, -58.6
1.500, -13.1, -61.8
1.562, -13.3, -61.4
1.625, -14.5, -59.8
1.688, -15.4, -60.5
1.750, -17.2, -63.1
1.812, -19.3, -63.5
1.875, -22.7, -61.6
1.938, -28.8, -61.9
2.000, -39.6, -63.7

単一のRRCパルスはナイキストパルスではありません。隣接シンボル時刻での最悪のサンプルは、ピークの7.7%です。2つ目のRRCを通すと最悪の残留は0.12%になります。これはフィルタを12シンボルで打ち切ったことによるものです。スパンを長くすればタップ数は増えますが、この値をさらに下げられます。

帯域幅については、スペクトラムが次のことを示しています。β = 0.35では、電力の99%がシンボルレートの1.16倍の範囲に収まり、理想的な端から端までの幅1.35の内側にあります。矩形パルスは、同じ99%を得るのにシンボルレートの少なくとも6.88倍を必要とします。この値は、シミュレーションがシンボルレートの8倍の範囲しか扱っていないことで頭打ちになっています。中心からシンボルレートの0.75倍の位置では、RRCのスペクトラムは42.1 dB下がっていますが、矩形の方は10.4 dBしか下がっていません。矩形のスペクトラムのシンボルレートの約1.4倍付近のサイドローブは、sincから予想されるとおり−13.3 dBにあります。

RRCの曲線は、どこまでも下がり続けるのではなく、−55 dBから−65 dBの間で平坦になります。インパルス応答を12シンボルで打ち切ると阻止域に下限ができ、スペクトラム推定の窓関数もわずかな漏れを加えます。スペクトラムマスクがさらに大きな抑圧を求める場合は、フィルタを長くしてください。

符号間干渉

符号間干渉（ISI）とは、現在のシンボルをサンプリングする時刻に存在する、他のシンボルからのエネルギーです。よくある原因は3つあります。第一に、単一のRRCのように、パルス自体がナイキスト基準を満たしていない場合です。第二に、受信機が誤った時刻でサンプリングし、隣接パルスの裾が干渉になる場合です。第三に、チャネルがパルスを広げる場合で、これはレッスン7でマルチパスとして扱います。

ISIは、1つの重要な点で雑音と異なります。雑音の相対的な電力は、送信電力を上げれば下がります。ISIは信号に比例して大きくなるので、電力を上げても取り除けません。ISIで制限されたリンクにはエラーフロアが現れます。あるSNRを超えると、ビット誤り率がそれ以上改善しなくなります。

整合フィルタ

受信機には、送信されたパルスに白色雑音が加わった信号が届きます。各シンボルを判定する前に、受信機は信号をフィルタに通し、1シンボルに1回サンプリングします。サンプリング時刻でのSNRを最大にするフィルタが整合フィルタです。そのインパルス応答は、送信パルスを時間反転した複素共役です。エネルギーEのパルスとスペクトル密度N0/2の白色雑音に対して、サンプルで達成できる最良のSNRは2E/N0です。これに到達できるのは整合フィルタだけです。

RRCパルスは対称な実数値なので、その整合フィルタは同じRRCです。レイズドコサインを分割するのはこのためです。両端に1つずつRRCを置くと最良のSNRが得られます。さらに2つの積がレイズドコサインになるので、ISIもゼロになります。1つのフィルタ設計で両方の性質が得られます。

2つ目のブロックでは、Eb/N0が6 dBの条件で、400,000ビットのBPSKを使ってこれを確認します。受信機Aは整合フィルタを適用し、1シンボルに1回サンプリングします。受信機Bはフィルタを使わず、雑音を含む波形をパルスのピークでサンプリングします。同じブロックで次に雑音を取り除き、3つのロールオフ率について、タイミングオフセットに対するアイの開口を測定します。

import numpy as np
from scipy.special import erfc

def rrc_taps(beta: float, sps: int, span: int) -> np.ndarray:
    """Root-raised-cosine impulse response, span symbols long, unit energy."""
    t = np.arange(-span * sps // 2, span * sps // 2 + 1) / sps
    h = np.zeros_like(t)
    for i, ti in enumerate(t):
        if np.isclose(ti, 0.0):
            h[i] = 1.0 - beta + 4 * beta / np.pi
        elif beta > 0 and np.isclose(abs(ti), 1 / (4 * beta)):
            h[i] = (beta / np.sqrt(2)) * ((1 + 2 / np.pi) * np.sin(np.pi / (4 * beta))
                                          + (1 - 2 / np.pi) * np.cos(np.pi / (4 * beta)))
        else:
            num = np.sin(np.pi * ti * (1 - beta)) + 4 * beta * ti * np.cos(np.pi * ti * (1 + beta))
            h[i] = num / (np.pi * ti * (1 - (4 * beta * ti) ** 2))
    return h / np.sqrt(np.sum(h**2))

rng = np.random.default_rng(7)
sps, span, beta = 8, 12, 0.35
h = rrc_taps(beta, sps, span)
delay = len(h) - 1                      # total delay of two filters, in samples (a multiple of sps)

# BPSK, real baseband, sample rate = 8 x symbol rate. Each pulse has unit energy, so Eb = 1.
n_sym = 400_000
bits = rng.integers(0, 2, n_sym)
sym = 2.0 * bits - 1.0
up = np.zeros(n_sym * sps)
up[::sps] = sym
tx = np.convolve(up, h)

ebn0_db = 6.0
n0 = 1.0 / 10 ** (ebn0_db / 10)
rx = tx + rng.normal(0.0, np.sqrt(n0 / 2), tx.size)   # real white noise, two-sided PSD N0/2 per sample

# Receiver A: matched filter (the same RRC), then sample once per symbol at the peak.
mf = np.convolve(rx, h)[delay : delay + n_sym * sps : sps]
ber_mf = np.mean((mf > 0) != bits)

# Receiver B: no filter, sample the noisy waveform directly at the pulse peaks.
raw = rx[delay // 2 : delay // 2 + n_sym * sps : sps]
ber_raw = np.mean((raw > 0) != bits)

ber_theory = 0.5 * erfc(np.sqrt(10 ** (ebn0_db / 10)))
loss_db = -20 * np.log10(h.max())       # SNR lost by skipping the filter
print(f"Eb/N0 = {ebn0_db} dB, {n_sym} bits")
print(f"BER with matched filter:  {ber_mf:.2e}  (theory {ber_theory:.2e})")
print(f"BER without any filter:   {ber_raw:.2e}  (SNR penalty {loss_db:.1f} dB)")

# Eye opening against sampling offset, noise-free, for three roll-offs.
# Opening = (smallest |sample| over all symbols) / (ideal amplitude). Zero or below means the eye is closed.
print("offset (T), eye opening for beta = 0.2, 0.35, 0.5")
sps_e = 16                              # finer time grid for the eye
sym_short = sym[:4000]
offsets = range(-sps_e // 2, sps_e // 2 + 1)
eye = {}
for b in (0.2, 0.35, 0.5):
    hb = rrc_taps(b, sps_e, span)
    up_s = np.zeros(sym_short.size * sps_e)
    up_s[::sps_e] = sym_short
    y = np.convolve(np.convolve(up_s, hb), hb)
    d = len(hb) - 1
    eye[b] = []
    for off in offsets:                 # skip the first and last span symbols (filter transients)
        s = y[d + off + span * sps_e : d + off + (sym_short.size - span) * sps_e : sps_e]
        eye[b].append(np.min(s * sym_short[span : span + s.size]))
for k, off in enumerate(offsets):
    print(f"{off / sps_e:+.4f}, {eye[0.2][k]:.3f}, {eye[0.35][k]:.3f}, {eye[0.5][k]:.3f}")

Eb/N0 = 6.0 dB, 400000 bits
BER with matched filter:  2.30e-03  (theory 2.39e-03)
BER without any filter:   1.40e-01  (SNR penalty 8.2 dB)
offset (T), eye opening for beta = 0.2, 0.35, 0.5
-0.5000, -0.793, -0.485, -0.300
-0.4375, -0.621, -0.316, -0.128
-0.3750, -0.424, -0.133, 0.048
-0.3125, -0.207, 0.058, 0.226
-0.2500, 0.026, 0.253, 0.400
-0.1875, 0.268, 0.448, 0.567
-0.1250, 0.515, 0.640, 0.725
-0.0625, 0.760, 0.825, 0.870
+0.0000, 0.980, 0.995, 0.998
+0.0625, 0.757, 0.826, 0.869
+0.1250, 0.512, 0.643, 0.725
+0.1875, 0.267, 0.452, 0.567
+0.2500, 0.025, 0.257, 0.400
+0.3125, -0.207, 0.062, 0.226
+0.3750, -0.424, -0.130, 0.049
+0.4375, -0.620, -0.313, -0.128
+0.5000, -0.792, -0.484, -0.299

整合フィルタを使った場合、測定したビット誤り率（BER）は2.30 × 10⁻³で、BPSKの理論式の値は2.39 × 10⁻³です。この差は統計的なもので、400,000ビットでの誤りは約920個です。フィルタなしではBERは1.40 × 10⁻¹です。フィルタを通さないサンプルは、8倍オーバーサンプリングした全帯域の雑音を集めます。一方、信号のサンプルはパルスのピークタップ1つ分しかありません。出力されたペナルティ8.2 dBは、パルス全体のエネルギーとピークタップの2乗の比です。整合フィルタを省いたり、誤ったフィルタを使ったりする受信機は、このペナルティの一部を負うことになります。

アイパターンの読み方

アイパターンは、整合フィルタ出力を2シンボル長ずつ切り出した多数の区間を、すべてシンボルクロックにそろえて重ね描きしたものです。BPSKでは軌跡が+1と−1の間で交差し、中央に開いた領域が残ります。これがアイです。サンプリング時刻でのアイの高さは雑音マージンです。アイの幅はタイミングマージンで、雑音がなくても判定を誤るまでにサンプリング時刻がどれだけずれてよいかを表します。

ブロックの後半では、各サンプリングオフセットでの最悪ケースのアイ開口を測定します。4,000シンボルにわたって、サンプルに正しいシンボルを掛けた値の最小値をとります。したがって1.0は完全なサンプルを意味します。0以下は、雑音がまったくなくても少なくとも1つのシンボルがすでに誤っていることを意味します。

正しい時刻では、3つのアイはすべてほぼ完全に開いています。β = 0.2で0.980、β = 0.35で0.995、β = 0.5で0.998です。β = 0.2のアイの高さがわずかに低いのは、裾の減衰が遅く、12シンボルのフィルタで打ち切られる部分が多いためです。4分の1シンボルずれると、差は大きくなります。β = 0.2のアイは約0.025まで縮み、ほぼ閉じています。β = 0.35のアイは0.253、β = 0.5のアイは0.400を保ちます。ロールオフを大きくすると、帯域幅と引き換えにタイミング耐性が得られます。βを選ぶときの主なトレードオフは、フィルタ長と並んでこれです。

実際の受信機では、アイはチャネルとハードウェアについての情報も与えます。縦につぶれたアイは、雑音やISIを示します。交差点で軌跡が太くなっていれば、タイミングジッタがあります。時間とともにゆっくり回転したり閉じたりするアイは、受信機が補正していない周波数オフセットを示唆します。タイミング再生ループは、サンプリング時刻をアイの最も広い点に向けて動かすことで動作します。

全体の構成

どの実装にも当てはまる実用上の注意がいくつかあります。フィルタのスパン（ここでは12シンボル）は、阻止域の下限と遅延を決めます。各フィルタは長さの半分の遅延を加えるので、縦続接続では信号が12シンボル遅れます。アップサンプリング比は(1 + β)の帯域幅を収める余裕を残す必要があります。そのため1シンボルあたり2サンプルが実用上の最小値で、ソフトウェアでは4〜8サンプルが一般的です。ハードウェアでは、同じフィルタは通常ポリフェーズFIR構造で実装されます。その作り方はFPGA Guide DSP hubで解説しています。多くのシングルキャリア規格は、RRCフィルタとそのロールオフを規定しています。LTEや5G NRのようなマルチキャリア方式は、ISIに別の方法で対処します。レッスン8で作るサイクリックプレフィックスを使う方法です。

演習

1. 最初のブロックで、スパンを12シンボルから6シンボルに変えてください。縦続接続での新しい最悪ISIサンプルと、シンボルレートの1.0倍でのRRCスペクトラムのレベルを報告し、両方の変化を説明してください。
2. 2つ目のブロックで、受信機Aに2サンプル（4分の1シンボル）のサンプリングオフセットを与え、β = 0.2とβ = 0.5で再実行してください。BERをグラフのアイ開口と比較してください。
3. RRC送信機はそのままにして、受信側のRRCを8タップの矩形フィルタ（積分・ダンプ受信機）に置き換えてください。6 dBでのBERを測定し、それが受信機AとBの間に入る理由を説明してください。

まとめ

• 矩形パルスのスペクトラムはsincで、サイドローブは約13 dBしか下がりません。レイズドコサインパルスは、スペクトラムをシンボルレートの(1 + β)倍に収めます。
• 他のすべてのシンボル時刻でゼロを横切るパルスは、ISIがゼロです。レイズドコサインはこれを満たします。単一のルートレイズドコサインは満たしませんが、2つを縦続接続すると満たします。
• 整合フィルタは、サンプリング時刻でのSNRを最大にします。RRC送信機に対する整合フィルタは同じRRCです。これが、レイズドコサインを送信機と受信機に分ける理由です。
• 1シンボルあたり8サンプルのBPSKの試験では、整合フィルタを省くと8.2 dBのSNRを失いました。
• ロールオフを大きくするとアイが広がり、より大きなタイミング誤差に耐えられますが、その分帯域幅が必要です。

さらに読む

• H. Nyquist, “Certain topics in telegraph transmission theory”, Transactions of the AIEE, vol. 47, no. 2, 1928：ナイキスト基準の名前の由来となった論文。
• J. G. Proakis and M. Salehi, Digital Communications, 5th edition, McGraw-Hill, 2008：「Optimum Receivers for AWGN Channels」と「Digital Communication Through Band-Limited Channels」の章を参照。
• A. Goldsmith, Wireless Communications, Cambridge University Press, 2005：「Digital Modulation and Detection」の章を参照。