複素ベースバンドとIQサンプル

ソフトウェア無線（SDR）はどれも複素数の列を出力します。DSPの教科書はどれも信号を複素指数関数で書きます。本レッスンではその理由を説明します。これらの数の実部と虚部が物理的に何を意味するか、ミキサが信号を搬送波とゼロ周波数の間でどう移すか、その結果にどのサンプルレートが必要かを示します。最後には、レッスン01のアップコンバータとダウンコンバータをnumpyで作り、往復の変換が正しく動くことを確かめます。

コードはnumpyのみを使います。各ブロックは単独で実行でき、1〜2秒で終わります。

1つの実信号になぜ2つの数が必要か

電波は実数の電圧です。ある瞬間の値は1つです。では、なぜ無線機はそれを2つの数で表すのでしょうか。

変調された搬送波が何であるかを見てみます。搬送波周波数fcの近くでは、どの狭帯域信号も、振幅と位相がゆっくり変化する余弦波として書けます。

x(t) = A(t) cos(2π fc t + φ(t))

情報はA(t)とφ(t)にあります。和の余弦を展開すると、90°ずれた2つの搬送波が得られます。それぞれが、ゆっくり変化する独自の振幅を持ちます。

x(t) = I(t) cos(2π fc t) − Q(t) sin(2π fc t)

with  I(t) = A(t) cos φ(t)   (in-phase)
      Q(t) = A(t) sin φ(t)   (quadrature)

同じ周波数の余弦と正弦は直交しています。多数の搬送波周期にわたって積を平均するとゼロになるため、受信機はIとQを再び分離できます。したがって、独立した2つの実信号が、1つの搬送波と同じスペクトラムの幅を共有します。1つの実数のメッセージを余弦だけに乗せた場合、同じスペクトラムを占有しながら、運べる情報は半分になります。

2つを1つの複素信号s(t) = I(t) + jQ(t)にまとめます。すると通過帯域信号は簡潔な形で書けます。

x(t) = Re{ s(t) · exp(j 2π fc t) }

複素信号s(t)は、x(t)の複素ベースバンド等価信号です。搬送波周波数そのもの以外、伝送に関するすべてを含んでいます。デジタル処理は、IQサンプルとして標本化したs(t)を扱います。s(t)はゆっくり変化するからです。3.5 GHzにある20 MHz幅の信号は、ベースバンドでは毎秒数千万サンプルで済み、数十億サンプルは要りません。

負の周波数

実信号のスペクトラムは対称です。100 kHzの実数の余弦波は、+100 kHzと−100 kHzに等しい成分を持ちます。cos(ωt) = (exp(jωt) + exp(−jωt)) / 2だからです。負の側は新しい情報を運びません。複素信号にはこの制約がありません。複素トーンexp(j2π·100 kHz·t)は、+100 kHzに1本だけのスペクトル線を持ちます。I成分は余弦、Q成分は正弦なので、フェーザは反時計回りに回転します。−100 kHzのトーンは時計回りに回転します。

ベースバンドでは、負の周波数は単純な物理的意味を持ちます。搬送波より下にある成分です。−260 kHzの複素ベースバンドトーンを3 MHzの搬送波にアップミキシングすると、空中では2.74 MHzに現れます。ベースバンドのスペクトラムは対称である必要がないため、搬送波の上側と下側で異なる情報を運べます。

以下のコードは、1 MS/sの複素ベースバンド信号を2つのトーンで作ります。振幅1の+100 kHzと、振幅0.5の−260 kHzです。複素信号のFFTで最も強いビンを求め、次に実部だけについて同じことをします。

import numpy as np

fs = 1e6                       # complex sample rate, Hz
n = np.arange(4000)            # 4 ms of samples
t = n / fs

# Two complex tones: one above the centre frequency, one below it
s = 1.0 * np.exp(2j * np.pi * 100e3 * t) + 0.5 * np.exp(-2j * np.pi * 260e3 * t)
s = s.astype(np.complex64)
print(f"dtype {s.dtype}, {s.size} samples at {fs / 1e6:.0f} MS/s")
print(f"first sample: I = {s[0].real:.3f}, Q = {s[0].imag:.3f}")

def peak_frequencies(x: np.ndarray, fs: float, count: int) -> np.ndarray:
    """Frequencies of the strongest FFT bins, sorted."""
    spec = np.fft.fftshift(np.abs(np.fft.fft(x)))
    freqs = np.fft.fftshift(np.fft.fftfreq(x.size, d=1 / fs))
    return np.sort(freqs[np.argsort(spec)[-count:]])

print("complex signal, peaks at (kHz):", peak_frequencies(s, fs, 2) / 1e3)
print("real part only, peaks at (kHz):", peak_frequencies(s.real, fs, 4) / 1e3)

dtype complex64, 4000 samples at 1 MS/s
first sample: I = 1.500, Q = 0.000
complex signal, peaks at (kHz): [-260.  100.]
real part only, peaks at (kHz): [-260. -100.  100.  260.]

複素信号には、−260 kHzと+100 kHzのちょうど2本の線が現れます。Qを捨てると、実部には4本が現れます。各トーンに鏡像ができ、エネルギーが搬送波のどちら側から来たかを区別できなくなります。受信機にIとQの両方が必要なのはこのためです。サンプルはnumpy.complex64で、1サンプルあたり32 ビットの浮動小数点数2個です。これはGNU Radioのgr_complex型と同じ配置です。gr_complexの実体はstd::complex<float>です。

アップコンバートとダウンコンバート

レッスン01では、各無線機にミキサがあることを示しました。数学的には、ミキサは乗算器です。アップコンバータはI cos(2π fc t) − Q sin(2π fc t)を作ります。これは上の式そのものです。ダウンコンバータは、受信した実信号に、同じ局部発振器（LO）からの余弦と負の正弦をそれぞれ掛け、各積をローパスフィルタに通します。

なぜこれで動くのでしょうか。x(t)に2 exp(−j2π fc t)を掛けます。これは、余弦の系統と、j倍した負の正弦の系統の和です。x(t)の簡潔な形を使うと、積はs(t)と第2の項の和になります。第2の項は、s(t)の複素共役を−2fcへ周波数シフトしたものです。ローパスフィルタが第2の項を除去し、s(t)が残ります。係数2は元の振幅を戻すためのものです。

次のブロックは、アナログの世界を16 MS/sで模擬し、2トーン信号を3 MHzの搬送波に乗せてから、ダウンコンバータを動かします。フィルタは129タップの窓付きsinc FIRで、カットオフは1 MHz付近です。フィルタ後に16サンプルごとに1つを残し、レートを1 MS/sに戻します。

import numpy as np

fs_rf = 16e6                   # sample rate of the simulated "analogue" world, Hz
fc = 3e6                       # carrier frequency, Hz
decim = 16                     # 16 MS/s / 16 = 1 MS/s complex baseband
t = np.arange(64000) / fs_rf   # 4 ms

def baseband(t: np.ndarray) -> np.ndarray:
    return 1.0 * np.exp(2j * np.pi * 100e3 * t) + 0.5 * np.exp(-2j * np.pi * 260e3 * t)

s = baseband(t)

# Up-conversion: one real signal on the air, x = I cos(wt) - Q sin(wt)
x = s.real * np.cos(2 * np.pi * fc * t) - s.imag * np.sin(2 * np.pi * fc * t)
spec = np.abs(np.fft.rfft(x))
f_axis = np.fft.rfftfreq(x.size, d=1 / fs_rf)
print("passband peaks (MHz):", np.sort(f_axis[np.argsort(spec)[-2:]]) / 1e6)

# Down-conversion: multiply by 2 exp(-j w t), then low-pass filter and decimate
mixed = x * 2 * np.exp(-2j * np.pi * fc * t)
taps = np.sinc(2 * 1e6 / fs_rf * np.arange(-64, 65)) * np.hamming(129)
taps /= taps.sum()             # unity gain at 0 Hz, cut-off near 1 MHz
filtered = np.convolve(mixed, taps, mode="same")
r = filtered[::decim]          # 1 MS/s complex baseband
ref = s[::decim]

core = slice(10, -10)          # skip the filter's edge transient
err = r[core] - ref[core]
print(f"recovered {r.size} samples at {fs_rf / decim / 1e6:.0f} MS/s")
print(f"error power relative to signal: "
      f"{10 * np.log10(np.mean(np.abs(err)**2) / np.mean(np.abs(ref[core])**2)):.1f} dB")

passband peaks (MHz): [2.74 3.1 ]
recovered 4000 samples at 1 MS/s
error power relative to signal: -50.6 dB

空中では、2つのトーンは3.1 MHzと2.74 MHzにあります。搬送波の100 kHz上と260 kHz下で、ベースバンドでの符号が示したとおりです。ダウンコンバート後、復元したサンプルは元のサンプルと一致し、誤差は信号より50.6 dB低くなっています。筆者らはこの残差の発生源を確認しました。6 MHzの項は、これよりはるかに低く抑えられています。残差のほぼすべては、100 kHzと260 kHzでのフィルタのわずかな利得誤差です。そのため、通過帯域がより平坦なフィルタを使えば、残差はさらに下がります。

実際のハードウェアでは、このシミュレーションが省いた誤差が加わります。IとQの系統の利得が完全に等しく、位相差がちょうど90°になることはありません。このため、各トーンの弱い鏡像が生じます。LOがミキサ出力に漏れると、0 Hzにスパイクが加わります。さらに、送信機と受信機の発振器が完全に同じ周波数で動くことはないため、復元信号はゆっくり回転します。この最後の効果が搬送波周波数オフセットで、実際の受信機が最初に推定して除去しなければならないものです。

サンプルレートとナイキスト限界

どれだけ速く標本化すればよいでしょうか。最高周波数がBの実信号は、毎秒2Bより速く取ったサンプルで完全に決まります。これが標本化定理で、2Bがナイキストレートです。シャノンは、雑音下の通信に関する1949年の論文の定理1としてこれを述べました。Wより高い周波数成分を持たない信号は、1/(2W)秒間隔のサンプルで完全に決まる、というものです。シャノンは、この間隔の重要性を認識した功績を、電信に関するナイキストの1928年の研究に帰しています。

複素信号では、この条件が緩みます。各複素サンプルは2つの実数で、スペクトラムは対称である必要がありません。複素サンプルレートfsは、−fs/2から+fs/2までの帯域、つまり幅fs Hzをカバーします。したがって、20 MHz幅の信号には少なくとも20 MS/sの複素サンプルレートが必要です。実際には、フィルタのロールオフの余裕を見てもう少し高くします。毎秒の実数値の総数は4000万で、ナイキストレートでの実数標本化と同じです。何かを得するわけではありません。データが処理に適した形に並ぶだけです。

SDRを中心周波数fc、サンプルレートfsに設定すると、fc − fs/2からfc + fs/2までの帯域が見えます。ベースバンドの0 Hzが中心周波数に当たり、負のベースバンド周波数は帯域の下半分です。

エイリアシング

標本化された信号は、fsの整数倍だけ異なる2つの周波数を区別できません。1 MS/sでは、+640 kHzのトーンは−360 kHzのトーンとまったく同じサンプルを生みます。640 kHz − 1000 kHz = −360 kHzだからです。±fs/2の窓の外にあるエネルギーがADCに届くと、窓の内側に折り返し、そこにあるものに重なります。これがエイリアシングです。対策は、ADCの前に置くアナログのアンチエイリアシングフィルタしかありません。標本化した後では、損傷を元に戻せないからです。

最後のブロックは、前の2つのトーンに、+640 kHz・振幅0.25の第3のトーンと複素ガウス雑音を加えます。そして平均化したスペクトラムを計算します。条件は、50サンプルのフレーム400個、Hann窓、1ビンあたり20 kHzです。グラフはこのスペクトラムをプロットしたものです。

import numpy as np

rng = np.random.default_rng(1)
fs = 1e6                         # complex sample rate, Hz
nfft, frames = 50, 400           # 20 kHz bins, 400 frames averaged
n = np.arange(nfft * frames)

# +640 kHz lies outside the -500 to +500 kHz range that 1 MS/s complex sampling can hold
alias = np.exp(2j * np.pi * 640e3 * n / fs)
same = np.allclose(alias, np.exp(-2j * np.pi * 360e3 * n / fs))
print(f"+640 kHz and -360 kHz give identical samples at 1 MS/s: {same}")

s = (1.0 * np.exp(2j * np.pi * 100e3 * n / fs)
     + 0.5 * np.exp(-2j * np.pi * 260e3 * n / fs)
     + 0.25 * alias)
noise = (rng.standard_normal(n.size) + 1j * rng.standard_normal(n.size)) * np.sqrt(1e-3 / 2)
x = (s + noise).astype(np.complex64)

# Averaged periodogram with a Hann window, normalised to the strongest bin
w = np.hanning(nfft)
blocks = x.reshape(frames, nfft) * w
psd = np.fft.fftshift(np.mean(np.abs(np.fft.fft(blocks, axis=1))**2, axis=0))
freqs = np.fft.fftshift(np.fft.fftfreq(nfft, d=1 / fs))
psd_db = 10 * np.log10(psd / psd.max())

for k in range(1, nfft - 1):
    if psd_db[k] > -30 and psd_db[k] >= psd_db[k - 1] and psd_db[k] >= psd_db[k + 1]:
        print(f"peak at {freqs[k] / 1e3:5.0f} kHz: {psd_db[k]:6.1f} dB")
print(f"median noise floor: {np.median(psd_db):.1f} dB")
print("chart:", ", ".join(f"[{f / 1e3:.0f}, {p:.1f}]" for f, p in zip(freqs, psd_db)))

+640 kHz and -360 kHz give identical samples at 1 MS/s: True
peak at  -360 kHz:  -12.0 dB
peak at  -260 kHz:   -6.0 dB
peak at   100 kHz:    0.0 dB
median noise floor: -44.9 dB
chart: [-500, -45.1], [-480, -45.2], ... , [480, -45.2]

chart:の行全体はコードが出力します。3つのピークは0 dB、−6.0 dB、−12.0 dBにあり、振幅比1、0.5、0.25をデシベルで表した値です。+640 kHzで入れたトーンは−360 kHzに現れ、偽物であることを示す手がかりは何もありません。各ピークが3ビン幅なのは、Hann窓が周波数分解能と引き換えにビン間の漏れを小さくするためです。−45 dB付近のフロアは、50ビンに広がった雑音です。

演習

1. 最初のブロックの−260 kHzのトーンを+260 kHzに変えてください。実部のみのスペクトラムには何が現れるでしょうか。また、なぜ2つの場合を区別できなくなるのでしょうか。
2. ミキシングのブロックで、フィルタを33タップに短くし、次に257タップに長くしてください。そのたびに誤差電力をdBで記録してください。次に、ミキシングせずにsを直接フィルタに通して比較してください。誤差を決めているのは6 MHzの項でしょうか、それとも通過帯域でしょうか。
3. 受信機の指数関数でfc + 200 Hzを使い、ダウンコンバータに小さなLO周波数誤差を加えてください。r[core] / ref[core]の位相を時間に対してプロットしてください。どれだけの速さで回転するでしょうか。
4. Q系統に1.05を掛けて、受信機にIQインバランスを加えてください。スペクトラムで+100 kHzのトーンの鏡像を見つけ、トーンに対するレベルを測ってください。
5. 2.45 GHzにある15 MHz幅の信号を受信したいとします。この信号を収められる最小の複素サンプルレートはいくつでしょうか。アンチエイリアシングフィルタは、どの周波数帯域を除去する必要があるでしょうか。

まとめ

• 狭帯域の通過帯域信号はx(t) = I(t) cos(2π fc t) − Q(t) sin(2π fc t)です。複素信号s(t) = I(t) + jQ(t)はその等価ベースバンド信号で、すべての情報を持っています。
• 複素信号は、正の周波数と負の周波数で異なる内容を持てます。ベースバンドでは、負の周波数は搬送波より下の成分です。
• アップコンバートでは搬送波を掛けます。ダウンコンバートでは2 exp(−j2π fc t)を掛け、ローパスフィルタに通し、デシメーションします。上のnumpyでの往復では、信号より50.6 dB低い誤差で信号を復元しました。
• 複素サンプルレートfsは、−fs/2から+fs/2までの帯域を収めます。その外にあるものはすべて帯域内にエイリアスするため、ADCの前でフィルタにより除去する必要があります。

さらに読む

• H. Nyquist, “Certain topics in telegraph transmission theory”, 「Transactions of the AIEE」, vol. 47, no. 2, pp. 617-644, 1928: シャノンがナイキスト間隔の功績を認めている電信の論文です。
• C. E. Shannon, “Communication in the presence of noise”, 「Proceedings of the IRE」, vol. 37, no. 1, pp. 10-21, 1949: 定理1が本レッスンで使った標本化定理です。定理2はレッスン04で使う通信路容量の式です。
• J. G. Proakis and M. Salehi, 「Digital Communications」, 5th edition, McGraw-Hill, 2008: 第2章「Deterministic and Random Signal Analysis」は、本レッスンの背景にある信号理論の教科書的な解説です。
• D. Tse and P. Viswanath, 「Fundamentals of Wireless Communication」, Cambridge University Press, 2005: 第2章「The wireless channel」は、これらのベースバンド信号が通過するチャネルをモデル化しています。