ASK/FSK変調の理論 — デジタル振幅・周波数変調の原理と実装

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
fc = 10.0        # 搬送波周波数 [Hz]
Rs = 1.0         # シンボルレート [symbol/s]
Ts = 1.0 / Rs    # シンボル周期 [s]
fs = 1000.0      # サンプリング周波数 [Hz]
N_bits = 10      # ビット数

# ランダムビット列
np.random.seed(42)
bits = np.random.randint(0, 2, N_bits)

# 時間軸
t_total = N_bits * Ts
t = np.arange(0, t_total, 1/fs)

# ベースバンド信号（矩形パルス列）
baseband = np.zeros(len(t))
for i, b in enumerate(bits):
    idx = (t >= i * Ts) & (t < (i + 1) * Ts)
    baseband[idx] = b

# OOK変調信号
carrier = np.cos(2 * np.pi * fc * t)
ook_signal = baseband * carrier

# パワースペクトル密度
from numpy.fft import fft, fftfreq
N_fft = len(ook_signal)
freq = fftfreq(N_fft, 1/fs)
spectrum = np.abs(fft(ook_signal)) ** 2 / N_fft

fig, axes = plt.subplots(3, 1, figsize=(12, 9))

# ビット列とベースバンド
axes[0].step(np.arange(N_bits), bits, where='mid', linewidth=2, color='tab:blue')
axes[0].set_ylabel('Bit value')
axes[0].set_title('Data Bits')
axes[0].set_xlim(-0.5, N_bits - 0.5)
axes[0].set_ylim(-0.2, 1.4)
axes[0].grid(True, alpha=0.3)
for i, b in enumerate(bits):
    axes[0].text(i, b + 0.1, str(b), ha='center', fontsize=10)

# OOK変調波形
axes[1].plot(t, ook_signal, linewidth=0.8, color='tab:orange')
axes[1].set_ylabel('Amplitude')
axes[1].set_title('OOK Modulated Signal')
axes[1].set_xlim(0, t_total)
axes[1].grid(True, alpha=0.3)

# パワースペクトル
mask = (freq >= 0) & (freq <= 30)
axes[2].plot(freq[mask], 10 * np.log10(spectrum[mask] / np.max(spectrum[mask]) + 1e-12),
             linewidth=1, color='tab:green')
axes[2].set_xlabel('Frequency [Hz]')
axes[2].set_ylabel('PSD [dB]')
axes[2].set_title('Power Spectral Density of OOK Signal')
axes[2].set_ylim(-60, 5)
axes[2].grid(True, alpha=0.3)
axes[2].axvline(x=fc, color='r', linestyle='--', alpha=0.5, label=f'fc = {fc} Hz')
axes[2].legend()

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

import numpy as np
import matplotlib.pyplot as plt

fig, axes = plt.subplots(2, 3, figsize=(15, 8))

fc = 10.0
fs = 1000.0
N_symbols = 8

for idx, M in enumerate([2, 4, 8]):
    m = int(np.log2(M))

    # コンスタレーション点
    constellation = np.arange(-(M-1), M, 2)  # {-(M-1), ..., -1, 1, ..., (M-1)}

    # ランダムシンボル
    np.random.seed(idx)
    symbol_indices = np.random.randint(0, M, N_symbols)
    symbols = constellation[symbol_indices]

    # コンスタレーション図
    ax = axes[0, idx]
    ax.scatter(constellation, np.zeros(M), s=100, c='tab:blue', zorder=3)
    for i, c in enumerate(constellation):
        bits_label = format(i, f'0{m}b')
        ax.annotate(bits_label, (c, 0.05), ha='center', fontsize=9)
    ax.axhline(y=0, color='gray', linewidth=0.5)
    ax.axvline(x=0, color='gray', linewidth=0.5)
    ax.set_xlim(-(M+1), M+1)
    ax.set_ylim(-0.5, 0.5)
    ax.set_title(f'{M}-ASK Constellation')
    ax.set_xlabel('In-phase')
    ax.grid(True, alpha=0.3)

    # 変調波形
    Ts = 1.0
    t_total = N_symbols * Ts
    t = np.arange(0, t_total, 1/fs)
    signal = np.zeros(len(t))
    for i, s in enumerate(symbols):
        mask = (t >= i * Ts) & (t < (i + 1) * Ts)
        signal[mask] = s * np.cos(2 * np.pi * fc * t[mask])

    ax = axes[1, idx]
    ax.plot(t, signal, linewidth=0.8, color='tab:orange')
    ax.set_xlabel('Time [s]')
    ax.set_ylabel('Amplitude')
    ax.set_title(f'{M}-ASK Modulated Signal')
    ax.set_xlim(0, t_total)
    ax.grid(True, alpha=0.3)

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

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
fc = 10.0         # 搬送波中心周波数 [Hz]
delta_f = 2.0     # 周波数偏移 [Hz]（h = delta_f * Ts = 2.0）
f0 = fc - delta_f / 2  # ビット0の周波数
f1 = fc + delta_f / 2  # ビット1の周波数
Rs = 1.0
Ts = 1.0 / Rs
fs = 1000.0
N_bits = 10

# ランダムビット列
np.random.seed(42)
bits = np.random.randint(0, 2, N_bits)

# 時間軸
t_total = N_bits * Ts
t = np.arange(0, t_total, 1/fs)

# 非連続位相FSK
fsk_signal = np.zeros(len(t))
for i, b in enumerate(bits):
    mask = (t >= i * Ts) & (t < (i + 1) * Ts)
    freq = f1 if b == 1 else f0
    fsk_signal[mask] = np.cos(2 * np.pi * freq * t[mask])

# CPFSK（連続位相）
cpfsk_signal = np.zeros(len(t))
phase = 0.0
dt = 1.0 / fs
for i in range(len(t)):
    symbol_idx = min(int(t[i] / Ts), N_bits - 1)
    freq = f1 if bits[symbol_idx] == 1 else f0
    cpfsk_signal[i] = np.cos(2 * np.pi * fc * t[i] + phase)
    phase += 2 * np.pi * (freq - fc) * dt

# パワースペクトル
from numpy.fft import fft, fftfreq
N_fft = len(fsk_signal)
freq_axis = fftfreq(N_fft, 1/fs)
spectrum_disc = np.abs(fft(fsk_signal)) ** 2 / N_fft
spectrum_cont = np.abs(fft(cpfsk_signal)) ** 2 / N_fft

fig, axes = plt.subplots(3, 1, figsize=(12, 10))

# ビット列
axes[0].step(np.arange(N_bits), bits, where='mid', linewidth=2, color='tab:blue')
axes[0].set_ylabel('Bit value')
axes[0].set_title('Data Bits')
axes[0].set_xlim(-0.5, N_bits - 0.5)
axes[0].set_ylim(-0.2, 1.4)
axes[0].grid(True, alpha=0.3)
for i, b in enumerate(bits):
    axes[0].text(i, b + 0.1, str(b), ha='center', fontsize=10)

# CPFSK波形
axes[1].plot(t, cpfsk_signal, linewidth=0.8, color='tab:orange')
axes[1].set_ylabel('Amplitude')
axes[1].set_title(f'CPFSK Signal (fc={fc} Hz, Δf={delta_f} Hz, h={delta_f*Ts:.1f})')
axes[1].set_xlim(0, t_total)
axes[1].grid(True, alpha=0.3)

# シンボル境界を描画
for i in range(1, N_bits):
    axes[1].axvline(x=i * Ts, color='gray', linestyle=':', alpha=0.5)

# パワースペクトル比較
mask = (freq_axis >= 0) & (freq_axis <= 25)
axes[2].plot(freq_axis[mask],
             10 * np.log10(spectrum_disc[mask] / np.max(spectrum_disc[mask]) + 1e-12),
             linewidth=1, color='tab:blue', label='Discontinuous FSK', alpha=0.7)
axes[2].plot(freq_axis[mask],
             10 * np.log10(spectrum_cont[mask] / np.max(spectrum_cont[mask]) + 1e-12),
             linewidth=1, color='tab:red', label='CPFSK', alpha=0.7)
axes[2].set_xlabel('Frequency [Hz]')
axes[2].set_ylabel('PSD [dB]')
axes[2].set_title('Power Spectral Density Comparison')
axes[2].set_ylim(-60, 5)
axes[2].grid(True, alpha=0.3)
axes[2].axvline(x=f0, color='green', linestyle='--', alpha=0.5, label=f'f0 = {f0} Hz')
axes[2].axvline(x=f1, color='purple', linestyle='--', alpha=0.5, label=f'f1 = {f1} Hz')
axes[2].legend()

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

import numpy as np
import matplotlib.pyplot as plt

# パラメータ設定
fc = 10.0
Rs = 1.0
Ts = 1.0 / Rs
h = 0.5            # MSKの変調指数
delta_f = h * Rs   # 周波数偏移 = 0.5 Hz
fs = 1000.0
N_bits = 12

# ランダムビット列
np.random.seed(123)
bits = np.random.randint(0, 2, N_bits)
symbols = 2 * bits - 1  # {-1, +1}

# MSK信号の生成（CPFSK方式）
t_total = N_bits * Ts
t = np.arange(0, t_total, 1/fs)
dt = 1.0 / fs

msk_signal = np.zeros(len(t))
phase = np.zeros(len(t))
cum_phase = 0.0

for i in range(len(t)):
    symbol_idx = min(int(t[i] / Ts), N_bits - 1)
    a = symbols[symbol_idx]
    # 瞬時周波数 = fc + a * delta_f/2 ではなく fc + a * delta_f
    # h = 0.5 なので delta_f = 0.5, 瞬時周波数偏移 = +/-0.25 Hz
    freq_dev = a * delta_f / 2
    msk_signal[i] = np.cos(2 * np.pi * fc * t[i] + cum_phase)
    phase[i] = cum_phase
    cum_phase += 2 * np.pi * freq_dev * dt

# 位相をアンラップして正規化
phase_normalized = phase / np.pi

fig, axes = plt.subplots(3, 1, figsize=(12, 10))

# ビット列
axes[0].step(np.arange(N_bits), bits, where='mid', linewidth=2, color='tab:blue')
axes[0].set_ylabel('Bit value')
axes[0].set_title('Data Bits')
axes[0].set_xlim(-0.5, N_bits - 0.5)
axes[0].set_ylim(-0.2, 1.4)
axes[0].grid(True, alpha=0.3)

# MSK波形
axes[1].plot(t, msk_signal, linewidth=0.8, color='tab:orange')
axes[1].set_ylabel('Amplitude')
axes[1].set_title(f'MSK Signal (h=0.5, fc={fc} Hz)')
axes[1].set_xlim(0, t_total)
axes[1].grid(True, alpha=0.3)

# シンボル境界
for i in range(1, N_bits):
    axes[1].axvline(x=i * Ts, color='gray', linestyle=':', alpha=0.5)

# 位相トラジェクトリ
axes[2].plot(t, phase_normalized, linewidth=2, color='tab:red')
axes[2].set_xlabel('Time [s]')
axes[2].set_ylabel('Phase [x π rad]')
axes[2].set_title('MSK Phase Trajectory')
axes[2].set_xlim(0, t_total)
axes[2].grid(True, alpha=0.3)

# シンボル境界とπ/2の倍数のグリッド
for i in range(1, N_bits):
    axes[2].axvline(x=i * Ts, color='gray', linestyle=':', alpha=0.5)

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

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erfc

def simulate_ook_ber(snr_db_range, n_bits=100000):
    """OOKのコヒーレント検波BERシミュレーション"""
    ber_list = []
    for snr_db in snr_db_range:
        snr_lin = 10 ** (snr_db / 10)
        # OOK: Eb = A^2*Ts/4 (平均), sigma^2 = N0/2
        A = 2.0
        Eb = A ** 2 / 4  # 平均ビットエネルギー（Ts=1と正規化）
        N0 = Eb / snr_lin
        sigma = np.sqrt(N0 / 2)

        bits = np.random.randint(0, 2, n_bits)
        # 整合フィルタ出力
        z = bits * A + sigma * np.random.randn(n_bits)
        # 最適閾値はA/2
        bits_hat = (z > A / 2).astype(int)
        ber = np.mean(bits != bits_hat)
        ber_list.append(max(ber, 1e-7))
    return np.array(ber_list)

def simulate_bfsk_coherent_ber(snr_db_range, n_bits=100000):
    """コヒーレントBFSKのBERシミュレーション"""
    ber_list = []
    for snr_db in snr_db_range:
        snr_lin = 10 ** (snr_db / 10)
        Eb = 1.0
        N0 = Eb / snr_lin
        sigma = np.sqrt(N0 / 2)

        bits = np.random.randint(0, 2, n_bits)
        # 2つの相関器出力（直交信号）
        # ビット=0: s0成分=sqrt(Eb), s1成分=0
        # ビット=1: s0成分=0, s1成分=sqrt(Eb)
        z0 = (1 - bits) * np.sqrt(Eb) + sigma * np.random.randn(n_bits)
        z1 = bits * np.sqrt(Eb) + sigma * np.random.randn(n_bits)
        # 判定: z1 > z0 なら 1
        bits_hat = (z1 > z0).astype(int)
        ber = np.mean(bits != bits_hat)
        ber_list.append(max(ber, 1e-7))
    return np.array(ber_list)

def simulate_bfsk_noncoherent_ber(snr_db_range, n_bits=100000):
    """非コヒーレントBFSKのBERシミュレーション"""
    ber_list = []
    for snr_db in snr_db_range:
        snr_lin = 10 ** (snr_db / 10)
        Eb = 1.0
        N0 = Eb / snr_lin
        sigma = np.sqrt(N0 / 2)

        bits = np.random.randint(0, 2, n_bits)
        # 包絡線検波: 各フィルタの I/Q 出力
        z0_I = (1 - bits) * np.sqrt(Eb) + sigma * np.random.randn(n_bits)
        z0_Q = sigma * np.random.randn(n_bits)
        z1_I = bits * np.sqrt(Eb) + sigma * np.random.randn(n_bits)
        z1_Q = sigma * np.random.randn(n_bits)
        # 包絡線
        env0 = np.sqrt(z0_I**2 + z0_Q**2)
        env1 = np.sqrt(z1_I**2 + z1_Q**2)
        bits_hat = (env1 > env0).astype(int)
        ber = np.mean(bits != bits_hat)
        ber_list.append(max(ber, 1e-7))
    return np.array(ber_list)

# シミュレーション実行
snr_db = np.arange(0, 16, 1)

ber_ook = simulate_ook_ber(snr_db, n_bits=200000)
ber_bfsk_coh = simulate_bfsk_coherent_ber(snr_db, n_bits=200000)
ber_bfsk_ncoh = simulate_bfsk_noncoherent_ber(snr_db, n_bits=200000)

# 理論曲線
snr_fine = np.linspace(0, 15, 200)
snr_lin_fine = 10 ** (snr_fine / 10)

# BPSK理論値
ber_bpsk_theory = 0.5 * erfc(np.sqrt(snr_lin_fine))
# OOK理論値
ber_ook_theory = 0.5 * erfc(np.sqrt(snr_lin_fine / 2))
# コヒーレントBFSK理論値
ber_bfsk_coh_theory = 0.5 * erfc(np.sqrt(snr_lin_fine / 2))
# 非コヒーレントBFSK理論値
ber_bfsk_ncoh_theory = 0.5 * np.exp(-snr_lin_fine / 2)

fig, ax = plt.subplots(figsize=(10, 7))

# シミュレーション結果
ax.semilogy(snr_db, ber_ook, 'rs', markersize=6, label='OOK (sim)')
ax.semilogy(snr_db, ber_bfsk_coh, 'b^', markersize=6, label='Coherent BFSK (sim)')
ax.semilogy(snr_db, ber_bfsk_ncoh, 'gD', markersize=5, label='Non-coherent BFSK (sim)')

# 理論曲線
ax.semilogy(snr_fine, ber_bpsk_theory, 'k-', linewidth=2, label='BPSK/MSK (theory)')
ax.semilogy(snr_fine, ber_ook_theory, 'r--', linewidth=1.5, label='OOK (theory)')
ax.semilogy(snr_fine, ber_bfsk_coh_theory, 'b--', linewidth=1.5, label='Coherent BFSK (theory)')
ax.semilogy(snr_fine, ber_bfsk_ncoh_theory, 'g--', linewidth=1.5, label='Non-coherent BFSK (theory)')

ax.set_xlabel('Eb/N0 [dB]')
ax.set_ylabel('BER')
ax.set_title('BER Comparison: ASK and FSK Modulation Schemes')
ax.legend(loc='lower left')
ax.grid(True, which='both', alpha=0.3)
ax.set_ylim(1e-6, 1)
ax.set_xlim(0, 15)

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

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erfc

def simulate_mfsk_ber(M, snr_db_range, n_symbols=100000):
    """M-FSKのコヒーレント検波BERシミュレーション"""
    m = int(np.log2(M))
    ber_list = []

    for snr_db in snr_db_range:
        snr_lin = 10 ** (snr_db / 10)
        Eb = 1.0
        Es = Eb * m  # シンボルエネルギー
        N0 = Eb / snr_lin
        sigma = np.sqrt(N0 / 2)

        # ランダムシンボル
        symbol_indices = np.random.randint(0, M, n_symbols)
        bit_errors = 0
        total_bits = 0

        # M個の相関器出力をベクトル化
        # 送信シンボルiの場合: z_i ~ N(sqrt(Es), sigma^2), z_j ~ N(0, sigma^2) (j != i)
        z = sigma * np.random.randn(n_symbols, M)
        for n in range(n_symbols):
            z[n, symbol_indices[n]] += np.sqrt(Es)

        # 最大値の相関器を選択
        detected = np.argmax(z, axis=1)

        # ビット誤り数をカウント（Gray符号は一般にM-FSKには適用しにくいので、
        # シンボル誤りあたりの平均ビット誤り数 = m*M/(2*(M-1)) を使う）
        symbol_errors = np.sum(detected != symbol_indices)
        # シンボル誤りをビット誤りに変換（M直交信号の場合）
        ber = symbol_errors / n_symbols * (M / 2) / (M - 1)
        ber_list.append(max(ber, 1e-7))

    return np.array(ber_list)

snr_db = np.arange(0, 18, 1)

fig, ax = plt.subplots(figsize=(10, 7))

for M in [2, 4, 8, 16]:
    ber = simulate_mfsk_ber(M, snr_db, n_symbols=200000)
    ax.semilogy(snr_db, ber, 'o-', markersize=5, label=f'{M}-FSK (sim)')

# BPSK理論値
snr_fine = np.linspace(0, 17, 200)
snr_lin_fine = 10 ** (snr_fine / 10)
ber_bpsk = 0.5 * erfc(np.sqrt(snr_lin_fine))
ax.semilogy(snr_fine, ber_bpsk, 'k--', linewidth=2, label='BPSK (theory)')

# コヒーレントBFSK理論値
ber_bfsk_theory = 0.5 * erfc(np.sqrt(snr_lin_fine / 2))
ax.semilogy(snr_fine, ber_bfsk_theory, 'b--', linewidth=1.5, label='BFSK (theory)')

ax.set_xlabel('Eb/N0 [dB]')
ax.set_ylabel('BER')
ax.set_title('BER of M-FSK (Coherent Detection)')
ax.legend()
ax.grid(True, which='both', alpha=0.3)
ax.set_ylim(1e-6, 1)
ax.set_xlim(0, 17)

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

import numpy as np
import matplotlib.pyplot as plt

# パラメータ
fc = 50.0       # 搬送波周波数
Rb = 10.0       # ビットレート
Ts = 1.0 / Rb   # ビット周期
fs = 10000.0    # サンプリング周波数
N_bits = 1000   # ビット数

np.random.seed(42)
bits = np.random.randint(0, 2, N_bits)
symbols_bipolar = 2 * bits - 1

t_total = N_bits * Ts
t = np.arange(0, t_total, 1/fs)

# OOK信号
ook = np.zeros(len(t))
for i in range(N_bits):
    mask = (t >= i * Ts) & (t < (i + 1) * Ts)
    ook[mask] = bits[i] * np.cos(2 * np.pi * fc * t[mask])

# BPSK信号（参考）
bpsk = np.zeros(len(t))
for i in range(N_bits):
    mask = (t >= i * Ts) & (t < (i + 1) * Ts)
    bpsk[mask] = symbols_bipolar[i] * np.cos(2 * np.pi * fc * t[mask])

# BFSK信号（h=1）
delta_f_bfsk = Rb  # h = delta_f * Ts = 1
f0_bfsk = fc - delta_f_bfsk / 2
f1_bfsk = fc + delta_f_bfsk / 2
bfsk = np.zeros(len(t))
for i in range(N_bits):
    mask = (t >= i * Ts) & (t < (i + 1) * Ts)
    freq = f1_bfsk if bits[i] == 1 else f0_bfsk
    bfsk[mask] = np.cos(2 * np.pi * freq * t[mask])

# MSK信号（h=0.5, CPFSK）
delta_f_msk = 0.5 * Rb
msk = np.zeros(len(t))
phase_msk = 0.0
dt = 1.0 / fs
for i in range(len(t)):
    sym_idx = min(int(t[i] / Ts), N_bits - 1)
    freq_dev = symbols_bipolar[sym_idx] * delta_f_msk / 2
    msk[i] = np.cos(2 * np.pi * fc * t[i] + phase_msk)
    phase_msk += 2 * np.pi * freq_dev * dt

# スペクトル計算
from numpy.fft import fft, fftfreq
N_fft = len(t)
freq_axis = fftfreq(N_fft, 1/fs)

signals = {'OOK': ook, 'BPSK': bpsk, 'BFSK (h=1)': bfsk, 'MSK (h=0.5)': msk}
colors = {'OOK': 'tab:red', 'BPSK': 'black', 'BFSK (h=1)': 'tab:blue', 'MSK (h=0.5)': 'tab:green'}

fig, ax = plt.subplots(figsize=(12, 7))
mask_freq = (freq_axis >= 30) & (freq_axis <= 70)

for name, sig in signals.items():
    spec = np.abs(fft(sig)) ** 2 / N_fft
    # 移動平均でスムージング
    window = 50
    spec_smooth = np.convolve(spec, np.ones(window)/window, mode='same')
    psd_db = 10 * np.log10(spec_smooth[mask_freq] / np.max(spec_smooth[mask_freq]) + 1e-15)
    ax.plot(freq_axis[mask_freq], psd_db, linewidth=1.5, color=colors[name],
            label=name, alpha=0.8)

ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Normalized PSD [dB]')
ax.set_title(f'Spectral Comparison of Modulation Schemes (fc={fc} Hz, Rb={Rb} bit/s)')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_ylim(-40, 5)
ax.axvline(x=fc, color='gray', linestyle=':', alpha=0.5)

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