波形整形技術

ルートレイズドコサインフィルタ

sample.py

import numpy as np
import matplotlib.pyplot as plt

def rrc_filter(beta, span, sps):
    """
    ルートレイズドコサインフィルタ（RRC）のインパルス応答を生成
    
    Parameters
    ----------
    beta : float
        ロールオフ係数 (0 <= beta <= 1)
    span : int
        フィルタ長（シンボル数単位）
    sps : int
        1シンボルあたりのサンプル数 (samples per symbol)

    Returns
    -------
    h : ndarray
        フィルタ係数
    t : ndarray
        時間軸（単位シンボル）
    """
    # 時間軸（シンボル単位で正規化）
    N = span * sps
    t = np.arange(-N/2, N/2 + 1) / sps

    h = np.zeros_like(t, dtype=float)

    for i in range(len(t)):
        if t[i] == 0.0:
            h[i] = (1.0 - beta + 4*beta/np.pi) / np.sqrt(1.0)
        elif beta != 0 and np.isclose(abs(t[i]), 1/(4*beta)):
            h[i] = (beta/np.sqrt(2.0)) * (
                (1 + 2/np.pi) * np.sin(np.pi/(4*beta)) +
                (1 - 2/np.pi) * np.cos(np.pi/(4*beta))
            )
        else:
            numerator = (np.sin(np.pi*t[i]*(1-beta)) +
                         4*beta*t[i]*np.cos(np.pi*t[i]*(1+beta)))
            denominator = (np.pi*t[i]*(1-(4*beta*t[i])**2))
            h[i] = numerator / denominator

    # 正規化（エネルギーを1に）
    h = h / np.sqrt(np.sum(h**2))
    return h, t

# パラメータ例
beta = 0.25   # ロールオフ係数
span = 8      # フィルタ長 (シンボル数)
sps = 8       # サンプル/シンボル

h, t = rrc_filter(beta, span, sps)

# プロット
plt.plot(t, h, 'b-')
plt.title(f'Root Raised Cosine Filter (beta={beta})')
plt.xlabel('Time [symbol]')
plt.ylabel('Amplitude')
plt.grid(True)
plt.show()

処理フロー
1. 送信側
• QPSKシンボル生成
• アップサンプリング
• RRCパルス整形
• 搬送波 (fc=2.6 GHz) で変調 → 実数RF波形
2. チャネル
• （ここでは理想的にそのまま伝送、オプションで雑音を加える）
3. 受信側
• 搬送波でダウンコンバート（cos/sin を掛けて I/Q を復調）
• ローパスフィルタ（RRC マッチドフィルタ）
• シンボルタイミングでダウンサンプリング
• QPSK判定（最近傍コンステレーション）

sample.py

import numpy as np
import matplotlib.pyplot as plt
from numpy import pi

# --- パラメータ ---
fc = 2.6e9        # carrier frequency [Hz]
fs = 10e9         # sampling frequency [Hz]
Rs = 100e6        # symbol rate [symbols/sec]
sps = int(fs / Rs)  # samples per symbol
span = 8          # RRC filter span [symbols]
beta = 0.25       # roll-off
num_sym = 200     # number of symbols

# --- QPSKシンボル生成 ---
bits = np.random.randint(0, 4, size=num_sym)
sym_map = {0: 1+1j, 1: -1+1j, 2: -1-1j, 3: 1-1j}
symbols = np.array([sym_map[b] for b in bits]) / np.sqrt(2)

# --- RRCフィルタ設計 ---
def rrc_filter(beta, span, sps):
    N = span * sps
    t = np.arange(-N/2, N/2 + 1) / sps
    h = np.zeros_like(t, dtype=float)
    for i, ti in enumerate(t):
        if np.isclose(ti, 0.0):
            h[i] = 1.0 - beta + 4*beta/pi
        elif beta != 0 and np.isclose(abs(ti), 1/(4*beta)):
            h[i] = (beta/np.sqrt(2.0)) * (
                (1 + 2/pi) * np.sin(pi/(4*beta)) +
                (1 - 2/pi) * np.cos(pi/(4*beta))
            )
        else:
            num = np.sin(pi*ti*(1-beta)) + 4*beta*ti*np.cos(pi*ti*(1+beta))
            den = pi*ti*(1 - (4*beta*ti)**2)
            h[i] = num / den
    h /= np.sqrt(np.sum(h**2))
    return h, t

h, t_h = rrc_filter(beta, span, sps)

# --- アップサンプリング ---
upsampled = np.zeros(len(symbols) * sps, dtype=complex)
upsampled[::sps] = symbols

# --- パルス整形（送信側RRC） ---
s_bb = np.convolve(upsampled, h, mode="same")

# --- RF変調 ---
n = np.arange(len(s_bb))
t = n / fs
carrier = np.exp(1j * 2 * pi * fc * t)
s_rf_real = np.real(s_bb * carrier)  # 実数信号として送信

# === チャネル（ノイズ追加可） ===
rx = s_rf_real  # 今回は理想チャネル（雑音なし）

# === 受信側 ===
# I/Q ダウンコンバート
rx_I = rx * np.cos(2*pi*fc*t)
rx_Q = -rx * np.sin(2*pi*fc*t)
rx_bb = rx_I + 1j*rx_Q  # 複素ベースバンド

# RRCマッチドフィルタ
rx_matched = np.convolve(rx_bb, h, mode="same")

# シンボル判定用サンプル位置（遅延を考慮）
delay = span * sps // 2
rx_down = rx_matched[delay::sps][:num_sym]

# --- コンステレーション描画 ---
plt.figure(figsize=(6,6))
plt.scatter(np.real(rx_down), np.imag(rx_down), color='b')
plt.axhline(0, color='k', linestyle='--')
plt.axvline(0, color='k', linestyle='--')
plt.title("Received Constellation (QPSK)")
plt.xlabel("In-phase")
plt.ylabel("Quadrature")
plt.grid(True)
plt.axis('equal')
plt.show()

# --- 判定（最近傍シンボルにマッピング） ---
def decision(z):
    return (np.sign(np.real(z)) + 1j*np.sign(np.imag(z))) / np.sqrt(2)

decisions = np.array([decision(z) for z in rx_down])

# --- BER計算 ---
tx_bits = symbols  # 送信側
errors = np.sum(tx_bits != decisions)
print(f"Symbol errors: {errors}/{num_sym}, SER={errors/num_sym:.3e}")