多重化技術 — FDM・TDM・CDMによる通信路の共有

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import butter, filtfilt

# パラメータ
fs = 10000          # サンプリング周波数 [Hz]
T = 0.5             # 信号長 [s]
t = np.arange(0, T, 1/fs)
N = len(t)

# 3チャンネルのメッセージ信号（異なる周波数）
m1 = np.cos(2 * np.pi * 10 * t)  # 10 Hz 正弦波
m2 = 0.5 * (np.sign(np.cos(2 * np.pi * 8 * t)))  # 8 Hz 矩形波（近似）
m3 = np.cos(2 * np.pi * 15 * t) + 0.3 * np.cos(2 * np.pi * 5 * t)  # 複合波

# 搬送波周波数（十分に離す: ガードバンドあり）
fc1, fc2, fc3 = 500, 800, 1100  # Hz
guard_band = 100  # Hz

# FDM変調
s1 = m1 * np.cos(2 * np.pi * fc1 * t)
s2 = m2 * np.cos(2 * np.pi * fc2 * t)
s3 = m3 * np.cos(2 * np.pi * fc3 * t)
s_fdm = s1 + s2 + s3

# スペクトル
freqs = np.fft.rfftfreq(N, 1/fs)
S_fdm = np.abs(np.fft.rfft(s_fdm)) / N * 2

# FDM復調（帯域通過フィルタ + 復調）
def bandpass_demod(signal, fc, bw, fs, t):
    """帯域通過フィルタ → 搬送波乗算 → ローパスフィルタ"""
    # 帯域通過フィルタ
    low = (fc - bw) / (fs/2)
    high = (fc + bw) / (fs/2)
    b_bp, a_bp = butter(4, [low, high], btype='band')
    filtered = filtfilt(b_bp, a_bp, signal)

    # 搬送波乗算（同期検波）
    demod = 2 * filtered * np.cos(2 * np.pi * fc * t)

    # ローパスフィルタ
    b_lp, a_lp = butter(4, bw / (fs/2), btype='low')
    recovered = filtfilt(b_lp, a_lp, demod)
    return recovered

bw = 80  # 帯域通過フィルタの半幅 [Hz]
m1_recv = bandpass_demod(s_fdm, fc1, bw, fs, t)
m2_recv = bandpass_demod(s_fdm, fc2, bw, fs, t)
m3_recv = bandpass_demod(s_fdm, fc3, bw, fs, t)

# プロット
fig, axes = plt.subplots(4, 2, figsize=(14, 14))

# 左列: 送信側
axes[0, 0].plot(t[:1000], m1[:1000], 'b-', linewidth=1.5)
axes[0, 0].set_title('Ch1: Message signal (10 Hz sine)')
axes[0, 0].set_ylabel('Amplitude')
axes[0, 0].grid(True, alpha=0.3)

axes[1, 0].plot(t[:1000], m2[:1000], 'g-', linewidth=1.5)
axes[1, 0].set_title('Ch2: Message signal (8 Hz square)')
axes[1, 0].set_ylabel('Amplitude')
axes[1, 0].grid(True, alpha=0.3)

axes[2, 0].plot(t[:1000], m3[:1000], 'r-', linewidth=1.5)
axes[2, 0].set_title('Ch3: Message signal (composite)')
axes[2, 0].set_ylabel('Amplitude')
axes[2, 0].grid(True, alpha=0.3)

# FDMスペクトル
axes[3, 0].plot(freqs, 20*np.log10(S_fdm + 1e-10), 'k-', linewidth=0.8)
axes[3, 0].set_xlim([0, 1500])
axes[3, 0].set_ylim([-60, 0])
axes[3, 0].set_title('FDM Composite Spectrum')
axes[3, 0].set_xlabel('Frequency [Hz]')
axes[3, 0].set_ylabel('Magnitude [dB]')
for fc, label in zip([fc1, fc2, fc3], ['Ch1', 'Ch2', 'Ch3']):
    axes[3, 0].axvline(x=fc, color='r', linestyle='--', alpha=0.5)
    axes[3, 0].text(fc, -5, label, ha='center', fontsize=9)
axes[3, 0].grid(True, alpha=0.3)

# 右列: 受信側（復調後）
axes[0, 1].plot(t[:1000], m1_recv[:1000], 'b-', linewidth=1.5, label='Recovered')
axes[0, 1].plot(t[:1000], m1[:1000], 'k--', linewidth=1, alpha=0.5, label='Original')
axes[0, 1].set_title('Ch1: Recovered')
axes[0, 1].set_ylabel('Amplitude')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)

axes[1, 1].plot(t[:1000], m2_recv[:1000], 'g-', linewidth=1.5, label='Recovered')
axes[1, 1].plot(t[:1000], m2[:1000], 'k--', linewidth=1, alpha=0.5, label='Original')
axes[1, 1].set_title('Ch2: Recovered')
axes[1, 1].set_ylabel('Amplitude')
axes[1, 1].legend()
axes[1, 1].grid(True, alpha=0.3)

axes[2, 1].plot(t[:1000], m3_recv[:1000], 'r-', linewidth=1.5, label='Recovered')
axes[2, 1].plot(t[:1000], m3[:1000], 'k--', linewidth=1, alpha=0.5, label='Original')
axes[2, 1].set_title('Ch3: Recovered')
axes[2, 1].set_ylabel('Amplitude')
axes[2, 1].legend()
axes[2, 1].grid(True, alpha=0.3)

# 復調精度
for idx, (orig, recv, name) in enumerate(zip([m1, m2, m3], [m1_recv, m2_recv, m3_recv], ['Ch1', 'Ch2', 'Ch3'])):
    mid = slice(1000, 4000)
    corr = np.corrcoef(orig[mid], recv[mid])[0, 1]
    axes[idx, 1].text(0.02, 0.05, f'Correlation: {corr:.4f}', transform=axes[idx, 1].transAxes, fontsize=9)

axes[3, 1].axis('off')
axes[3, 1].text(0.1, 0.5, 'FDM: Each channel occupies\na separate frequency band.\nBPF + coherent detection\nrecovers the original signals.',
                transform=axes[3, 1].transAxes, fontsize=11, va='center')

plt.tight_layout()
plt.savefig('fdm_demo.png', dpi=150, bbox_inches='tight')
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# パラメータ
N_channels = 4       # チャンネル数
N_frames = 20        # フレーム数
samples_per_slot = 50 # スロットあたりサンプル数
guard_samples = 5     # ガードタイム（サンプル数）

slot_total = samples_per_slot + guard_samples
frame_length = N_channels * slot_total

# 各チャンネルのデータ信号
np.random.seed(42)
channels = []
channel_labels = ['Voice A', 'Voice B', 'Data C', 'Video D']
channel_colors = ['blue', 'green', 'red', 'purple']

for ch in range(N_channels):
    freq = 3 + ch * 2  # 異なる周波数
    t_ch = np.arange(N_frames * samples_per_slot) / 1000
    signal = np.sin(2 * np.pi * freq * t_ch) + 0.3 * np.random.randn(len(t_ch))
    channels.append(signal)

# TDM多重化
tdm_signal = np.zeros(N_frames * frame_length)
for frame in range(N_frames):
    for ch in range(N_channels):
        src_start = frame * samples_per_slot
        src_end = src_start + samples_per_slot
        dst_start = frame * frame_length + ch * slot_total
        dst_end = dst_start + samples_per_slot
        tdm_signal[dst_start:dst_end] = channels[ch][src_start:src_end]
        # ガードタイムはゼロ（デフォルト）

# TDM分離
recovered = [np.zeros(N_frames * samples_per_slot) for _ in range(N_channels)]
for frame in range(N_frames):
    for ch in range(N_channels):
        src_start = frame * frame_length + ch * slot_total
        src_end = src_start + samples_per_slot
        dst_start = frame * samples_per_slot
        dst_end = dst_start + samples_per_slot
        recovered[ch][dst_start:dst_end] = tdm_signal[src_start:src_end]

# プロット
fig, axes = plt.subplots(3, 1, figsize=(14, 10))

# 各チャンネルの元信号（一部）
for ch in range(N_channels):
    t_plot = np.arange(3 * samples_per_slot) / 1000
    axes[0].plot(t_plot, channels[ch][:3*samples_per_slot],
                color=channel_colors[ch], linewidth=1.5, label=channel_labels[ch])
axes[0].set_title('Original Channel Signals (first 3 frames)')
axes[0].set_ylabel('Amplitude')
axes[0].legend(loc='upper right')
axes[0].grid(True, alpha=0.3)

# TDM多重化信号（3フレーム分）
t_tdm = np.arange(3 * frame_length)
tdm_3frames = tdm_signal[:3 * frame_length]
# カラーリング
for frame in range(3):
    for ch in range(N_channels):
        start = frame * frame_length + ch * slot_total
        end = start + samples_per_slot
        t_slot = np.arange(start, end)
        axes[1].fill_between(t_slot, -2, 2, alpha=0.1, color=channel_colors[ch])
        axes[1].plot(t_slot, tdm_signal[start:end],
                    color=channel_colors[ch], linewidth=1)
    # ガードタイムを表示
    for ch in range(N_channels):
        guard_start = frame * frame_length + ch * slot_total + samples_per_slot
        guard_end = guard_start + guard_samples
        axes[1].axvspan(guard_start, guard_end, alpha=0.3, color='gray')

axes[1].set_title(f'TDM Signal (3 frames, {N_channels} slots/frame, guard time = {guard_samples} samples)')
axes[1].set_ylabel('Amplitude')
axes[1].set_xlabel('Sample index')
axes[1].grid(True, alpha=0.3)

# 復元信号
for ch in range(N_channels):
    t_recv = np.arange(3 * samples_per_slot) / 1000
    axes[2].plot(t_recv, recovered[ch][:3*samples_per_slot],
                color=channel_colors[ch], linewidth=1.5, label=channel_labels[ch])
axes[2].set_title('Recovered Channel Signals')
axes[2].set_ylabel('Amplitude')
axes[2].set_xlabel('Time [s]')
axes[2].legend(loc='upper right')
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('tdm_demo.png', dpi=150, bbox_inches='tight')
plt.show()

# 復元精度
print("=== TDM復元精度 ===")
for ch in range(N_channels):
    mse = np.mean((channels[ch] - recovered[ch])**2)
    print(f"  {channel_labels[ch]}: MSE = {mse:.2e}")

import numpy as np
import matplotlib.pyplot as plt

def hadamard_matrix(n):
    """n次アダマール行列を生成（nは2のべき乗）"""
    if n == 1:
        return np.array([[1]])
    H_half = hadamard_matrix(n // 2)
    return np.block([[H_half, H_half], [H_half, -H_half]])

# パラメータ
N = 8  # ウォルシュ符号長（= 最大ユーザー数）
K = 4  # アクティブユーザー数
N_bits = 10  # 各ユーザーのビット数

np.random.seed(42)

# ウォルシュ符号
H = hadamard_matrix(N)
walsh_codes = H[1:K+1]  # 行0はオール1なので避ける

# 各ユーザーのデータ
data = 2 * np.random.randint(0, 2, (K, N_bits)) - 1

# CDM送信: 各ビットをウォルシュ符号で拡散して加算
cdm_signal = np.zeros((N_bits, N))
for bit_idx in range(N_bits):
    for user in range(K):
        cdm_signal[bit_idx] += data[user, bit_idx] * walsh_codes[user]

# AWGN雑音を付加
snr_db = 10
noise_power = K * N * 10**(-snr_db/10)  # K*Nは信号電力の目安
noise = np.sqrt(noise_power/2) * np.random.randn(N_bits, N)
cdm_received = cdm_signal + noise

# CDM復調: ウォルシュ符号との相関
decoded = np.zeros((K, N_bits), dtype=int)
correlation_values = np.zeros((K, N_bits))
for user in range(K):
    for bit_idx in range(N_bits):
        corr = np.dot(cdm_received[bit_idx], walsh_codes[user]) / N
        correlation_values[user, bit_idx] = corr
        decoded[user, bit_idx] = 1 if corr > 0 else -1

# 結果表示
fig, axes = plt.subplots(3, 1, figsize=(14, 10))

# ウォルシュ符号の直交性マトリクス
cross_corr = np.zeros((K, K))
for i in range(K):
    for j in range(K):
        cross_corr[i, j] = np.dot(walsh_codes[i], walsh_codes[j]) / N

im = axes[0].imshow(cross_corr, cmap='RdBu_r', vmin=-1, vmax=1, aspect='auto')
axes[0].set_title(f'Walsh Code Cross-Correlation Matrix ({K} users, length {N})')
axes[0].set_xlabel('User index')
axes[0].set_ylabel('User index')
for i in range(K):
    for j in range(K):
        axes[0].text(j, i, f'{cross_corr[i,j]:.1f}', ha='center', va='center', fontsize=12)
plt.colorbar(im, ax=axes[0])

# CDM合成信号の1ビット分
bit_show = 0
chip_indices = np.arange(N)
bar_width = 0.35
for user in range(K):
    axes[1].bar(chip_indices + user*bar_width/K, data[user, bit_show] * walsh_codes[user],
               width=bar_width/K, label=f'User {user+1} (d={data[user,bit_show]:+d})',
               alpha=0.7)
axes[1].bar(chip_indices + 0.5, cdm_signal[bit_show], width=0.15, color='black',
           label='CDM sum', alpha=0.8)
axes[1].set_title(f'CDM Signal Composition (bit index = {bit_show})')
axes[1].set_xlabel('Chip index')
axes[1].set_ylabel('Amplitude')
axes[1].legend(loc='upper right', fontsize=8)
axes[1].grid(True, alpha=0.3)

# 復調結果の比較
user_labels = [f'User {i+1}' for i in range(K)]
x_pos = np.arange(N_bits)
bar_w = 0.8 / K
for user in range(K):
    offset = (user - K/2 + 0.5) * bar_w
    colors = ['green' if d == o else 'red' for d, o in zip(decoded[user], data[user])]
    axes[2].bar(x_pos + offset, correlation_values[user], width=bar_w,
               color=colors, alpha=0.6, edgecolor='black', linewidth=0.5)
axes[2].axhline(y=0, color='k', linewidth=0.5)
axes[2].set_title(f'Correlation Values after Despreading (SNR = {snr_db} dB, green=correct, red=error)')
axes[2].set_xlabel('Bit index')
axes[2].set_ylabel('Correlation')
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('cdm_orthogonality.png', dpi=150, bbox_inches='tight')
plt.show()

# BER計算
total_errors = np.sum(decoded != data)
total_bits = K * N_bits
print(f"=== CDM復調結果 ===")
print(f"ユーザー数: {K}, 符号長: {N}, SNR: {snr_db} dB")
print(f"総ビット数: {total_bits}, エラー数: {total_errors}, BER: {total_errors/total_bits:.4f}")
for user in range(K):
    errors = np.sum(decoded[user] != data[user])
    print(f"  User {user+1}: 送信 {data[user]}, 復調 {decoded[user]}, エラー: {errors}")